# VISAGE: Video Instance Segmentation with Appearance-Guided Enhancement

Hanjung Kim1 Jaehyun Kang1 Miran Heo1 Sukjun Hwang2
Seoung Wug Oh3 Seon Joo Kim1
1Yonsei University 2Carnegie Mellon University 3Adobe Research

###### Abstract

In recent years, online Video Instance Segmentation (VIS) methods have shown
remarkable advancement with their powerful query-based detectors. Utilizing
the output queries of the detector at the frame level, these methods achieve
high accuracy on challenging benchmarks. However, we observe the heavy
reliance of these methods on the location information that leads to incorrect
matching when positional cues are insufficient for resolving ambiguities.
Addressing this issue, we present VISAGE that enhances instance association by
explicitly leveraging appearance information. Our method involves a generation
of queries that embed appearances from backbone feature maps, which in turn
get used in our suggested simple tracker for robust associations. Finally,
enabling accurate matching in complex scenarios by resolving the issue of
over-reliance on location information, we achieve competitive performance on
multiple VIS benchmarks. For instance, on YTVIS19 and YTVIS21, our method
achieves 54.5 AP and 50.8 AP. Furthermore, to highlight appearance-awareness
not fully addressed by existing benchmarks, we generate a synthetic dataset
where our method outperforms others significantly by leveraging the appearance
cue. Code will be made available at https://github.com/KimHanjung/VISAGE.

Figure 1:  Qualitative results across challenging scenarios. Predicted results
using query-propagation [16, 37], query-matching [17, 36, 29], and our
appearance-guided methods. The first row illustrates a shot change across
consecutive frames, a scenario where previous methods fail to maintain
consistent tracking. The second and third rows demonstrate trajectory
intersections, leading to id-switching with previous methods. Unlike previous
methods, our method successfully tracks objects without switching or losses.
Best viewed in color.

##  1 Introduction

Video Instance Segmentation (VIS) is a challenging task that requires
classification, segmentation, and tracking of distinct instances throughout a
video sequence [31]. Current studies in VIS can be primarily categorized into
two approaches: online and offline, based on whether a video is processed in a
per-frame or per-clip manner. Recently, the advancement of frame-level object
detectors has resulted in online methods becoming increasingly dominant in the
VIS field.

As detectors directly impact the accuracy in the video domain, recent online
models are primarily built using the powerful query-based detectors [38, 5].
Spatially decoding image information, the detectors are designed to represent
object-wise information using the queries. Therefore, the online VIS methods
reuse these queries from the detectors and have achieved substantial
improvements in multiple challenging benchmarks [25, 33] by mostly adopting
either propagation [16, 23] or matching [17, 29, 36] strategies. However,
tracking under complex scenarios such as shot changes or trajectory
intersections (Fig. 1) are yet imperfect, resulting in the degradation of the
overall accuracy.

Examining these failure cases, we observe that object-wise information of the
queries are significantly imbalanced: heavy reliance on positional cues, and
less reflection of appearances. As demonstrated in Fig. 1, previous query-
based methods [16, 36, 17, 29, 37], display a tendency to maintain the spatial
order of previous predictions in their current predictions. To support this
argument, we conduct additional analyses by horizontally flipping images to
generate two-frame pseudo videos. The existing models manifest association
errors despite the distinct exterior patterns of objects (Fig. 2 Top), which
highlight the dependence on object locations. As there exist multiple
scenarios that cannot be fully handled with the imbalanced information, such a
phenomenon necessitates the models to take object appearances into
consideration.

We introduce VISAGE (Video Instance Segmentation with Appearance-Guided
Enhancement), a method that leverages appearance cues as a crucial indicator
for distinguishing instances. First, we present a lightweight branch for
generating appearance queries that encapsulate the visual characteristics of
an instance. Specifically, given mask predictions of object queries from the
detector, this branch explicitly embeds appearance by performing mask pooling
operations to backbone feature maps. Then, our method separately uses the
contrastive loss [26, 3] to pairs of appearance queries and object queries,
which improve the association of instances across a video by alleviating the
position-appearance imbalance.

Additionally, we present a simple tracker that employs appearance cues as an
additional indicator and leverages temporal information. Previous methods [29,
36] employ multiple refinement steps for removing redundant mask proposals and
adopt handcrafted threshold values to improve the accuracy. While these
intricate tracker modules raise the tracking quality, they become dependent on
multiple hyperparameters that are heuristically optimized for each dataset. In
order to alleviate such dependence, our method streamlines the tracker and
dramatically reduces the number of hyperparameters, eliminating thresholding
and instead incorporating appearance cues. Nonetheless, a lack of temporal
information still exists, inherent in query-matching online methods, as they
are only aware of adjacent frames. We address this limitation by using a
simple memory bank to facilitate temporal awareness.

Our approach VISAGE has many desirable properties. Enhanced by the appearance
guidance, VISAGE achieves robust tracking in complex scenarios, as illustrated
at Fig. 1 (a) and bottom of Fig. 2. Furthermore, with its simplified tracker
that effectively utilizes the past history of both object and appearance
queries, VISAGE demonstrates competitive performance across all benchmark
datasets. We validate our methodology through extensive evaluations on the
four standard VIS benchmarks: YouTubeVIS-19/21/22[31, 32, 33] and OVIS [25].
Unfortunately, these datasets are mostly composed of videos with general
scenarios, thus not suitable for precise analyses of the effectiveness of
leveraging appearance embeddings. Therefore, to validate our method’s
effectiveness in leveraging appearance cues, we develop a large-scale pseudo
dataset and conduct comparisons of its performance with other methods.
Finally, VISAGE demonstrates superior performance over previous methods [16,
36] in the comparison.

Figure 2:  Proof of concept demonstrated with a flipped image. Previous
methods [17, 36, 37, 16, 29] struggle with instance matching in flipped
images, showing a dependency on location. Our method, VISAGE, addresses this
by emphasizing appearance, enabling accurate instance matching even with image
flipping.

##  2 Related Works

###  2.1 Online Video Instance Segmentation

Contrary to traditional online approaches [31, 34, 11, 19], modern online
methods utilize query-based detectors with an emphasis on query association
strategy. Query-based online methods can be divided into two main strategies:
query-matching and query-propagation.

Query-matching approaches focus on dynamically building tracklets in an online
manner. For instance, MinVIS [17] utilizes Mask2Former [5] as its query-based
detector trained exclusively at the image level. During inference, MinVIS
independently predicts outputs at the frame level and then associates them
through bipartite matching between frame queries. To enhance the
distinctiveness of queries, IDOL [29] introduces contrastive learning,
inspired by QDTrack [8]. They also adopt same inference process with QDTrack,
additionally propose temporal weighted softmax with memory bank for leveraging
temporal information. However, those memory bank is not used at training
stage. CTVIS [36] further improves upon IDOL by removing the discrepancy
between training and inference caused by the use of the memory bank. By
utilizing the memory bank during training, CTVIS eliminates this discrepancy
and incorporates historical instances. All these methods rely on the outputs
of query-based detectors for their matching processes. In contrast, our
approach VISAGE, introduces additional query to ensure more stable matching.

Figure 3:  Overview of VISAGE. (a) The proposed VISAGE’s architecture which
generate object embedding and appearance embedding. (b) Overall inference
pipeline of VISAGE: At time step $t-1$, the memory bank is updated with both
the appearance embedding and the object embedding. Then, at time step $t$, the
memory embedding is read from the memory bank and used for matching. (c).
Details of the matching process: In that scenario, using only object
embeddings leads to incorrect matching. However, when guided by the appearance
embedding, the matching process can be corrected. Best viewed in color.

Query-propagation methods utilize output queries from previous frames to
detect objects in the current frame. InsPro [12] propagates their instance
query and proposal for association and proposes duplicate removal strategy.
CAROQ [6] extends Mask2Former [5] for application in VIS methods and
propagates its output queries along with relative object queries. GenVIS [16]
employs an offline VIS model [15], and adapts it to a query-propagation
approach, notably omitting heuristic association strategies. This method
reduces computational overhead by freezing the backbone network and query-
based detector, thus effectively handling longer video sequences. Similarly,
DVIS [37] introduces referring cross-attention, which is designed to capture
the correlation between the current and previous frames using their queries.
Yet, these methods are limited by their local matching strategy, which focuses
on aligning tracklets with ground truths. This approach often leads to
unstable tracking outcomes, as it does not consider the full video context.
TCOVIS [23] addresses this issue by shifting from the local matching strategy
of GenVIS to a global matching approach, thereby achieving more robust
tracking across videos.

###  2.2 Offline Video Instance Segmentation

Previous offline VIS methods [1, 30] treat video clip input as a spatio-
temporal volume. However, with the advent of query-based detectors, recent
offline methods have shifted to treating videos through queries. Recent
offline methods mostly leverage query-based detectors. VisTR [27] builds upon
DETR [2] to create an end-to-end framework. However, the self-attention
mechanism in VisTR is computationally heavy. To address this, IFC [18]
introduces memory tokens that allow communication between frames, reducing the
computational overhead. SeqFormer [28] restructures the decoder to process
each frame separately, enhancing the model’s sensitivity to instance movement
by capturing each location. Mask2Former-VIS [4] extends the Mask2Former [5]
for video-level tasks, achieving notable performance improvements. TeViT [35]
presents a novel backbone that improves temporal information processing,
leveraging the capabilities of ViT [7]. VITA [15] effectively aggregates video
information by learning the relationships between the frame-level queries.

##  3 Method

We propose a novel online VIS framework that leverages appearance information
to reduce the heavy reliance on location information. Our newly introduced
appearance-aware query guides the association procedure, ensuring correct
matching. Finally, we introduce a simple inference method that aligns with our
design principles, without heuristics.

###  3.1 Query-based Detector

In most recent query-based object detection frameworks, the architecture is
typically segmented into three integral components: the backbone, the
transformer encoder, and the transformer decoder. Initially, the backbone
generates the low-level image feature maps, which encapsulate the essential
visual information. Subsequent to this, the transformer encoder employs self-
attention mechanisms to enhance the feature representation, as detailed in [2,
38]. Finally, in the transformer decoder, the identified objects are decoded
into $N$ learnable queries.

Our method uses the well-established query-based detection framework [5],
maintaining the original structure with the backbone, transformer encoder, and
transformer decoder. Furthermore, we utilize the feature maps from the
backbone to create the appearance query, and we use the output of the
transformer decoder as the object query.

###  3.2 Appearance-Guided Enhancement

Before the emergence of query-based tracking approaches, video tracking
methods [31, 8, 20] generate instance features from backbone feature maps.
This is achieved using operations like RoIPool [10] and RoIAlign [14],
followed by associating them across the video. In line with these traditional
methods, we use average pooling to extract appearance feature from the
backbone feature maps, guided by the predicted masks as shown in Fig. 3 (a).
Therefore, in addition to the object queries already in use, we further employ
the appearance query as an additional indicator. We transform both appearance
queries and object queries, leading to the creation of appearance embeddings
$\mathbf{e}_{a}\in\mathbb{R}^{N\times C}$ and object embeddings
$\mathbf{e}_{i}\in\mathbb{R}^{N\times C}$, respectively.

Our appearance embeddings enhance the matching process when using only object
embeddings leads to incorrect matches, as illustrated in Fig. 3 (c). Relying
solely on the similarity of object embeddings may lead to ambiguity due to
their similar positions across subsequent frames. However, this ambiguity can
be resolved by also considering the similarity of appearance embeddings, as
their distinct appearances provide additional discriminative information. With
this guidance, we leverage both appearance embeddings $\mathbf{e}_{a}$ and
object embeddings $\mathbf{e}_{i}$ to identify the optimal match between
queries across subsequent frames as shown in Algorithm 1 (lines 4-10).

To ensure more distinctive embeddings for better association, we apply
contrastive loss to both object and appearance embeddings, similar to the
approach used in previous methods [29, 8]. We sample key frame and reference
frame, then calculate the loss for each set of embeddings separately.
Consequently, our final loss is a weighted summation of original query-based
detector’s loss and the contrastive loss for both the appearance embeddings
and object embeddings.

Algorithm 1 Inference pipeline of VISAGE.

1:model $\mathcal{F}$, frames $\\{x^{t}\\}_{t=1}^{T}$

2:predictions $\mathcal{P}$

3:$\mathcal{P}\leftarrow\\{\\}$

4:$\mathcal{M}\leftarrow\\{\\}$

5:$\texttt{idx}\leftarrow[0,N-1]$

6:for $t$ in $[1,T]$:

7:
${\mathbf{p}}^{t},\mathbf{e}_{i}^{t},\mathbf{e}_{a}^{t}\leftarrow\mathcal{F}(x^{t})$

8:
$\mathcal{M}\leftarrow\mathcal{M}+(\mathbf{e}_{i}^{t}[\texttt{idx}],\mathbf{e}_{a}^{t}[\texttt{idx}])$

9:
$\mathbf{m}_{i}^{t},\mathbf{m}_{a}^{t}\leftarrow\mathcal{M}\texttt{.read_memory()}$

10: $\mathcal{P}\leftarrow\mathcal{P}+{\mathbf{p}}^{t}[\texttt{idx}]$

11:
$\mathbf{s}\leftarrow\texttt{cos_sim}(\mathbf{e}_{i}^{t},\mathbf{m}_{i}^{t})+\texttt{cos_sim}(\mathbf{e}_{a}^{t},\mathbf{m}_{a}^{t})$

12: $\texttt{idx}\leftarrow\texttt{linear_sum_assignment}(\mathbf{s})$

13:end for

14:return $\mathcal{P}$

###  3.3 Inference with Appearance

As shown in Algorithm 1, our tracker employs a simple yet effective matching
process that aligns the current object embeddings and appearance embeddings
with their respective counterparts from the previous frame. We compute the
similarity scores for each object-appearance embedding pair using Cosine
Similarity, following a method similar to that in [17]. Then, employing the
sum of object and appearance similarities, we utilize the Hungarian algorithm
[21] (linear_sum_assignment in Algorithm 1) to achieve optimal assignment.

In addition, we incorporate a simple memory bank to compensate for the lack of
temporal information, a limitation inherent in online processes, as used in
previous methods [29, 36]. We stack the states of previous queries from the
most recent frames within our memory bank, which has a size of $W$. From this
bank, we read a memory embedding $m\in\mathbb{R}^{N\times C}$ using the
read_memory() function in Algorithm 1. Specifically, this function is
implemented by temporally weighting the embeddings to put more emphasis on
recent queries while utilizing the confidence scores for selective weighting.
For the current memory embedding $m^{t}$, the calculation of the weighting for
each embedding at the previous time step $w\in[1,W]$ can be formally expressed
as follows:

$\mathbf{m}^{t}=\sum\limits_{w=1}^{W}\left(\mathbf{e}^{t-w}s^{t-w}\times\frac{W}{w}\right),$
(1)

where $s$ denotes the confidence score. This memory embedding represents the
object and appearance, denoted as $m^{t}_{i}$ and $m^{t}_{a}$ respectively, in
Algorithm 1.

Finally, we sequentially associate the predictions ($\mathbf{p}$ in Algorithm
1) from each frame using the obtained assignment. By employing a simple
inference pipeline coupled with an efficient memory bank, we introduce an
effective approach.

##  4 Experiments

Method | Setting | YouTube-VIS 2019 | YouTube-VIS 2021
---|---|---|---
AP | AP50 | AP75 | AR1 | AR10 | AP | AP50 | AP75 | AR1 | AR10
IFC [18] | offline | 41.2 | 65.1 | 44.6 | 42.3 | 49.6 | 35.2 | 55.9 | 37.7 | 32.6 | 42.9
Mask2Former-VIS [4] | offline | 46.4 | 68.0 | 50.0 | - | - | 40.6 | 60.9 | 41.8 | - | -
MinVIS [17] | online | 47.4 | 69.0 | 52.1 | 45.7 | 55.7 | 44.2 | 66.0 | 48.1 | 39.2 | 51.7
IDOL [29] | online | 49.5 | 74.0 | 52.9 | 47.7 | 58.7 | 43.9 | 68.0 | 49.6 | 38.0 | 50.9
VITA [15] | offline | 49.8 | 72.6 | 54.5 | 49.4 | 61.0 | 45.7 | 67.4 | 49.5 | 40.9 | 53.6
GenVIS [16] | online | 50.0 | 71.5 | 54.6 | 49.5 | 59.7 | 47.1 | 67.5 | 51.5 | 41.6 | 54.7
GenVIS [16] | offline | 51.3 | 72.0 | 57.8 | 49.5 | 60.0 | 46.3 | 67.0 | 50.2 | 40.6 | 53.2
DVIS [37] | online | 51.2 | 73.8 | 57.1 | 47.2 | 59.3 | 46.4 | 68.4 | 49.6 | 39.7 | 53.5
DVIS [37] | offline | 52.6 | 76.5 | 58.2 | 47.4 | 60.4 | 47.4 | 71.0 | 51.6 | 39.9 | 55.2
TCOVIS [23] | online | 52.3 | 73.5 | 57.6 | 49.8 | 60.2 | 49.5 | 71.2 | 53.8 | 41.3 | 55.9
CTVIS [36] | online | 55.1 | 78.2 | 59.1 | 51.9 | 63.2 | 50.1 | 73.7 | 54.7 | 41.8 | 59.5
VISAGE | online | 54.5 | 77.2 | 60.3 | 50.9 | 61.7 | 50.8 | 72.9 | 55.6 | 43.3 | 58.8
Table 1:  Comparisons on the YouTube-VIS 2019 and 2021 validation sets.
Methods are denoted as online or offline, indicated by the text color. We
highlight the best performance in bold and the second-best performance with an
underline.

###  4.1 Datasets

We evaluate VISAGE on four VIS benchmarks: Youtube-VIS (YTVIS) 2019/2021/2022
[31, 32, 33] and Occluded VIS (OVIS) [25]. The YTVIS datasets contain 40
predefined categories in their videos. YTVIS 2019 is the first and largest
dataset for video instance segmentation. It includes 2,238 videos for
training, 302 for validation, and 343 for the testing. YTVIS 2021 expands upon
the YTVIS 2019 dataset, containing 2,985 videos for training, 421 for
validation, and 453 for testing, while refining annotations and modifying some
categories. Additionally, to represent challenging scenarios with longer
sequences, YTVIS 2022 is added. It adds 71 extra videos to the validation set
of YTVIS 2021.

OVIS stands out for its more complex and longer videos, featuring 25
categories and comprising 2,238 videos for training, 302 for validation, and
343 for testing. In comparison to YTVIS, it contains a higher number of
instances per video, averaging 5.8, and a total of 296k masks. Additionally,
the average length of videos in OVIS is approximately 12 seconds.

###  4.2 Implementation Details

We adopt the Mask2Former [5] as our query-based detector. All of our models
are initialized with parameters pre-trained on COCO dataset [24] with
ResNet-50 [13] backbone. We also adopt COCO joint training following previous
works [28, 15, 16, 36, 23]. Our batch includes 16 videos. In our experimental
setup, we set the weights for the contrastive losses applied to both the
appearance and object embeddings at $2.0$. We follow our detector’s loss
weight configurations as detailed in [5] for other losses in our setup. The
window size of memory bank $W$ is set to 5 during inference. Finally, our
method is trained using 4 NVIDIA A6000 GPUs.

Method | AP | AP50 | AP75 | AR1 | AR10
---|---|---|---|---|---
MinVIS [17] | 23.3 | 47.9 | 19.3 | 20.2 | 28.0
VITA [15] | 32.6 | 53.9 | 39.3 | 30.3 | 42.6
GenVIS [16] | 37.5 | 61.6 | 41.5 | 32.6 | 42.2
TCOVIS [23] | 38.6 | 59.4 | 41.6 | 32.8 | 46.7
VISAGE | 36.7 | 60.0 | 41.1 | 35.4 | 45.1

Table 2:  Comparison on the YTVIS 2022 long videos set. Bold and underline
highlight the highest and second-highest performances, respectively. For
methods that are both online and semi-online, the online version is used.

###  4.3 Main Results

We compare our methods with state-of-the-art methodsg on four VIS benchmarks:
YTVIS 2019/2021/2022 and OVIS. The results are reported in Tab. 1, Tab. 2 and
Tab. 3.

#### Youtube-VIS 2019 & 2021.

As shown in Tab. 1, we compare our VISAGE with state-of-the-art methods. On
the YTVIS 2019 benchmark, VISAGE demonstrates performance comparable to the
highest-performing method [36]. Notably, both our VISAGE method and the best-
performing method [36], surpass previous methods by a large margin. In the
YTVIS 2021 benchmark, an improved version of YTVIS 2019, VISAGE outperforms
other existing methods. We achieve competitive performance on both the YTVIS
2019 and 2021 datasets, by incorporating appearance-guided enhancement.

#### Youtube-VIS 2022.

The results for the YTVIS 2022 benchmark are presented in Tab. 2. This
benchmark contains longer videos compared to previous benchmarks (YTVIS
2019/21). Our VISAGE exhibits slightly lower performance compared to GenVIS
[16] and TCOVIS [23]. These methods utilize the entire history of previous
information by propagating the queries. Nevertheless, our method employs only
a fixed number of previous frames and shows comparable performance.

Method | AP | AP50 | AP75 | AR1 | AR10
---|---|---|---|---|---
VITA [15] | 19.6 | 41.2 | 17.4 | 11.7 | 26.0
MinVIS [17] | 25.0 | 45.5 | 24.0 | 13.9 | 29.7
IDOL [29] | 30.2 | 51.3 | 30.0 | 15.0 | 37.5
GenVIS [16] | 35.8 | 60.8 | 36.2 | 16.3 | 39.6
DVIS [37] | 31.0 | 54.8 | 31.9 | 15.2 | 37.6
TCOVIS [23] | 35.3 | 60.7 | 36.6 | 15.7 | 39.5
CTVIS [36] | 35.5 | 60.8 | 34.9 | 16.1 | 41.9
VISAGE | 34.8 | 59.8 | 35.3 | 16.2 | 39.4

Table 3:  Comparison on the OVIS validation set. Bold and underline indicate
the best and second-best performances, respectively. We use the online version
for GenVIS [16] and DVIS [37].

Feature | AP | AP50 | AP75 | AR1 | AR10
---|---|---|---|---|---
Transformer Encoder | 52.8 | 73.8 | 58.1 | 49.4 | 60.2
Backbone | 54.5 | 77.2 | 60.3 | 50.9 | 61.7

Table 4:  Ablation study of the target feature maps which generate appearance
query

#### OVIS.

Our method, VISAGE, is evaluated on the OVIS dataset. In Tab. 3, VISAGE
demonstrates competitive performance when compared to other state-of-the-art
methods [16, 23, 36]. It is important to note that these methods incorporate
sequential information during training. Consequently, these methods shows
marginally better performance than VISAGE, which is expected given that the
OVIS dataset is characterized by its longer video sequences.

###  4.4 Ablation Studies

In this section, we provide the ablation studies for our proposed method and
discuss its effect. All ablation experiments are conducted on YTVIS 2019 [31]
validation set.

#### Appearance feature.

In Tab. 4, we analyze the impact of feature maps which create the appearance
query. When we make the appearance query from the transformer encoder feature
maps, there is a slight drop in performance. This indicates that backbone
feature maps contain rich visual information compare to transformer encoder
feature maps.

#### Appearance guidance.

Tab. 5 shows effectiveness of our appearance-guided enhancement. We omit
appearance information by not computing the appearance similarity in line 9 of
Algorithm 1, resulting in a matching process that relies solely on object
similarity. As indicated by rows 1 and 3 in Tab. 5, the absence of appearance
information leads to a degradation in performance. Further analysis of the
effectiveness of appearance guidance is discussed in Sec. 4.5.

#### Memory bank.

In Tab. 5, we demonstrate the necessity of the memory bank for compensating
for the lack of temporal information. We evaluate performance differences
between scenarios with and without the use of the memory bank. Rows 1 and 2 in
Tab. 5 show decreased performance compared to rows 3 and 4, respectively. If
there is no memory bank, our method lacks access to historical information
beyond the frame right before the current one, since it measures similarities
solely between the current frame and the one immediately preceding it. By
implementing a memory bank, our method gains awareness of previous frames,
which leads to improved performance.

Appearance |  Memory Bank | YouTube-VIS 2019
---|---|---
AP | AP50 | AP75 | AR1 | AR10
|  | 49.9 | 71.4 | 54.7 | 47.0 | 58.7
✓ |  | 50.5 | 71.6 | 55.5 | 47.4 | 60.6
| ✓ | 53.4 | 76.8 | 58.7 | 49.8 | 61.2
✓ | ✓ | 54.5 | 77.2 | 60.3 | 50.9 | 61.7

Table 5:  Ablation study on appearance guidance and memory bank utilization in
YTVIS 2019 validation Set, with memory bank window size $W=5$.

Tracker | AP | AP50 | AP75 | AR1 | AR10
---|---|---|---|---|---
[29, 36] (w/o NMS) | 52.7 | 75.3 | 57.6 | 49.3 | 61.6
[29, 36] (w/ NMS) | 54.4 | 77.9 | 59.0 | 49.9 | 62.8
ours (w/o NMS) | 54.5 | 77.2 | 60.3 | 50.9 | 61.7

Table 6:  Ablation study of the tracker type on YTVIS 2019 validation set.

#### Simplified tracker.

Previous query-matching based methods [29, 36], extend the tracker from [8],
which employs heuristic post-processing techniques such as NMS. In contrast,
our tracker relies solely on a cosine similarity metric and the Hungarian
algorithm for matching. To analyze the impact of such handcrafted design, we
modify our inference pipeline to align with the approaches used in previous
methods [29, 36]. As demonstrated in Tab. 6, both the previous strategies and
our simplified version show comparable performance. However, removing the
heuristic design element NMS from the previous approach results in performance
degradation. This outcome implies that previous approaches rely on handcrafted
design, unlike our simpler tracker.

###  4.5 Analysis of Appearance-Guided Enhancement

#### T-SNE Visualization.

As illustrated in Fig. 4, we analyze the impact of appearance-guided
enhancement in VISAGE on its association capabilities using the OVIS dataset
by visualizing the query embeddings for an entire video, with consistent
colors indicating the same instance. The visualization includes three types of
query embeddings: object-only, appearance-only, and a combination of both
object and appearance queries. The embeddings that combine object and
appearance query result in a well-clustered visualization. While the
appearance-only query embeddings exhibit reasonable clustering, the
integration of both types of queries leads to more distinct and pronounced
clustering. This enhanced clustering indicates that our appearance-guided
enhancement leads to better association.

Figure 4:  T-SNE visualization on the OVIS dataset. Each row representing
three different videos. Each column corresponds to the type of query embedding
utilized. Points plotted in the same color indicate the same instance across
the dataset. Best viewed in color.

#### Pseudo Dataset.

Traditional datasets do not fully cover the complex scenarios that our
approach is designed to address. To validate VISAGE in a more intuitive
manner, we construct a pseudo dataset consisting of synthetic videos.
Instances in pseudo videos are created by compositing objects from the COCO
dataset [24] with Copy-and-Paste augmentation [9]. The background images for
the pseudo dataset are randomly sourced from BG-20k [22]. Additionally, we
employ Bezier Curves to simulate the movement of objects. It includes two
types of videos: track, where objects move following arbitrary bezier curves,
and swap, where objects’ locations are randomly swapped along their
trajectories. Except for the movement of instances, both types of datasets are
generated under the same conditions.

As shown at the top of Fig. 5, instances in the track type pseudo video move
along a Bezier curve. Consequently, complex scenarios, such as intersections
between instances or movements out of the frame, occur naturally. On the other
hand, the swap type presents a challenge for methods that primarily rely on
location information, as this dependence results in incorrect matching. This
is illustrated at the bottom of Fig. 5, where the positions of instances in
the pseudo video are suddenly swapped. Such scenarios verify the method’s
awareness of appearance cues when location information is no longer a reliable
indicator.

Figure 5:  Visualization of the pseudo dataset. In track type videos,
instances move along random bezier curves. On the other hand, the swap type
refers to a scenario where the positions of each instance are exchanged in the
middle of the video. The colored dot above each instance represents the
corresponding instance in the swapped frame.

Method | AP | AP50 | AP75 | AR1 | AR10
---|---|---|---|---|---
Track | GenVIS [16] | 54.8 | 70.6 | 58.8 | 61.2 | 65.9
CTVIS [36] | 63.7 | 78.7 | 69.7 | 70.3 | 74.8
VISAGE (w/o app) | 65.2 | 80.9 | 71.3 | 70.2 | 75.2
VISAGE | 65.8 | 81.1 | 71.7 | 70.8 | 76.3
Swap | GenVIS [16] | 41.8 | 59.7 | 44.0 | 50.8 | 57.0
CTVIS [36] | 53.1 | 71.7 | 56.9 | 61.7 | 65.7
VISAGE (w/o app) | 51.8 | 72.1 | 54.4 | 57.7 | 62.8
VISAGE | 63.0 | 79.5 | 69.3 | 67.8 | 73.5

Table 7:  Comparisons on Pseudo dataset. Bold denote the highest accuracy.

In Tab. 7, we conduct an evaluation of various online VIS methods [16, 36]
including VISAGE on the pseudo dataset using published ResNet-50 backbone
weights, which are trained on the YouTube-VIS 2019 dataset. All of these
methods use COCO joint training, ensuring a fair comparison. The overall trend
is similar to that observed on the YTVIS 2019 dataset. In contrast, VISAGE
outperforms other methods in swap type videos, mainly due to the challenges
posed by the dataset in associating instances using only location information.
Consequently, other methods demonstrate degraded performance on swap type
videos compared to their performance on track types. Furthermore, the absence
of appearance information in the matching process notably impacts VISAGE’s
performance on swap type videos, although it remains relatively unaffected for
track type videos. It proves that appearance cue serves as an additional
indicator, providing robustness in scenarios with swapped object positions.

Figure 6:  Qualitative results of VISAGE. Videos are sourced from OVIS [25]
and YTVIS 2021 [32] datasets. These videos represent complex scenarios,
characterized by intersections and reappearances. Best viewed in color.

###  4.6 Qualitative Results

Fig. 6 displays the visualization results of our VISAGE across various
challenging scenarios. In the first row, several cats are positioned closely
together and the red-colored cat crosses another cat. In such cases, relying
solely on location information can often result in identity-switching errors.
However, our VISAGE tracks each cat accurately despite their adjacent
positions and sequential intersections. In the second row, we present a
similar scenario involving different classes: a person and a cow.

In the third row, we showcase another challenging scenario where an object
completely disappears and then reappears: a cat disappears and then reappears.
Our memory bank, which stores the previous queries of the cat, enables us to
reestablish the association when the cat reappears. The fourth row
demonstrates a similar scenario, where our memory bank proves effective: a
small sedan is completely occluded by the truck. Even though the large truck
in the foreground and the small, disappeared car are located closely, their
distinct appearances ensure accurate matching.

##  5 Limitations and Future Works

VISAGE addresses the insufficiency of existing query-based online VIS methods,
specifically their lack of appearance awareness. However, our approach has two
limitations. First, it assumes equal importance for appearance and location
information for its simplicity. In some scenarios, appearance may be more
critical than location, while in others, location may be prioritized.
Designing a method to adaptively weight the importance of appearance and
location information for stable tracking will be a future direction. Second,
VISAGE focuses exclusively on the query-matching approach, and we observe that
query-propagation methods also tend to rely on location information. Another
future direction involves incorporating appearance information awareness into
query-propagation approaches.

##  6 Conclusion

In this paper, we explore the importance of appearance in tracking objects, an
aspect often taken for granted yet overlooked by current online VIS methods.
We simply leverage the appearance cue when tracking objects, enabling more
effective distinction of objects in challenging scenarios. Our method, while
simple yet effective, achieves performance comparable to previous methods
across various VIS benchmarks. However, as existing VIS benchmarks do not
focus on appearance-requiring scenarios, we generate a synthetic dataset that
necessitates the use of appearance information. This dataset serves to
validate our appearance-aware approach and our approach surpasses other
methods with very large margin. We believe that recognizing and leveraging the
importance of appearance can lead to progress of VIS.

## References

  * Athar et al. [2020] Ali Athar, Sabarinath Mahadevan, Aljoša Ošep, Laura Leal-Taixé, and Bastian Leibe.  Stem-seg: Spatio-temporal embeddings for instance segmentation in videos.  In _ECCV_ , 2020.
  * Carion et al. [2020] Nicolas Carion, Francisco Massa, Gabriel Synnaeve, Nicolas Usunier, Alexander Kirillov, and Sergey Zagoruyko.  End-to-end object detection with transformers.  In _ECCV_ , 2020.
  * Chen et al. [2020] Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton.  A simple framework for contrastive learning of visual representations.  In _ICML_ , 2020.
  * Cheng et al. [2021] Bowen Cheng, Anwesa Choudhuri, Ishan Misra, Alexander Kirillov, Rohit Girdhar, and Alexander G Schwing.  Mask2former for video instance segmentation.  _arXiv preprint arXiv:2112.10764_ , 2021.
  * Cheng et al. [2022] Bowen Cheng, Ishan Misra, Alexander G Schwing, Alexander Kirillov, and Rohit Girdhar.  Masked-attention mask transformer for universal image segmentation.  In _CVPR_ , 2022.
  * Choudhuri et al. [2023] Anwesa Choudhuri, Girish Chowdhary, and Alexander G. Schwing.  Context-aware relative object queries to unify video instance and panoptic segmentation.  In _CVPR_ , 2023.
  * Dosovitskiy et al. [2021] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby.  An image is worth 16x16 words: Transformers for image recognition at scale.  In _ICLR_ , 2021.
  * Fischer et al. [2023] Tobias Fischer, Thomas Huang, Jiangmiao Pang, Linlu Qiu, Haofeng Chen, Trevor Darrell, and Fisher Yu.  Qdtrack: Quasi-dense similarity learning for appearance-only multiple object tracking.  _TPAMI_ , 2023.
  * Ghiasi et al. [2021] Golnaz Ghiasi, Yin Cui, Aravind Srinivas, Rui Qian, Tsung-Yi Lin, Ekin D. Cubuk, Quoc V. Le, and Barret Zoph.  Simple copy-paste is a strong data augmentation method for instance segmentation.  In _CVPR_ , 2021.
  * Girshick [2015] Ross Girshick.  Fast r-cnn.  In _ICCV_ , 2015.
  * Han et al. [2022] Su Ho Han, Sukjun Hwang, Seoung Wug Oh, Yeonchool Park, Hyunwoo Kim, Min-Jung Kim, and Seon Joo Kim.  Visolo: Grid-based space-time aggregation for efficient online video instance segmentation.  In _CVPR_ , 2022.
  * He et al. [2022] Fei He, Haoyang Zhang, Naiyu Gao, Jian Jia, Yanhu Shan, Xin Zhao, and Kaiqi Huang.  Inspro: Propagating instance query and proposal for online video instance segmentation.  In _NeurIPS_ , 2022.
  * He et al. [2016] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _CVPR_ , 2016.
  * He et al. [2017] Kaiming He, Georgia Gkioxari, Piotr Dollar, and Ross Girshick.  Mask r-cnn.  In _ICCV_ , 2017.
  * Heo et al. [2022] Miran Heo, Sukjun Hwang, Seoung Wug Oh, Joon-Young Lee, and Seon Joo Kim.  Vita: Video instance segmentation via object token association.  In _NeurIPS_ , 2022.
  * Heo et al. [2023] Miran Heo, Sukjun Hwang, Jeongseok Hyun, Hanjung Kim, Seoung Wug Oh, Joon-Young Lee, and Seon Joo Kim.  A generalized framework for video instance segmentation.  In _CVPR_ , 2023.
  * Huang et al. [2022] De-An Huang, Zhiding Yu, and Anima Anandkumar.  Minvis: A minimal video instance segmentation framework without video-based training.  In _NeurIPS_ , 2022.
  * Hwang et al. [2021] Sukjun Hwang, Miran Heo, Seoung Wug Oh, and Seon Joo Kim.  Video instance segmentation using inter-frame communication transformers.  In _NeurIPS_ , 2021.
  * Ke et al. [2021] Lei Ke, Xia Li, Martin Danelljan, Yu-Wing Tai, Chi-Keung Tang, and Fisher Yu.  Prototypical cross-attention networks for multiple object tracking and segmentation.  In _NeurIPS_ , 2021.
  * Kim et al. [2020] Dahun Kim, Sanghyun Woo, Joon-Young Lee, and In So Kweon.  Video panoptic segmentation.  In _CVPR_ , 2020.
  * Kuhn [1955] Harold W Kuhn.  The hungarian method for the assignment problem.  _NRL_ , 1955.
  * Li et al. [2022] Jizhizi Li, Jing Zhang, Stephen J. Maybank, and Dacheng Tao.  Bridging composite and real: Towards end-to-end deep image matting.  _IJCV_ , 2022.
  * Li et al. [2023] Junlong Li, Bingyao Yu, Yongming Rao, Jie Zhou, and Jiwen Lu.  Tcovis: Temporally consistent online video instance segmentation.  In _ICCV_ , 2023.
  * Lin et al. [2014] Tsung-Yi Lin, Michael Maire, Serge Belongie, James Hays, Pietro Perona, Deva Ramanan, Piotr Dollár, and C. Lawrence Zitnick.  Microsoft coco: Common objects in context.  In _ECCV_ , 2014.
  * Qi et al. [2022] Jiyang Qi, Yan Gao, Yao Hu, Xinggang Wang, Xiaoyu Liu, Xiang Bai, Serge Belongie, Alan Yuille, Philip HS Torr, and Song Bai.  Occluded video instance segmentation: A benchmark.  _IJCV_ , 2022.
  * Sohn [2016] Kihyuk Sohn.  Improved deep metric learning with multi-class n-pair loss objective.  In _NeurIPS_ , 2016.
  * Wang et al. [2020] Yuqing Wang, Zhaoliang Xu, Xinlong Wang, Chunhua Shen, Baoshan Cheng, Hao Shen, and Huaxia Xia.  End-to-end video instance segmentation with transformers.  In _CVPR_ , 2020.
  * Wu et al. [2022a] Junfeng Wu, Yi Jiang, Wenqing Zhang, Xiang Bai, and Song Bai.  Seqformer: a frustratingly simple model for video instance segmentation.  In _ECCV_ , 2022a.
  * Wu et al. [2022b] Junfeng Wu, Qihao Liu, Yi Jiang, Song Bai, Alan Yuille, and Xiang Bai.  In defense of online models for video instance segmentation.  In _ECCV_ , 2022b.
  * Wu et al. [2022c] Jialian Wu, Sudhir Yarram, Hui Liang, Tian Lan, Junsong Yuan, Jayan Eledath, and Gerard Medioni.  Efficient video instance segmentation via tracklet query and proposal.  In _CVPR_ , 2022c.
  * Yang et al. [2019] Linjie Yang, Yuchen Fan, and Ning Xu.  Video instance segmentation.  In _ICCV_ , 2019.
  * Yang et al. [2021a] Linjie Yang, Yuchen Fan, and Ning Xu.  The 3rd large-scale video object segmentation challenge - video instance segmentation track, 2021a.
  * Yang et al. [2022a] Linjie Yang, Yuchen Fan, and Ning Xu.  The 4th large-scale video object segmentation challenge - video instance segmentation track, 2022a.
  * Yang et al. [2021b] Shusheng Yang, Yuxin Fang, Xinggang Wang, Yu Li, Chen Fang, Ying Shan, Bin Feng, and Wenyu Liu.  Crossover learning for fast online video instance segmentation.  In _ICCV_ , 2021b.
  * Yang et al. [2022b] Shusheng Yang, Xinggang Wang, Yu Li, Yuxin Fang, Jiemin Fang, Wenyu Liu, Xun Zhao, and Ying Shan.  Temporally efficient vision transformer for video instance segmentation.  In _CVPR_ , 2022b.
  * Ying et al. [2023] Kaining Ying, Qing Zhong, Weian Mao, Zhenhua Wang, Hao Chen, Lin Yuanbo Wu, Yifan Liu, Chengxiang Fan, Yunzhi Zhuge, and Chunhua Shen.  Ctvis: Consistent training for online video instance segmentation.  In _ICCV_ , 2023.
  * Zhang et al. [2023] Tao Zhang, Xingye Tian, Yu Wu, Shunping Ji, Xuebo Wang, Yuan Zhang, and Pengfei Wan.  Dvis: Decoupled video instance segmentation framework.  In _ICCV_ , 2023.
  * Zhu et al. [2021] Xizhou Zhu, Weijie Su, Lewei Lu, Bin Li, Xiaogang Wang, and Jifeng Dai.  Deformable detr: Deformable transformers for end-to-end object detection.  In _ICLR_ , 2021.


# Detecting entanglement of unknown states by violating the Clauser-Horne-
Shimony-Holt inequality

J. Cortés-Vega<EMAIL_ADDRESS>Instituto Milenio de Investigación en
Óptica, Universidad de Concepción, Concepción, Chile  Facultad de Ciencias
Físicas y Matemáticas, Departamento de Física, Universidad de Concepción,
Concepción, Chile  J. F. Barra  Instituto Milenio de Investigación en Óptica,
Universidad de Concepción, Concepción, Chile  Facultad de Ciencias Físicas y
Matemáticas, Departamento de Física, Universidad de Concepción, Concepción,
Chile  L. Pereira  Instituto de Física Fundamental IFF-CSIC, Calle Serrano
113b, Madrid 28006, Spain  A. Delgado  Instituto Milenio de Investigación en
Óptica, Universidad de Concepción, Concepción, Chile  Facultad de Ciencias
Físicas y Matemáticas, Departamento de Física, Universidad de Concepción,
Concepción, Chile

###### Abstract

Entangled states play a fundamental role in Quantum Mechanics and are at the
core of many contemporary applications, such as quantum communication and
quantum computing. Therefore, determining whether a state is entangled or not
is an important task. Here, we propose a method to detect the entanglement of
unknown two-qubit quantum states. Our method is based on the violation of the
Clauser-Horne-Shimony-Holt inequality. This maximizes the value of the
inequality even when it contains an unknown quantum state. The method
iteratively generates local measurement settings that lead to increasing
values of the inequality. We show by numerical simulations for pure and mixed
states that our algorithm exceeds the classical limit of 2 after a few
iterations.

###### pacs:

03.67.-a, 0365.-w, 02.60.Pn

##  I Introduction

Quantum mechanics predicts the existence of quantum states of composite
systems that cannot be written as products of states of their individual
components [1]. These are the so called entangled states. Today, these states
play a central role in quantum information theory [3, 2] and in many
applications, such as, for instance, quantum cryptography [4], quantum
teleportation [5, 6], frequency standards improvement [7, 8, 9], one-way
quantum computing [10], clock synchronization [11], and entanglement assisted
orientation in space [12], among many others. Interestingly, entangled states
play a key role in the argument put forward by Einstein, Podolsky, and Rosen
[13]. This was aimed at ascribing objective values to measurable quantities,
that is, values that exist prior to and independently of measurements. Bell’s
inequality [14] shows that precisely the existence of entangled states
precludes such conception of reality.

In view of the foundational significance of entangled states and their many
applications, theoretical and experimental characterization and detection of
entangled states are important research subjects. One of the first criteria
employed to study the entanglement of quantum states is the violation of the
Clauser-Horne-Shimony-Holt inequality [15, 16] (CHSH), which is the
generalization of Bell’s inequality to two observers each having the choice of
two measurement settings with two outcomes. In this scenario, the violation of
the CHSH inequality indicates the presence of entanglement. This approach has
also been studied in the context of the theory of entanglement witnesses [17,
18]. These are observables with positive expectation values with respect to
the complete set of separable states that for at least one entangled state
provide a negative expectation value. Thus, a negative expectation value
signals the presence of entanglement. It has been shown that the CHSH
inequality can be related to an entanglement witness [18, 19].

Here, we study the detection of entanglement of unknown states via the
violation of the CHSH inequality. Since the majority of the entanglement
measures and entanglement detectors are based on the knowledge of the quantum
state, the unknown character of the state increases the difficulty of the
problem. The presence of unknown quantum states is common in quantum
communication [20, 21, 22] and quantum computing [23, 24, 25], where an
objective entangled state is prepared, but it is modified by the action of the
environment. Entanglement detection of unknown quantum states has been
previously studied from the point of view of quantum tomography [26, 27] by
means of an adaptive scheme [28, 29], employing a succession of measurements
of witness operators [30, 31], via the measurement of the energy observable
[32], via local parity measurements on two-fold copies of the unknown state
[33], series of local random measurements from which entanglement witnesses
are constructed [34], and variational determination of geometrical
entanglement [35], among many others. We follow a different approach. For a
given known state, the maximal violation of the CHSH inequality is obtained by
maximizing the inequality onto the set of 4-tuples of dichotomic observables.
This procedure is typically carried out by means of semidefinite programming
(SDP) techniques. If the state is unknown, then the function to be optimized,
that is, the target function, contains unknown fix parameters and SDP cannot
be employed to find the measurements leading to the maximal violation.
Analogously, the use of an entanglement witness also requires the knowledge
about the state. To overcome this problem we employ a recently developed
optimization algorithm [36], the Complex simultaneous perturbation stochastic
approximation (CSPSA), which can handle functions with unknown parameters.
CSPSA works natively within the field of the complex numbers. Thereby, no
parameterization of the complex arguments onto the real numbers is necessary.
Also, this algorithm has exhibited an improved convergence rate in certain
applications such as, for instance, the estimation of unknown quantum pure
states [37]. CSPSA uses a stochastic approximation of the complex Wirtinger
gradient of the target function, that is, the function to be optimized, which
requires the value of the target function at two different points in the
optimization space. In the case at hand, these two values can be obtained
experimentally, regardless of whether the state remains unknown. CSPSA
iteratively generates a sequence of sets with four local measurement settings
with increasing values of the CHSH function until reaching the highest
possible violation of the inequality.

We first study via numerical simulations the performance of the method here
proposed when applied to unknown pure 2-qubit states. In this case, the
maximal value achieved by the CHSH function depends on the Schmidt coefficient
of the state. Thereby, the performance of the method can be compared with an
analytical bound. We show that for the set formed by states that have the same
set of local Schmidt bases, the method leads in tens of iterations to a value
close to the maximum of the CHSH function for each value of the Schmidt
coefficient. We also consider sets of states that have the same concurrence
value but different local Schmidt bases. In this case, the method also
approaches the corresponding maximum value of the CHSH inequality in tens of
iterations. However, the higher the concurrence value, the fewer iterations
are required for a violation of the CHSH inequality. Also, all states with the
same concurrence value exhibit a very similar behavior of the CHSH function as
a function of the number of iterations, that is, CSPSA produces results that
are nearly independent of the particular set of local Schmidt bases. We also
consider the average behavior of the method on the Hilbert space of two
qubits. In this case, the method reaches a CHSH function value greater than 2
after 17 iterations for an ensemble size of $10^{2}$. After 25 iterations the
interquartile range is also above 2, which indicates that for 75% of the
simulated states the method reached a violation of the CHSH inequality. A
further increase of the ensemble size leads to a reduction in the number of
iterations required to achieve a violation of the CHSH inequality. In order to
study the accuracy achieved by our method we employ the squared error. We show
that the mean and median squared error on the 2-qubit Hilbert space are nearly
indistinguishable. After 25 iterations the mean square error achieves a value
in the order of $10^{-1}$ for an ensemble size of $10^{2}$. A further increase
of the ensemble size to $10^{3}$ leads to a decrease in mean square error in
the order of half order of magnitude. Thereafter, we study the case of two-
qubit mixed states. Unlike the case of pure states, there is no known
analytical formula for the maximum value of the CHSH function for an arbitrary
mixed state. However, in the particular case of Werner states, that is, a
maximally entangled state affected by white noise, it is possible to obtain
the maximum value of the CHSH function in terms of the mixing parameter. We
show that CSPSA is capable of achieving a value close to the maximum violation
of the CHSH inequality for all Werner states. As the ensemble size increases
the value of the function provided by CSPSA becomes closer to the maximal
violation. Finally, we analyze the results achieved by CSPSA for unknown mixed
states. For these states there is no analytical expression for maximal
violation, so we calculate this value via SDP. After generating $10^{6}$
density matrices, a subset of $8\times 10^{3}$ density matrices that violate
the CHSH inequality is identified. These states have a small value of the
negativity, a well-known entanglement measure. Within this subset, the mean
and median values of the CHSH function provided by CSPSA achieve a value close
to the theoretical maximal violation after approximately 75 iterations.

Our results show that the maximization of the CHSH function via the CSPSA
method allows detecting the entanglement of unknown states, pure or mixed,
with a high degree of accuracy. Furthermore, the highest value of the CHSH
function can also be achieved. Our approach requires the ability to adapt
local measurements, which are carried out on single copies of the unknown
state. This can be implemented in various experimental platforms [16, 38, 39,
40, 41, 42]. We stress the fact that no a priori information about the unknown
state, such as purity, Schmidt coefficient, or Schmidt bases, has been
employed to optimize the performance of CSPSA.

##  II CHSH inequality and CSPSA optimization algorithm

The target function to be optimized is the Clauser-Horne-Shimony-Holt function
$S$ defined by the expression [15]

$\displaystyle S({\bm{z}},{\bm{z}}^{*})$ $\displaystyle=$ $\displaystyle
E({\bm{z}}_{a},{\bm{z}}_{b})+E({\bm{z}}_{a},{\bm{z}}^{\prime}_{b})+E({\bm{z}}^{\prime}_{a},{\bm{z}}_{b})$
(1) $\displaystyle-$ $\displaystyle
E({\bm{z}}^{\prime}_{a},{\bm{z}}^{\prime}_{b}),$

where the expectation value $E({\bm{z}}_{a},{\bm{z}}_{b})$ is given by the
average of the products of the outcomes of two locally performed dichotomic
measurements $A({\bm{z}}_{a})$ and $B({\bm{z}}_{b})$ defined by the settings
${\bm{z}}_{a}$ and ${\bm{z}}_{b}$, respectively. The vector ${\bm{z}}$
contains the settings of the four local measurements, that is,
$\bm{z}=({\bm{z}}_{a},{\bm{z}}^{\prime}_{a},{\bm{z}}_{b},{\bm{z}}^{\prime}_{b})$.
The CHSH inequality adopts the form $|S|\leq 2$.

A quantum mechanical dichotomic observable $A({\bm{z}}_{a})$ is defined as the
one having $\pm 1$ eigenvalues, that is, an observable with the spectral
decomposition

$A({\bm{z}}_{a})=|\psi({\bm{z}}_{a})\rangle\langle\psi({\bm{z}}_{a})|-|\psi^{\perp}({\bm{z}}_{a})\rangle\langle\psi^{\perp}({\bm{z}}_{a})|,$
(2)

where $|\psi({\bm{z}}_{a})\rangle$ is an arbitrary two-dimensional quantum
state

$|\psi({\bm{z}}_{a})\rangle=\frac{z_{a,1}|0\rangle+z_{a,2}|1\rangle}{\sqrt{|z_{a,1}|^{2}+|z_{a,2}|^{2}}}.$
(3)

The state $|\psi^{\perp}({\bm{z}}_{a})\rangle$ is orthogonal to
$|\psi({\bm{z}}_{a})\rangle$ and the components $z_{a,1}$ and $z_{a,2}$ of the
vector ${\bm{z}}_{a}$ are complex numbers. Thereby, the expectation value
$E({\bm{z}}_{a},{\bm{z}}_{b})$ is given by the expression

$\displaystyle E({\bm{z}}_{a},{\bm{z}}_{b})$ $\displaystyle=$ $\displaystyle
Tr(\rho|\psi({\bm{z}}_{a})\rangle\langle\psi({\bm{z}}_{a})|\otimes|\psi({\bm{z}}_{b})\rangle\langle\psi({\bm{z}}_{b})|)$
(4) $\displaystyle+$ $\displaystyle
Tr(\rho|\psi^{\perp}({\bm{z}}_{a})\rangle\langle\psi^{\perp}({\bm{z}}_{a})|\otimes|\psi^{\perp}({\bm{z}}_{b})\rangle\langle\psi^{\perp}({\bm{z}}_{b})|)$
$\displaystyle-$ $\displaystyle
Tr(\rho|\psi({\bm{z}}_{a})\rangle\langle\psi({\bm{z}}_{a})|\otimes|\psi^{\perp}({\bm{z}}_{b})\rangle\langle\psi^{\perp}({\bm{z}}_{b})|)$
$\displaystyle-$ $\displaystyle
Tr(\rho|\psi^{\perp}({\bm{z}}_{a})\rangle\langle\psi^{\perp}({\bm{z}}_{a})|\otimes|\psi({\bm{z}}_{b})\rangle\langle\psi({\bm{z}}_{b})|),$

where $\rho$ is a fixed known two-qubit state.

The problem of violating the CHSH inequality consists in finding a complex
vector ${\bm{z}}$ such that for a given known state $\rho$ leads to a maximal
value of $|S({\bm{z}},{\bm{z}}^{*})|$ larger than the classical bound of 2.
This optimization problem can be solved by means of semidefinite programing or
other numerical optimization techniques. However, when the state $\rho$
entering in the function $S$ is unknown, the standard approaches to the
problem cannot be employed. The reason for this is that the function $S$ and
its derivatives cannot be evaluated.

Algorithm 1 CSPSA optimization of S($\rho,{\bm{z}}$)

Consider a state $\rho$. This plays the role of the unknown state.

Set initial guess $\hat{{\bm{z}}}_{0}$, and gain coefficients $a$, $A$, $s$,
$b$, and $r$.

for $k=1,\dots,k_{max}$ do

Set  $a_{k}=\frac{a}{(k+1+A)^{s}},\quad c_{k}=\frac{b}{(k+1)^{r}}.$

Choose $\Delta_{k,i}$ randomly in the set $\\{\pm 1,\pm i\\}$.

Calculate $\hat{{\bm{z}}}_{k\pm}=\hat{{\bm{z}}}_{k}\pm
c_{k}{\bm{\Delta}_{k}}$.

Estimate from experimentally acquired data or numerically simulate the values
$S(\rho,\hat{{\bm{z}}}_{k\pm})$ considering an ensemble of $N$ equally
prepared pairs of qubits in the state $\rho$.

Estimate the gradient as
$\hat{g}_{k,i}=\frac{S(\rho,\hat{{\bm{z}}}_{k+})-S(\rho,\hat{{\bm{z}}}_{k-})}{2c_{k}\Delta_{k,i}^{*}}.$

Actualize the guess
$\hat{{\bm{z}}}_{k+1}=\hat{{\bm{z}}}_{k}+a_{k}\hat{{\bm{g}}}_{k}$.

Normalize coefficients $\hat{{\bm{z}}}_{k+1}$

In order to overcome this problem, we resort to the recently introduced CSPSA
[36] optimization algorithm for real-valued functions of complex arguments.
This algorithm works natively on the field of the complex numbers, which make
unnecessary the use of real parameterizations of the complex arguments. For a
target function
$f({\bm{z}},{\bm{z}}^{*}):\mathbb{C}^{n}\times\mathbb{C}^{n}\rightarrow\mathbb{R}$,
CSPSA is defined by the iterative rule

$\hat{\bm{z}}_{k+1}=\hat{\bm{z}}_{k}+a_{k}\hat{\bm{g}}_{k}(\hat{\bm{z}}_{k},\hat{\bm{z}}_{k}^{*}),$
(5)

where $a_{k}$ is a positive gain coefficient and $\hat{\bm{z}}_{k}$ is the
estimate of the maximizer $\tilde{\bm{z}}$ of $f({\bm{z}},{\bm{z}}^{*})$ at
the k-th iteration. The iteration starts from an initial guess
$\hat{\bm{z}}_{0}$, which is randomly chosen. The function
$\hat{\bm{g}}_{k}(\hat{\bm{z}}_{k},\hat{\bm{z}}_{k}^{*})$ is an estimator for
the Wirtinger gradient [43] of $f({\bm{z}},{\bm{z}}^{*})$ whose components are
defined by

$\hat{g}_{k,i}=\frac{f(\hat{\bm{z}}_{k+},\hat{\bm{z}}_{k+}^{*})+\epsilon_{k,+}-(f(\hat{\bm{z}}_{k-},\hat{\bm{z}}_{k-}^{*})+\epsilon_{k,-})}{2c_{k}{\Delta}_{k,i}^{*}},$
(6)

with

$\hat{\bm{z}}_{k\pm}=\hat{\bm{z}}_{k}\pm c_{k}{\bm{\Delta}}_{k},$ (7)

where $c_{k}$ is a positive gain coefficient and $\epsilon_{k,\pm}$ describes
the presence of noise in the values of
$f(\hat{\bm{z}}_{k\pm},\hat{\bm{z}}_{k\pm}^{*})$. The components of the vector
${\bm{\Delta}}_{k}\in\mathbb{C}^{n}$ are identically and independently
distributed random variables in the set $\\{\pm 1,\pm i\\}$. The gain
coefficients $a_{k}$ and $c_{k}$ control the convergence of CSPSA and are
chosen as

$a_{k}=\frac{a}{(k+1+A)^{s}},~{}~{}c_{k}=\frac{b}{(k+1)^{r}}.$ (8)

The values of $a,A,s,b$ and $r$ are adjusted to optimize the rate of
convergence depending on the target function. We use the values: $a=1.0$,
$b=0.25$, $s=1.0$, $r=1/6$, and $A=0$.

Two main properties of CSPSA are: (i) it converges asymptotically in mean to
the maximizer $\tilde{\bm{z}}$ of $f({\bm{z}},{\bm{z}}^{*})$ and (ii)
$\hat{\bm{g}}_{k}$ is an asymptotically unbiased estimator of the Wirtinger
gradient. With proper conditions, these properties are maintained even in the
presence of the noise terms $\epsilon_{k,\pm}$ entering in Eq. (6). CSPSA is
the generalization of the Simultaneous perturbation stochastic approach (SPSA)
[44, 45] from the field of real numbers to the field of complex numbers. SPSA
has been applied to the problem of estimating pure states [46, 47] and
experimentally realized [48].

Thus, the application of CSPSA to the maximization of the CHSH function
proceeds as follows: an initial guess $\hat{\bm{z}}_{0}$ for the vector
containing the measurement settings and a vector ${\bm{\Delta}}_{0}$ are
randomly generated. These two vectors are employed to calculate the vectors
$\hat{\bm{z}}_{0\pm}$ according to Eq. (7). Thereafter, the values
$S(\hat{\bm{z}}_{0\pm},\hat{\bm{z}}^{*}_{0\pm})$ of the CHSH function are
obtained, which involves the realization of measurements on a finite ensemble
of $N$ copies of the unknown state $\rho$. The values
$S(\hat{\bm{z}}_{0\pm},\hat{\bm{z}}^{*}_{0\pm})$ are then employed to
calculate the estimator for the Wirtinger gradient
$\hat{\bm{g}}_{0}(\hat{\bm{z}}_{0},\hat{\bm{z}}_{0}^{*})$ using Eq. (6).
Finally, a new estimate $\hat{\bm{z}}_{1}$ for the vector of settings is
obtained by means of Eq. (5). This process is iterated until achieving a
violation of the CHSH inequality or until reaching a predefined number of
iterations. Algorithm 1 shows a pseudocode for the optimization of the CHSH
function via CSPSA.

According to Eq. (6), the use of CSPSA to maximize the CHSH function
$S({\bm{z}},{\bm{z}}^{*})$ requires the capability of obtaining the values
$S(\hat{\bm{z}}_{k\pm},\hat{\bm{z}}_{k\pm}^{*})$ at each iteration, which in
turn requires an experimental platform capable of measuring the CHSH function
at any value of the setting vector ${\bm{z}}$. In photonic platforms where a
qubit is encoded in the polarization degree of freedom of a single photon, the
local measurements on a qubit are carried out by the interaction of the photon
with a sequence of half- and quarter-wave plates followed by a polarizing beam
splitter and single-photon detectors. In this case a setting vector is given
by the rotation angles of the wave plates. Thereby, it is possible to
implement any local measurement up to the angular resolution of the wave
plates. It is possible to achieve a high degree of control in other
experimental platforms, for instance in time-bin or energy-time encoded
qubits, where local measurements can be implemented introducing electronically
controlled phase shifts. Thus, we will assume that the CHSH function can be
measured for any value of the setting vector ${\bm{z}}$.

##  III Results

A single run of CSPSA starts with the choice of an initial guess
${\bm{z}}_{0}$ of the four local measurement bases and proceeds through the
choice of the vector ${\bm{\Delta}}_{k}$ at every iteration. Since there is no
a priori information about the initial state, the initial guess for each of
the local measurements, which are defined by Eqs. (2) and (3), is randomly
chosen according to a Haar-uniform distribution. The choice of
${\bm{\Delta}}_{k}$ is equally random. Thereby, CSPSA is an intrinsically
stochastic optimization algorithm. A third source of randomness is the value
of the CHSH function. This is obtained by means of probabilities that are
inferred from local measurements made on a set of equally prepared copies of
the unknown state. Since the size $N$ of the ensemble is finite, the inferred
probabilities are affected by finite statistic noise. Thereby, CSPSA exhibits
three different sources of randomness and, consequently, each run of CSPSA
will follow a different trajectory in the optimization space, that is, the
space of all four setting vectors. Here, we report the results of numerical
experiments for the cases of pure and mixed states considering the sources of
randomness affecting the performance of the proposed method.

To study the violation of the CHSH inequality with an unknown state $\rho$,
pure or mixed, we compute the expected value $\bar{S}(\rho)$ by sampling a
sufficiently large number of independent trajectories, each obtained through
the optimization of $S$ by CSPSA for $\rho$, as

$\bar{S}(\rho)=\frac{1}{K}\sum_{{\bm{z}}_{0},\\{{\bm{\Delta}}_{1},\dots,{\bm{\Delta}}_{k}\\}}S(\rho,{\bm{z}}_{0},\\{{\bm{\Delta}}_{1},\dots,{\bm{\Delta}}_{k}\\}),$
(9)

where $S(\rho,{\bm{z}}_{0},\\{{\bm{\Delta}}_{1},\dots,{\bm{\Delta}}_{k}\\})$
is the value of the CHSH function evaluated on a particular trajectory
generated by a single run of CSPSA and $K$ is the total number of simulated
trajectories.
$S(\rho,{\bm{z}}_{0},\\{{\bm{\Delta}}_{1},\dots,{\bm{\Delta}}_{k}\\})$ depends
on the unknown state $\rho$, the set ${\bm{z}}_{0}$ of complex numbers that
defines the initial guess for the four local measurements, and the particular
sequence of choices $\\{{\bm{\Delta}}_{1},\dots,{\bm{\Delta}}_{k}\\}$. The
mean $\bar{S}(\rho)$ will be studied as a function of the number $k$ of
iterations for a fixed ensemble size $N$.

Since we are interested in the overall behavior of the algorithm for unknown
states, we calculate the mean $\bar{S}_{C}$ of $\bar{S}(\rho)$ in a set
$\Omega_{C}$, that is,

$\bar{S}_{C}=\frac{1}{M}\sum_{\rho\in\Omega_{C}}\bar{S}(\rho),$ (10)

where $M$ is the number of states in $\Omega_{C}$ and $C$ is a parameter that
characterizes the states in the set. Alternatively, we calculate the median
$\tilde{S}_{C}$ of $\bar{S}(\rho)$ in the set $\Omega_{C}$ and the
interquartile range. This is done to determine whether the distribution of
$\bar{S}(\rho)$ in $\Omega_{C}$ exhibits a symmetric distribution or not and
the possible existence of outliers.

Figure 1: CHSH function $S(|\psi_{\lambda}\rangle)$ as a function of the
Schmidt coefficient $\lambda$ for two-qubit states with fixed local Schmidt
bases. Continuos green line represents the theoretical prediction given by Eq.
(12). Solid red circles (blue x’s) represent the mean
$\bar{S}(|\psi_{\lambda}\rangle)$ (median $\tilde{S}(|\psi_{\lambda}\rangle)$)
of $S(|\psi_{\lambda}\rangle)$ obtained via CSPSA considering $10^{4}$ initial
guesses for each state $|\psi_{\lambda}\rangle$, 200 iterations, and an
ensemble size $N=10^{2}$. Figure 2: CHSH function $S(|\psi_{\lambda}\rangle)$
as a function of the Schmidt coefficient $\lambda$ for two-qubit states with
fixed local Schmidt bases. Continuos green line represents the theoretical
prediction given by Eq. (12). Solid red circles (blue x’s) represent the mean
$\bar{S}(|\psi_{\lambda}\rangle)$ (median $\tilde{S}(|\psi_{\lambda}\rangle)$)
of $S(|\psi_{\lambda}\rangle)$ obtained via CSPSA considering $10^{4}$ initial
guesses for each state $|\psi_{\lambda}\rangle$, 200 iterations, and an
ensemble size $N=10^{4}$.

###  III.1 Unknown pure states

We start our analysis of the proposed algorithm by considering the violation
of the CHSH inequality for the set $\Omega_{\lambda}$ of two-qubit pure states
defined by the Schmidt decomposition

$|\psi(\lambda)\rangle=\sqrt{\lambda}|0\rangle_{1}|0\rangle_{2}+\sqrt{1-\lambda}|1\rangle_{1}|1\rangle_{2},$
(11)

where $\lambda\in[0,1/2]$ is the Schmidt coefficient and
$\\{|0\rangle_{1},|1\rangle_{1}\\}$ and $\\{|0\rangle_{2},|1\rangle_{2}\\}$
are fixed local Schmidt bases of each qubit. States in $\Omega_{\lambda}$ lead
to a value of the function $S$ given by

$S(\lambda)=2\sqrt{1+4\lambda(1-\lambda)}.$ (12)

In Fig. 1 we show $\bar{S}(\rho_{\lambda})$ for
$\rho_{\lambda}=|\psi_{\lambda}\rangle\langle\psi_{\lambda}|$ as a function of
$\lambda$ for $N=10^{2}$ after 200 iterations and $K=10^{4}$. Initial guesses
for the set of four local observables are randomly chosen. In particular,
information about the fixed bases in $|\psi_{\lambda}\rangle$ has not been
used to improve the performance of CSPSA. As is apparent from Fig. 1, CSPSA
provides mean and median of $S(|\psi_{\lambda}\rangle)$ that closely resemble
the theoretical prediction of Eq. (12) for any value of $\lambda$. A much
better agreement can be obtained by increasing the ensemble from $N=10^{2}$ to
$N=10^{4}$, which is illustrated in Fig. 2.

Figure 3: Mean $\bar{S}_{C}$ of $\bar{S}(\rho)$ in $\Omega_{C}$ as a function
of the number of iterations for several values of the concurrence $C$ in the
interval $[0.1,1.0]$, from bottom to top. The mean $\bar{S}(\rho)$ is
calculated with $10^{4}$ independent trajectories and each local measurement
is simulated with an ensemble size $N=10^{2}$. Upper and lower straight lines
represent the values 2$\sqrt{2}$ and 2, correspondingly. Figure 4: Median of
$\bar{S}(\rho)$ in $\Omega_{C}$ as a function of the number of iterations for
several values of the concurrence $C$ in the interval $[0.1,1.0]$, from bottom
to top. The mean $\bar{S}(\rho)$ is calculated with $10^{4}$ independent
trajectories and each local measurement is simulated with an ensemble size
$N=10^{2}$. Upper and lower straight lines represent the values $2\sqrt{2}$
and 2, correspondingly.

Next we analyze the case of pure states with a known value of the concurrence
$C$, which is given by the expression

$C(\lambda)=2\sqrt{\lambda}\sqrt{1-\lambda}.$ (13)

The local Schmidt bases of the state are unknown. In the simulations we choose
a fixed value $C$ of the concurrence, which in turn fixes the value of the
Schmidt coefficient. The local Schmidt bases are randomly chosen. As in the
previous simulations, the knowledge about the value of the concurrence is not
employed to improve the performance of CSPSA. Figure 3 shows the behavior of
$\bar{S}_{C}$, which is the mean of $\bar{S}(\rho)$ calculated on a set
$\Omega_{C}$ of pure states with a fixed value $C$ of the concurrence, as a
function of the number $k$ of iterations for several values of $C$. Each set
$\Omega_{C}$ contains 100 states chosen according to a Haar-uniform
distribution and $\bar{S}(\rho)$ is calculated with $10^{4}$ trajectories.
Each one of the four local measurements is simulated considering an ensemble
size of $N=10^{2}$. According to Fig. 3, the quantity $\bar{S}_{C}$ exhibits a
fast increase of the value of the CHSH function within the first tens of
iterations followed by a linear behavior, which asymptotically approaches the
maximal value of the function $S$ for the value $C$ of the concurrence. The
overall behavior of $\bar{S}_{C}$ does not depend on the value of $C$.

Figure 4 displays the median $\tilde{S}_{C}$ of $\bar{S}(\rho)$ in
$\Omega_{C}$ as a function of the number of iterations for several values of
the concurrence $C$. Shaded areas represent the interquartile range. Monte
Carlo experiments are carried out as in Fig. 3. As is apparent from this
figure, the median exhibits the same overall behavior as the mean
$\bar{S}_{C}$. Mean and median reach after a few tens interations values that
are nearly indistinguishable and contained within the interquartile range.
This indicates that the stochasticity of CSPSA does not lead to outliers in
the histogram of $\bar{S}(\rho)$ for all simulated sets $\Omega_{C}$. The
interquartile range, which is a quartile-based measure of variability,
decreases rapidly with the number of iterations and becomes a very narrow
fringe. This is an indication that the histogram of $\bar{S}(\rho)$ for a
particular $\Omega_{C}$ after a few tens iterations is highly concentrated
around the mean.

Thus, Figs. 3 and 4 clearly indicate that CSPSA can be employed to iteratively
increase the value of the CHSH function for unknown pure states and detect
entanglement. The greater the entanglement of the unknown state, the fewer
iterations will be required to obtain a violation of the CHSH inequality.
Furthermore, approximately 70 iterations are necessary to reach a value of the
CHSH function close to the maximal violation allowed by quantum mechanics.

Figure 5: Mean $\bar{S}_{C}$ of $\bar{S}(\rho)$ in $\Omega_{C}$ as a function
of the number of iterations for several values of the concurrence $C$ in the
interval $[0.1,1.0]$, from bottom to top. The mean $\bar{S}(\rho)$ is
calculated with $10^{4}$ independent trajectories and each local measurement
is simulated with an ensemble size $N=10^{4}$. Upper and lower straight lines
represent the values $2\sqrt{2}$ and 2, correspondingly. Figure 6: Median of
$\bar{S}(\rho)$ in $\Omega_{C}$ as a function of the number of iterations for
several values of the concurrence $C$ in the interval $[0.1,1.0]$, from bottom
to top. The mean $\bar{S}(\rho)$ is calculated with $10^{4}$ independent
trajectories and each local measurement is simulated with an ensemble size
$N=10^{4}$. Upper and lower straight lines represent the values $2\sqrt{2}$
and 2, correspondingly.

Figs. 5 and 6 depicts the mean $\bar{S}_{C}$ and the median $\tilde{S}_{C}$ of
$\bar{S}(\rho)$ in $\Omega_{C}$, correspondingly. In this case local
measurements are simulated with an ensemble size of $N=10^{4}$, that is, a
quadratic increase with respect to previous simulations. As is apparent from
Figs. 5 and 6, the overall behavior remains unchanged with respect to Figs. 3
and 4. In particular, both values of ensemble size, $N=10^{2}$ and $N=10^{4}$,
show small differences in the asymptotic linear regime. For instance, for
weakly entangled states, that is, $C=0.1$, after the total of iterations, in
Fig. 4 CSPSA is close to 2 but below. In Fig. 6, CSPSA is slightly above 2.
Similar differences can be observed for other values of $C$. Furthermore, a
small reduction in the number of iterations required to violated the CHSH
inequality can be observed. This reduction depends on the initial amount of
entanglement of the unknown state. Also, the increase in $N$ leads to narrower
interquartile ranges.

This is more clearly illustrated in Fig. 7, which shows the median
$\tilde{S}_{C}$ of $S(\rho)$ in $\Omega_{C}$ for $C=0.5$ and $C=0.9$ for three
values of ensemble size $N=10^{2},10^{3},10^{4}$. The interquartile range is
also depicted. As is apparent from Fig. 7, CSPSA provides very similar values
of $\tilde{S}_{C}$ almost independently of the size of the ensemble employed.
However, in the regime of a few tens of iterations, $N=10^{2}$ leads to lower
values of $\tilde{S}_{C}$, while $N=10^{3}$ and $10^{4}$ lead to very similar
values of $\tilde{S}_{C}$, which are higher than in the case $N=10^{2}$. This
has for consequence that higher values of $N$ lead to a decrease in the number
of iterations required to observe a violation of the CHSH inequality, but this
improvement is saturated for an enough large sample size.

Figure 7: Median $\tilde{S}_{C}$ of $\bar{S}(\rho)$ in $\Omega_{C}$ as a
function of the number of iterations for $C=0.5$ and $C=0.9$. Each local
measurement is simulated with an ensemble size $N=10^{4},10^{3},10^{2}$. The
median $\tilde{S}(\rho)$ for each value of $C$ is calculated with $10^{4}$
independent trajectories. Upper and lower straight lines represent the values
$2\sqrt{2}$ and 2, correspondingly.

This later effect is analyzed with the help of Fig. 8 that displays the number
of iterations $k_{S>2}$ required to obtain a violation of the inequality with
$75\%$ of the states generated for a given value of $C$ and with
$N=10^{2},10^{3},10^{4}$. Here, we observe that $N=10^{4}$ and $N=10^{3}$ lead
to a very similar behavior while $N=10^{2}$ requires the largest number of
iterations to reach a violation of the CHSH inequality. Also, the lower the
concurrency value, the greater the number of iterations required for the
violation. In fact, Fig. 8 suggests that $k_{S>2}$ decreases exponentially
with $C$. This figure also illustrates the interplay between $k_{S>2}$ and the
total ensemble size $N_{S>2}$ required for violating the CHSH inequality. For
example, in the case of $C=0.1$ and $N=10^{2}$, we have that approximately
$k_{S>2}=100$, which leads to $N_{S>2}=8\times 10^{4}$. For $N=10^{4}$ we have
that approximately $k_{S>2}=35$ and thus $N_{S>2}=280\times 10^{4}$. Clearly,
the reduction in the value of $k_{S>2}$ comes at the expense of using a much
larger total ensemble $N_{S>2}$. For states with a high value of concurrence
$C$, the reduction in the value of $k_{S>2}$ by increasing the value of $N$ is
marginal.

Figure 8: Number of iterations $k_{S>2}$ such that the interquartile range is
above $S=2$ as a function of the concurrence $C$ for $N=10^{2},10^{3}$, and
$10^{4}$, from top to bottom. Figure 9: Mean $\bar{S}_{\cal H}$ and median
$\tilde{S}_{\cal H}$ of $\bar{S}(|\psi\rangle\langle\psi|)$ with
$|\psi\rangle\in\cal H$ and interquartile range for $N=10^{2}$. Figure 10:
Mean $\bar{S}_{\cal H}$ and median $\tilde{S}_{\cal H}$ of
$\bar{S}(|\psi\rangle\langle\psi|)$ with $|\psi\rangle\in\cal H$ and
interquartile range for $N=10^{4}$.

So far, our study of the violation of CHSH inequality through CSPSA has been
done considering that the initial amount of entanglement is known. This was
done to show that CSPSA drives the value of the CHSH function $S$ close to the
maximum value regardless of the amount of entanglement. We now lift this
assumption and consider unknown pure states. In order to do this, we generate
a set $\Omega_{\cal H}$ with 100 pure states in the Hilbert space ${\cal
H}={\cal H}_{1}\otimes{\cal H}_{2}$ of two qubits according to a Haar-uniform
distribution and calculate the mean $\bar{S}_{\cal H}$ and the median
$\tilde{S}_{\cal H}$ of $\bar{S}(|\psi\rangle\langle\psi|)$ in $\Omega_{\cal
H}$, together with the corresponding interquartile range. These quantities are
depicted in Fig. 9 as a function of the number of iterations. The behavior
exhibited by the mean and media is very similar and characterized by a fast
increase within the first tens of iterations followed by an asymptotic linear
regime. Fig. 9 also shows the mean and media of the maximal theoretical values
of $S$ for each state in $\Omega_{\cal H}$, which are indicated as two
superposed straight lines. As can be seen from Fig. 9, CSPSA produces a mean
and a median that are very closely to the theoretical values. Also, the
expected number of iterations $k_{S>2}$ such that 75% of the simulated states
violates the CHSH inequality is about 25. Fig. 10 shows the same information
as Fig. 9 but with $N=10^{4}$. In this case, we see that the quadratic
increase in the ensemble size allows CSPSA to reach mean and media values that
are even closer to the theoretical values. Furthermore, there is a small
reduction in the number of iterations required to obtain a value of $S$
greater than two from 25 to 20.

Our previous simulations seem to indicate that the optimization of the CHSH
function for an unknown state through the CSPSA method provides maximum values
of the CHSH functional close to the theoretical maximum values. In order to
analyze this we employ the mean square error. For a given state
$\rho=|\psi\rangle\langle\psi|$ and a single realization of CSPSA we calculate
the square error $SE(\rho)$ as

$SE(\rho)=|S(\rho,{\bm{z}}_{0},\\{{\bm{\Delta}}_{1},\dots,{\bm{\Delta}}_{k}\\}-S_{max}(\rho)|^{2}.$
(14)

The mean square error $MSE(\rho)$ for a fixed unknown state $\rho$ with
respect to a large set of realizations is given by

$MSE(\rho)=\frac{1}{K}\sum_{{\bm{z}}_{0},\\{{\bm{\Delta}}_{1},\dots,{\bm{\Delta}}_{k}\\}}SE(\rho),$
(15)

which corresponds to an estimation accuracy metric. This is then used to
calculate the average of the mean square error $\overline{MSE}$ on the total
Hilbert space $\cal H$ as

$\overline{MSE}=\frac{1}{M}\sum_{\rho\in\Omega_{\cal H}}MSE(\rho).$ (16)

Figure 11: Mean square error $\overline{MSE}$ as a function of the number $k$
of iterations for $N=10^{2},10^{3}$, and $10^{4}$, from top to bottom. Shaded
areas represent interquartile range.

Figure 11 shows the mean $\overline{MSE}$ of the square error on the Hilbert
space as a function of the number of iterations for $N=10^{2},10^{3},10^{4}$.
For each value of ensemble size, $\overline{MSE}$ displays a fast decrease
followed by an approximately asymptotic lineal behavior. $N=10^{3}$ and
$N=10^{4}$ produce very similar values of the mean square error while an
ensemble size of $N^{2}$ produces a value that is almost half order of
magnitude higher. After 25 iterations the difference between the maximal
theoretical value and the value achieved by CSPSA is between $10^{-1}$ and
$10^{-2}$. Adding 50 more iterations this difference is approximately between
$10^{-2}$ and $10^{-3}$. Let us recall that after 75 iterations the lower
bound of the interquartile range of $\bar{S}(|\psi\rangle\langle\psi|)$ has an
approximate value of $2.12$, so that for 75% of states in the bipartite
Hilbert space we can ascertain its entangled nature and assign an accurate
value of the CHSH function. A further improvement in the accuracy achieved by
CSPSA can be obtained at the expense of a large increase in the number of
iterations, after adding 150 iterations we obtain a new decrease by one order
of magnitude, that is, the mean $\overline{MSE}$ of the square error on the
Hilbert space is approximately in the interval between $10^{-3}$ and
$10^{-4}$.

###  III.2 Unknown mixed states

In the previous section, we have studied the violation of the CHSH inequality
for unknown pure states by means of a CSPSA-driven sequence of local
measurements. Here, we study the case of mixed bipartite states.

We start by reproducing the value of the CHSH function on the set of the
Werner states, which are given by the expression

$\rho_{\lambda}=\lambda|\psi_{s}\rangle\langle\psi_{s}|+\frac{(1-\lambda)}{d}\emph{I},$
(17)

where $|\psi_{s}\rangle$ is the the maximally entangled singlet state defined
as

$|\psi_{s}\rangle=\frac{1}{\sqrt{2}}(|0\rangle|1\rangle-|1\rangle|0\rangle)$
(18)

and I is a 4-dimensional identity operator. This mixture of the singlet state
with white noise is separable if and only if $\lambda\leq 1/3$ and violates
the CHSH inequality if and only if $\lambda>1/\sqrt{2}$. The maximal value of
the CHSH function for a Werner state $\rho_{\lambda}$ is given by

$S(\rho_{\lambda})=2\sqrt{2}\lambda.$ (19)

Figure 12: Mean $\bar{S}(\rho_{\lambda})$ (solid red dots) and median
$\tilde{S}(\rho_{\lambda})$ (blue x’s) as a function of $\lambda$ for Werner
states. Continuous black line depicts the maximal value of the CHSH function
of Eq. (19). Local measurements are simulated with an ensemble size $N=10^{2}$
and 75 iterations are realized. Figure 13: Mean $\bar{S}(\rho_{\lambda})$
(solid red dots) and median $\tilde{S}(\rho_{\lambda})$ (blue x’s) as a
function of $\lambda$ for Werner states. Continuous black line depicts the
maximal value of the CHSH function of Eq. (19). Local measurements are
simulated with an ensemble size $N=10^{4}$ and 75 iterations are realized.

Figure 12 displays the mean $\bar{S}(\rho_{\lambda})$ and median
$\tilde{S}(\rho_{\lambda})$ as a function of $\lambda$ obtained via CSPSA for
an ensemble size $N=10^{2}$ after 75 iterations. With the exception of the
first 5 points, Figure 12 shows a very good agreement between the maximal
value of the CHSH function of Eq. (19) and the value achieved with the help of
CSPSA. Furthermore, mean and median exhibit values that also are very close
and the interquartile (not depicted) range is very narrow. Thus, within the
family of Werner states CSPSA drives the sequence of local measurement bases
very close to the optimal set. An increase in the ensemble size leads to even
better results. This is illustrated in Fig. 13, where local measurements are
simulated with an ensemble $N=10^{4}$. In this case all points are closer to
the maximal value of the CHSH inequality.

Next we proceed with the case of unknown mixed states. We randomly generated a
set of $10^{6}$ two-qubit mixed states. In order to determine whether a mixed
state violates or not the CHSH inequality we employ the $M$ quantity criterion
[49]. A mixed state $\rho$ acting on a Hilbert space
$\mathcal{H}=\mathcal{H}_{2}\otimes\mathcal{H}_{2}$ can be represented in the
form

$\displaystyle\rho$ $\displaystyle=$ $\displaystyle\frac{1}{4}\left(I\otimes
I+\sum_{i=1}^{3}r_{i}\sigma_{i}\otimes
I+I\otimes\sum_{i=1}^{3}s_{i}\sigma_{i}\right.$ (20) $\displaystyle+$
$\displaystyle\left.\sum_{n,m=1}^{3}t_{nm}\sigma_{n}\otimes\sigma_{m}\right),$

where $I$ represents the 2-dimensional identity operator,
$\\{\sigma_{n}\\}_{n=1}^{3}$ are the standard Pauli matrices, and the real
coefficients $r_{i},s_{i}$ and $t_{n,m}$ define the mixed state. The quantity
$M$ is defined by $M(\rho)=u+\tilde{u}$, where $u$ and $\tilde{u}$ denote the
greater positive eigenvalues of the matrix $U_{\rho}:=T_{\rho}^{T}T_{\rho}$
being the coefficients of the matrix $T(\rho)$ given by
$t_{nm}=\mbox{Tr}(\rho\sigma_{n}\otimes\sigma_{m})$. A state $\rho$ violates
the CHSH inequality if and only if the condition $M(\rho)~{}>~{}1$ holds [49].
Employing this criterium, the initial set of $10^{6}$ mixed states was reduced
to a set $\Omega$ containing $8\times 10^{3}$ mixed states with $M(\rho)>1$
that violate the CHSH inequality.

To analyze the values of the CHSH function obtained through CSPSA we use those
obtained through SDP. In the SDP case we need to fix the state that is used in
the maximization. However, let us recall that even when the states are fixed,
the maximization of $S$ remains to be a nonlinear problem. Therefore, to find
the maximum value of $S$ for each state in $\Omega$ we use the see-saw method
[50, 51] to iterate a SDP test [52, 53] where either observable A or B remain
fixed while optimizing in the other variable. The SDP that we solve is the
following

given $\displaystyle\rho_{\Omega},A(z_{a}),A(z_{a}^{\prime}),$ (21)
$\displaystyle\underset{B(z_{b}),B(z_{b}^{\prime})}{\mbox{max}}$
$\displaystyle
S(\rho_{\Omega},A(z_{a}),A(z_{a}^{\prime}),B(z_{b}),B(z_{b}^{\prime})),$ (22)

with the conditions

$\displaystyle|\Psi(z_{b})\rangle\langle\Psi(z_{b})|,|\Psi^{\perp}(z_{b})\rangle\langle\Psi^{\perp}(z_{b})|\geq
0\quad\forall\ z_{b},z_{b}^{\prime},$ (23)
$\displaystyle|\Psi(z_{b})\rangle\langle\Psi(z_{b})|+|\Psi^{\perp}(z_{b})\rangle\langle\Psi^{\perp}(z_{b})|=I\quad\forall\
z_{b},z_{b}^{\prime}.$ (24)

Notice that this SDP takes Alice’s observables $A(z_{a})$ and
$A(z_{a}^{\prime})$ as inputs and for a given mixed state from the $\Omega$
set, it finds Bob’s observables $B(z_{b})$ and $B(z_{b}^{\prime})$ that
maximally violate $S$. Then, we take the observables $B$ outputted by this SDP
as inputs in a new iteration to obtain optimal observables $A$. This procedure
is iterated until some suitable convergence condition is satisfied. We
performed this optimization for every mixed bipartite state in the set
$\Omega$, which allows us to find better lower bounds on $S$, together with
the optimal observables $A$ and $B$.

Figure 14: Mean $\bar{S}_{\Omega}$ (red solid line) and median
$\tilde{S}_{\Omega}$ (blue solid line) obtained via CSPSA on the set $\Omega$
of randomly generated mixed entangled states as a function of the number $k$
of iterations. Mean $\bar{S}_{\Omega}$ (yellow solid line) and median
$\tilde{S}_{\Omega}$ (green solid line) obtained via SDP on the set $\Omega$
of randomly generated mixed entangled states as a function of the number $k$
of iterations. Shaded areas correspond to interquartile range. CSPSA
simulations consider ensemble size $N=10^{4}$.

Figure 14 displays the behavior of the mean $\bar{S}_{\Omega}$, median
$\tilde{S}_{\Omega}$, and interquartile range as functions of the number of
iterations. This figure also displays the values of these quantities obtained
via SDP. As is apparent from this figure, the values of the mean and median
provided via CSPSA are very close and tend to agree with the values delivered
by SDP after tens of iterations. Also, the interquartile ranges tend to
overlap. However, in the case of mixed states the number of iterations needed
to obtain a violation of the CHSH inequality is much greater than in the case
of pure states. This is due to the fact that the mixed states in $\Omega$
typically have small values of the negativity, a well-known measure of
entanglement, and thus, as in the case of weakly-entangled pure states, need
more iterations to reach a violation of the CHSH inequality.

##  IV Conclusions

We have studied the problem of detecting the entanglement of unknown two-qubit
states, mixed or pure, by violating the Clauser-Horne-Shimony-Holt inequality.
Our approach to this problem is based on the maximization of the CHSH function
by means of a stochastic optimization method, the Complex simultaneous
perturbation stochastic approximation. This allows optimizing functions with
unknown parameters, which in our case correspond to the unknown quantum state.
CSPSA employs an iterative rule which requires at each iteration the value of
the target function, that is, the CHSH function, at two different points in
the optimization space. This is formed by vectors on the field of the complex
numbers containing the measurement settings of four observables. The values of
the CHSH function can be experimentally obtained even if the two-qubit state
remains unknown. Thereby, CSPSA generates a sequence of measurement settings
that in mean lead to increasing values of the CHSH inequality.

To analyze the characteristics of the proposed method, we carried out several
numerical experiments. In particular, due to the stochastic nature of CSPSA,
we employ random sampling to obtain estimates of the mean, median, and
interquartile range of the quantities of interest. We first note that for a
fixed unknown state, CSPSA provides very similar values of the mean and median
of the CHSH function and a very narrow interquartile range. This indicates
that CSPSA does not generates outliers, that is, for a given unknown state
different realizations of our method provide very close results. This feature
has been observed for each state in a universe of $5\times 10^{4}$ randomly
generated pure two-qubit states.

The typical behavior of the mean of the CHSH function, as a function of the
number of iterations, corresponds to a rapid increase followed by an
approximately linear asymptotic behavior, which approaches the maximal value
of the CHSH function. Unknown states characterized by the same concurrence
value exhibit a very similar behavior of the CHSH function. However, the rate
of convergence towards the maximum depends on the initial value of the
concurrence. The higher the concurrence value, the fewer iterations are
required to obtain a violation of the CHSH inequality and, consequently,
detect entanglement. For example, states with maximum concurrence need 13
iterations while states with a concurrence of 0.1 need approximately 75
iterations to reach a violation. The number of iterations required to detect
entanglement can be decreased by increasing the size of the ensemble of
identically prepared copies that is employed to estimate the expectation
values entering in the CHSH function. In our simulations, however, the effect
of increasing the ensemble size is more notorious in the case of highly
entangled states. We have studied the mean of the CHSH function on the 2-qubit
Hilbert space. In this case, for an ensemble size of $10^{2}$ the entanglement
of the randomly generated states is detected in mean by violating the CHSH
inequality after 17 iterations, while after 25 iterations 75% of the randomly
generated states violate the CHSH inequality. These figures can be reduced by
increasing the ensemble size. We have also studied the accuracy provided by
our method in the estimation of the maximum value of the CHSH function. As
accuracy metric we have used the mean squared error, which shows that after 25
iterations the difference between the maximal theoretical value and the value
achieved by CSPSA is between $10^{-1}$ and $10^{-2}$. After 75 iterations the
accuracy is approximately between $10^{-2}$ and $10^{-3}$. We have also
considered the case of mixed states. The proposed method is capable of
reproducing the maximal value of the CHSH function for Werner states and for
randomly chosen mixed states.

Therefore, the numerical simulations indicate that the maximization of the
CHSH function through CSPSA leads to the detection of the entanglement of
unknown states, pure or mixed. In mean, 25 iterations detect the entanglement
of 75% of the generated states. Also, it is possible to reach an accurate
value of the maximal violation.

There are some variations of the method here proposed that could reduce the
number of iterations used to detect entanglement. We implement CSPSA
considering the standard choice for the gain coefficients. However, these can
be optimized. This is in general a difficult problem. Nevertheless, some
simple heuristic prescriptions have been discussed in the study of various
proposals of variational quantum eigensolvers [54]. These are based on SPSA, a
version of CSPSA that works on the field of the real numbers. It seems
possible that the SPSA performance-enhancing prescriptions could also be used
to improve the CSPSA convergence rate, which would reduce the number of
iterations required to detect entanglement. The usage of second-order methods
or quantum natural gradient could also speed up the protocol [55, 56, 57, 58,
59]. These employ additional measurements of the objective function to
estimate its Hessian matrix, or fidelity to estimate the metric tensor.
Thereafter, these matrices are used to precondition the gradient in order to
improve the convergence rate, avoiding the need for tuning of some gain
coefficients. Another possibility arises when considering the large amount of
information generated by our method. At each iteration 4 local observables are
measured, which after several iterations provide a considerable amount of
information about the unknown state. Thus, we can obtain an estimate of the
unknown states by means of maximum likelihood [37]. This, together with the
estimate of the optimal measurement settings provided by CSPSA, can be used as
initial guesses in a SDP problem to optimize the CHSH function. The solution
of this problem can be used as the initial guess of the optimal measurement
settings in the next iteration of CSPSA. This procedure does not increases the
amount of measurements to be carried out but the computational cost. Besides,
the use of a priori information can be employed to further increase the CSPSA
convergence rate and achieve entanglement detection with a reduced number of
iterations.

We would like to remark that our approach based on CSPSA can be employed in
other interesting problems. The construction of entanglement witnesses is a
demanding computational task [30, 31], especially if the state is unknown, but
it could be done efficiently with our method. The search for the optimal
measurement settings to violate a multiqubit Bell inequality is challenging
[60, 61, 62]. This is because the dimension scales exponentially with the
number of qubits, so finding the optimal with quantum tomography and SPD is
unfeasible. Our approach could provide an advantage in this problem since its
resource scales with the number of iterations and not with the number of
qubits.

## Acknowledgments

This work was supported by ANID – Millennium Science Initiative Program –
ICN17-012 and by Fondo Nacional de Desarrollo Científico y Tecnológico
(FONDECYT) Grant No 1180558. J. C.-V. was supported by CONICYT-
PCHA/DoctoradoNacional/2018-21181692. J. F. B. acknowledges support from
FONDECYT Grant No 317030. L.P. was supported by ANID-PFCHA/DOCTORADO-BECAS-
CHILE/2019-772200275, the CSIC Interdisciplinary Thematic Platform (PTI+) on
Quantum Technologies (PTI-QTEP+), the CAM/FEDER Project No. S2018/TCS-4342
(QUITEMAD-CM), and the Proyecto Sinérgico CAM 2020 Y2020/TCS-6545 (NanoQuCo-
CM).

## Declarations

The code and the data simulated to generate the figures are available on
reasonable request.

## References

  * [1] J. von Neumann, Mathematische Grundlagen der Quantenmechanich (Springer-Verlag, Berlin, 1932).
  * [2] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Rev. Mod. Phys. 81, 865 (2009).
  * [3] G. Alber, T. Beth, M. Horodecki, P. Horodecki, R. Horodecki, M. Rötteler, H. Weinfurter, R. F. Werner, and A. Zeilinger, Quantum Information: An Introduction to Basic Theoretical Concepts and Experiments (Springer-Verlag, Heidelberg, 2001).
  * [4] A. K. Ekert, Phys. Rev. Lett. 67, 661 (1991).
  * [5] C. H. Bennett, G. Brassard, C. Crépeau, R. Jozsa, A. Peres, and W. K. Wootters, Phys. Rev. Lett. 70, 1895 (1993).
  * [6] D. Boschi, S. Branca, F. De Martini, L. Hardy, and S. Popescu, Phys. Rev. Lett. 80, 1121 (1998)
  * [7] D. J. Wineland, J. J. Bollinger, W. M. Itano, F. L. Moore, and D. J. Heinzen, Phys. Rev. A 46, R6797 (1992).
  * [8] S. F. Huelga, C. Macchiavello, T. Pellizzari, A. K. Ekert, M. B. Plenio, and J. I. Cirac, Phys. Rev. Lett 79, 3865 (1997).
  * [9] V. Giovannetti, S. Lloyd, and L. Maccone, Science 306, 1330 (2004).
  * [10] R. Raussendorf and H. J. Briegel, Phys. Rev. Lett 86, 5188 (2001).
  * [11] R. Jozsa, D. S. Abrams, J. P. Dowling, and C. P. Williams, Phys. Rev. Lett 85, 2010 (2000).
  * [12] Č. Bruckner, N. Paunković, T. Rudolph, and V. Vedral, Int. J. Quantum Inf. 4, 655 (2006).
  * [13] A. Einstein, B. Podolsky, and N. Rosen, Phys. Rev. 47, 777 (1935).
  * [14] J. C. Bell, Physics 1, 195 (1964).
  * [15] J. F. Clauser, M. A. Horne, A. Shimony, and R. A. Holt, Phys. Rev. Lett. 23, 880 (1969).
  * [16] A. Aspect, P. Grangier, and G. Roger, Phys. Rev. Lett. 49, 91 (1982).
  * [17] M. Horodecki, P. Horodecki, and R. Horodecki, Phys. Lett. A 223, 1 (1996).
  * [18] B. M. Terhal, Phys. Lett. A 271, 319 (2000).
  * [19] P. Hyllus, O. Gühne, D. Bruß, and M. Lewenstein, Phys. Rev. A 72, 012321 (2005).
  * [20] Ursin, R., Tiefenbacher, F., Schmitt-Manderbach, T. et al. Entanglement-based quantum communication over 144km. Nature Phys. 3, 481–486 (2007).
  * [21] Steinlechner, F., Ecker, S., Fink, M. et al. Distribution of high-dimensional entanglement via an intra-city free-space link. Nat. Commun. 8, 15971 (2017).
  * [22] Jin, XM., Ren, JG., Yang, B. et al. Experimental free-space quantum teleportation. Nature Photon. 4, 376–381 (2010).
  * [23] Häffner, H., Hänsel, W., Roos, C. et al. Scalable multiparticle entanglement of trapped ions. Nature 438, 643–646 (2005).
  * [24] Gary J Mooney et al, J. Phys. Commun. 5 095004 (2021).
  * [25] Pereira, L., Zambrano, L. and Delgado, A. Scalable estimation of pure multi-qubit states.  npj Quantum Inf 8, 57 (2022).
  * [26] G. M. D’Ariano, M. G. A. Paris, and M. F Sacchi, Adv. Imaging Electron Phys. 128, 205 (2003).
  * [27] M. G. A. Paris and J. Řeháček, Quantum State Estimation (Springer-Verlag, Heidelberg, 2004).
  * [28] T. O. Maciel and R. O. Vianna, Phys. Rev. A 80, 032325 (2009).
  * [29] G. Lima, E. S. Gómez, A. Vargas, R. O. Vianna, and C. Saavedra, Phys. Rev. A 82, 012302 (2010).
  * [30] H. Zhu, Y. S. Teo, and B.-G. Englert, Phys. Rev. A 81, 052339 (2010).
  * [31] J. Dai, Y. L. Len, Y. S. Teo, B.-G. Englert, and L. A. Krivitsky, Phys. Rev. Lett. 113, 170402 (2014).
  * [32] D. Cavalcanti and M. O. Terra-Cunha, Appl. Phys. Lett. 89, 084102 (2006).
  * [33] F. Mintert and A. Buchleitner, Phys. Rev. Lett. 98, 140505 (2007).
  * [34] J. Szangolies, H. Kampermann, and D. Bruß, New J. Phys. 17, 113051 (2015).
  * [35] A.D. Muñoz-Moller, L. Pereira, L. Zambrano, J. Cortés-Vega, and A. Delgado, Phys. Rev. Applied 18, 024048 (2022).
  * [36] A. Utreras-Alarcón, M. Rivera-Tapia, S. Niklitschek, and A. Delgado, Sci. Rep. 9, 16143 (2019).
  * [37] L. Zambrano, L. Pereira, S. Niklitschek, and A. Delgado, Sci. Rep. 10, 12781 (2020).
  * [38] B. Hensen et al., Nature 526, 682 (2015).
  * [39] G. Carvacho et al., Phys. Rev. Lett. 115, 030503 (2015).
  * [40] M. Giustina, M. A. M. Versteegh, S. Wengerowsky, J. Handsteiner, A. Hochrainer, K. Phelan et al., Phys. Rev. Lett. 115, 250401 (2015).
  * [41] C. Abellán et al., Nature 557, 212 (2018).
  * [42] F. Vedovato, C. Agnesi, M. Tomasin, M. Avesani, J. A. Larsson, G. Vallone, and P. Villoresi, Phys. Rev. Lett. 121, 190401 (2018).
  * [43] W. Wirtinger, Math. Ann. 97, 357 (1927).
  * [44] J. C. Spall, IEEE Trans. Autom. Control. 37, 332 (1992).
  * [45] J. C. Spall, Introduction to Stochastic Search and Optimization (Wiley, New Jersey, 2003).
  * [46] C. Ferrie, Phys. Rev. Lett. 113, 190404 (2014).
  * [47] C. Granade, C. Ferrie, and S. T. Flammia, New J. Phys. 19, 113017 (2017).
  * [48] R. J. Chapman, C. Ferrie, and A. Peruzzo, Phys. Rev. Lett. 117, 040402 (2016).
  * [49] R. Horodecki, P. Horodecki, and M. Horodecki, Phys. Lett. A 200, 340-344 (1995).
  * [50] K. F. Pál and T. Vértesi, Phys. Rev. A 82, 022116 (2010).
  * [51] D. Cavalcanti and P. Skrzypczyk, Rep. Prog. Phys. 80, 024001 (2017).
  * [52] L. Vandenberghe and S. Boyd, SIAM Rev. 38, 4 (1996).
  * [53] S. Boyd and L. Vandenberghe, Convex Optimization, (Cambridge University Press, Cambridge, 2004).
  * [54] A. Kandala et al., Nature 549, 242 (2017).
  * [55] J. Spall, Automatic Control, IEEE Transactions on 45, 1839 (2000).
  * [56] C. Wang, arXiv:2012.06952 (2020).
  * [57] J. Stokes, J. Izaac, N. Killoran, and G. Carleo, Quantum 4, 269 (2020).
  * [58] J. Gacon, C. Zoufal, G. Carleo, and S. Woerner, Quantum 5, 567 (2021).
  * [59] J. A. Gidi, B. Candia, A. D. Muñoz-Moller, A. Rojas, L. Pereira, M. Muñoz, L. Zambrano, A. Delgado, arXiv:2203.06044 (2022).
  * [60] G. Svetlichny, Phys. Rev. D 35, 3066 (1987).
  * [61] Daniel Collins, Nicolas Gisin, Sandu Popescu, David Roberts, and Valerio Scarani, Phys. Rev. Lett. 88, 170405 (2002)
  * [62] Michael Seevinck and George Svetlichny, Phys. Rev. Lett. 89, 060401 (2002).