#  ChangeChip: A Reference-Based Unsupervised Change Detection for PCB Defect
Detection

Yehonatan Fridman  Ben-Gurion University of the Negev
Israel Atomic Energy Commission
<EMAIL_ADDRESS>Matan Rusanovsky  Ben-Gurion University of the Negev
Israel Atomic Energy Commission
<EMAIL_ADDRESS>Gal Oren${}^{\textrm{{\char 0\relax}}}$ Technion –
Israel Institute of Technology
Nuclear Research Center – Negev
<EMAIL_ADDRESS>

###### Abstract

The usage of electronic devices increases, and becomes predominant in most
aspects of life. Surface Mount Technology (SMT) is the most common industrial
method for manufacturing electric devices in which electrical components are
mounted directly onto the surface of a Printed Circuit Board (PCB). Although
the expansion of electronic devices affects our lives in a productive way,
failures or defects in the manufacturing procedure of those devices might also
be counterproductive and even harmful in some cases. It is therefore desired
and sometimes crucial to ensure zero-defect quality in electronic devices and
their production. While traditional Image Processing (IP) techniques are not
sufficient to produce a complete solution, other promising methods like Deep
Learning (DL) might also be challenging for PCB inspection, mainly because
such methods require big adequate datasets which are missing, not available or
not updated in the rapidly growing field of PCBs. Thus, PCB inspection is
conventionally performed manually by human experts. Unsupervised Learning (UL)
methods may potentially be suitable for PCB inspection, having learning
capabilities on the one hand, while not relying on large datasets on the
other. In this paper, we introduce ChangeChip, an automated and integrated
change detection system for defect detection in PCBs, from soldering defects
to missing or misaligned electronic elements, based on Computer Vision (CV)
and UL. We achieve good quality defect detection by applying an unsupervised
change detection between images of a golden PCB (reference) and the inspected
PCB under various setting. In this work, we also present CD-PCB, a synthesized
labeled dataset of 20 pairs of PCB images for evaluation of defect detection
algorithms. The sources of ChangeChip, as well as CD-PCB, are available at:
https://github.com/Scientific-Computing-Lab-NRCN/ChangeChip.

###### Index Terms:

Change Detection, Unsupervised Machine Learning, PCA-Kmeans, PCBs, SMT Quality
Control

##  I Introduction

PCB manufacturing and Surface Mount Technology (SMT) are not zero-defect
soldering processes. The ever-changing technology in fabrication, placement,
and soldering of SMT electronic components have caused an increase in PCB
defects both in terms of numbers and types [1]. Defects in SMT processes such
as missing or misaligned components, and other soldering defects like solder
bridging and solder balling, might all cause PCBs to fail. The detection of
such defects is an important part in the control over PCB manufacturing and
SMT processes [2].

Generally, there are several main ways to detect defects in PCBs, each with a
unique disadvantage:

  * •

Manual Visual Inspection (MVI) by experts [3]: As PCBs have evolved to become
more complicated, multi-layered, and assembled with much tinier components,
manual visual inspection has become very tiring, costly, and subjective task.
Therefore, there is a growing need for automation in the PCB inspection field.

  * •

Classic Image Processing (IP) [4, 5, 1, 6, 7]: Seemingly not satisfactory to
produce a complete solution since it cannot understand abstract features from
the data.

  * •

Supervised (e.g. Deep) Learning [8, 9, 7, 10]: Since PCB architectures change
rapidly it is sometimes very hard to attain a proper and general enough
dataset for supervised learning tasks.

  * •

Unsupervised Learning [7, 11]: Usually unsupervised learning models are weaker
than supervised ones since they lack the supervision of the field experts, and
therefore it is harder for them to be optimized against the knowledge of the
experts on the training data.

Therefore, there is a need for algorithms that have learning capabilities that
successfully extract features of interest, and at the same time do not require
large datasets. In this work, we present ChangeChip, an automated change
detection system for defect detection for PCBs based on an unsupervised
learning approach. The objective of ChangeChip is to recommend experts on the
most salient changes between an inspected image and a ’perfect’ reference
image, without any prior knowledge on the PCB design and types of components.
This objective is crucial since the technology is in a constant state of
development toward smaller components which allows integration of more compact
designs at the PCB level, from small single-board microcontrollers to large
motherboards and GPUs.

Input (Reference and Target PCBs)Cropping (e.g. DEXTR)RegistrationExact
Histogram MatchingPCA-Kmeans Change DetectionClasses Analysis and
Outputoptionaloptional Figure 1: The workflow of ChangeChip.

Fig.1 presents the workflow of ChangeChip as a comprehensive integrated
system. ChangeChip performs a preprocessing stage on the input images to make
a correspondence between the two PCBs in cases of imaging discrepancies such
as lighting, imaging view, etc. The preprocessing stage is based on a deep
learning tool named Deep Extreme Cut (DEXTR) [12] that allows to extract the
PCBs and disregard the background, and on IP and CV methods to perform image
registration and histogram equalization. The preprocessing stage can be
interpreted as a pipeline: The PCBs are first extracted in order to ignore the
background if needed, then are registered to align with each other by applying
affine transformation on one of the images, and finally are matched by their
histograms via the Exact Histogram Specification algorithm [13] to reduce
illumination differences between the images. Subsequently, at the heart of the
system, ChangeChip uses PCA-Kmeans on the features of both images to cluster
pixels into classes of change significance. A post-processing stage is then
applied to map the clustering result into a meaningful clear detection of
defects using a Mean Squared Error (MSE) heuristic.

In Section II we mention previous work and discuss the contribution of our
suggested system for quality control of PCBs. In Section III we elaborate on
the algorithmic phases of ChangeChip. Then, in Section IV we present CD-PCB, a
database of PCB images that were created for the change detection problem to
evaluate defect detection algorithms for PCB inspection. We then evaluate the
performance of ChangeChip on that database. In section V we discuss the
conclusions of this work, limitations, and future work.

##  II Previous Work

In the last decades, many techniques and computer algorithms were implemented
in Automated Visual Inspection (AVI) systems for quality assurance for PCB
manufacturing [4, 14, 15, 5, 16, 17, 18, 19, 20]. The main class of such
algorithms is change detection algorithms which use a reference comparison
approach – by detecting changes between a reference image of a zero-defects
already-inspected PCB, with the inspected PCB image. Previous work has been
done on change detection for PCBs, and various state-of-the-art algorithms
were developed using Image Processing (IP), Computer Vision (CV), and Machine
Learning (ML) [1, 6]. The constantly evolving nature and scope of the PCB
assurance problem, pose challenges when using traditional IP and CV
approaches, which operate best in closed systems [7]. Hence, intelligent
algorithms capable of learning (i.e., machine learning and deep learning) are
necessary for achieving and maintaining high accuracy, while avoiding becoming
outdated [9, 19]. Yet, these methods mostly require a large and regularly
updated dataset of images labeled by Subject Matter Experts (SMEs). Although
there are several existing datasets of PCB components, there are few that are
publicly available, such as [10]. Moreover, many of the available datasets are
intended and formatted for a specific study, so it is difficult to combine and
re-purpose them effectively for a different study.

In this work, we present a change detection system for PCB inspection based on
IP and UL. The machine learning algorithm that is used in our system is PCA-
Kmeans, and is a variation of an implementation that was originally suggested
for satellite images in [21]. The usage of UL methods allows us to waive the
need for a labeled updated dataset on the one hand, and on the other hand to
perform a learning procedure that allows extraction of main features from the
data.

##  III The Algorithmic Stages of ChangeChip

###  III-A Preprocessing Stage

One of the key elements of an AVI system is the imaging setting. Naturally, if
the images were captured in identical conditions of lighting, positioning, and
image resolution, then it was easier to detect true changes between them while
avoiding false positives. By performing image registration and histogram
matching, ChangeChip, in effect, relaxes the necessity of both input images to
be captured under the identical conditions. However, insufficient efforts to
maintain identical imaging conditions can still lead to challenges that would
cause the system to fail. For example, capturing a 3D PCB from two different
directions can reveal different sides of the same elements, which can falsely
be detected as changes. Hence, we suggest capturing both input images in the
same imaging environment, maintaining the imaging conditions as identical as
possible. As described in Section IV-A, the CD-PCB database was constructed
under a fixed setting. We note that illumination changes and different object
orientations are still present, despite our efforts to maintain the same
settings.

####  III-A1 Deep Extreme Cut (DEXTR)

[12] is a deep learning tool that extracts objects from an image using 4
extreme points of each object (left-most, right-most, top, bottom pixels).
DEXTR is based on adding an extra channel to the image in the input of a
Convolutional Neural Network (CNN), which contains a Gaussian centered in each
of the extreme points. The CNN learns to transform this information into a
segmentation of an object that matches those extreme points. We use DEXTR in
the preprocessing stage of ChangeChip for PCBs extraction in order to ignore
the background, or even to extract a specific component if comparing only
parts of the PCBs is desired. We stress that ChangeChip is not limited to
DEXTR, and other cropping methods might be suitable too.

####  III-A2 Image Registration

Image registration is the process of estimating an optimal transformation
between two images, in order to bring them into one coordinate system, which
then makes it easier to compare them. We use a feature-based registration,
based on extraction of feature points along with descriptors and matching
them, then using the corresponding pairs to calculate a transformation matrix
that transforms the target image into the same orientation and position as the
reference image. In order to extract and describe feature points in PCBs we
use the Scale-Invariant Feature Transform (SIFT) [22] algorithm which is
rotation, scale, and luminance invariant. Pairs are then matched by finding
for each feature point in one image the best match among the feature points in
the other image (which minimizes the Euclidean distance). We then reject
matches in which the ratio of the nearest neighbor distance to the second-
nearest neighbor distance is greater then 0.8 (a value suggested by Lowe
[23]). According to Lowe, this rejection eliminates $\sim$90% of the false
matches while discarding less than $\sim$5% of the correct matches. Then, a
transformation matrix is determined to map one image to align with the other,
based on these distinct matches. Because the number of matches probably
exceeds the minimum required to define the appropriate transformation, and
they might not follow a single linear model, the iterative method RANSAC [24]
is used for a robust estimation of the transformation matrix.

Reference and Target PCBs
---
|
Extracting the PCBs with DEXTR
|
Matching feature points for registration
The registration and histogram matching result
|
The clustering map and the final detection
|

Figure 2: The workflow of ChangeChip on two optical images of PCBs

Reference and Target PCBs
---
|
Matching feature points for registration
The result of registration
|
The clustering map and the final detection
|

Figure 3: The workflow of ChangeChip on two radiography images of PCBs

####  III-A3 Reducing Illumination Changes

Various technologies are in use in the PCB inspection industry, ranging from
optical photography like still photography and micro-photography, to
radiography using X-rays, gamma rays, etc [4]. The use of radiography enables
the inspection of inner layers of PCBs, which is not possible with optical
imaging. The fact that optical images are easily influenced by illumination
makes it harder to compare them. In order to reduce illumination differences,
we apply Exact Histogram Specification [13] to match the two images’
histograms. This however does not provide a complete solution, especially
since optical images taken from different points of view or with different
lighting settings, might defer more significantly due to reflections on shiny
objects, that cannot be disregarded by applying histogram matching. Yet,
histogram equalization might help in the case of a fixed capture system, with
a slowly fading light source over time. We note that these issues of lighting
do not apply to radiography, and therefore this step of histogram matching is
optional.

###  III-B PCA-Kmeans Change Detection

After performing image registration and histogram matching, we use a variation
of the PCA-Kmeans change detection algorithm that was presented in [21]. This
algorithm requires both images to have the same dimensions
$[1,width]\times[1,height]$. In general, this algorithm uses a descriptor that
extracts for each pixel a set of features, that describe whether that pixel is
a ’change’ from one image to another. Each set of features (i.e. descriptor)
is reduced to a smaller dimension via PCA, and then all these reduced features
are clustered into classes of change degree using Kmeans. In [21], a grayscale
difference window is used to describe the change of each pixel, mainly because
that work attends grayscale images. In our adaptation for colored images, we
use the concatenation of all the RGB-channels difference windows, along with
the grayscale difference window, to preserve and learn meaningful information
about the color features of changes. The PCA-Kmeans algorithm is performed on
small windows around each pixel, and therefore the learning is window-based.
This kind of learning can eliminate further noise that exists in the input
images.

In more details, we define a descriptor for each pixel
$(i,j)\in[1,width]\times[1,height]$ as follows:

$\displaystyle des[i,j]=(des_{RGB}[i,j],des_{gray}[i,j])$ (1)

where

$\displaystyle des_{RGB}[i,j]=(|I_{1}(S)_{R}-I_{2}(S)_{R}|$
$\displaystyle,|I_{1}(S)_{G}-I_{2}(S)_{G}|,$ (2)
$\displaystyle|I_{1}(S)_{B}-I_{2}(S)_{B}|)$ $\displaystyle
des_{gray}[i,j]=(|I_{1}(S)_{gray}-I_{2}(S)_{gray}|)$ (3)

where $S=(i-\frac{h}{2}:i+\frac{h}{2},j-\frac{h}{2}:j+\frac{h}{2})$, a subset
of indices that constitutes a window of $h\times h$ pixels around the pixel
$(i,j)$. $I_{1}$ and $I_{2}$ represent the two images, and
$|I_{1}(S)_{R}-I_{2}(S)_{R}|$ stands for the absolute values of subtraction
between the Red channels of $I_{1}$ and $I_{2}$ restricted to the window $S$.
The same notation is also applied to the Green (G) and Blue (B) channels, and
also to the grayscale value (gray) which is calculated as the luminance of a
pixel according to the following relation: $gray=0.3R+0.59G+0.11B$. For
windows that cross the border, we pad $I_{1}$ and $I_{2}$ with zeros.

Pixels are partitioned into $h\times h$ non-overlapping windows (e.g. $5\times
5$). We denote this set of windows by $W$. For each $window\in W$ we calculate
$des_{RGB}[i,j]$ for the pixel $(i,j)$ in the center of it. $S_{RGB}\leq
h^{2}$ orthonormal eigenvectors are extracted through Principal Component
Analysis (PCA) of the set $\\{des_{RGB}[i,j]\\}_{window\in W}$ to create an
eigenvector space $Eigen_{RGB}$. PCA is also applied on the set
$\\{des_{gray}[i,j]\\}_{window\in W}$ to achieve an eigenvector space
$Eigen_{gray}$, consists of $S_{gray}\leq h^{2}$ vectors.

Each pixel $(i,j)\in[1,width]\times[1,height]$ can be represented by a
($S_{RGB}+S_{gray}$)-dimensional feature vector that is the projection of
$\\{des_{RGB}[i,j]\\}$ onto the generated eigenvector space $Eigen_{RGB}$,
concatenated with the projection of $\\{des_{gray}[i,j]\\}$ onto
$Eigen_{gray}$. According to [21], each $h\times h$ block of the difference
image could contain one of the three different types of data coming from the
difference image: 1) no-change data; 2) change data; and 3) mixture of change
and no-change data. The first two cases happen when the $h\times h$ block is
completely localized on the changed or unchanged areas of the difference
image. The last case happens when the $h\times h$ block is localized on the
boundaries between changed and unchanged regions of the difference image.
Then, if $S_{gray}$ is set to $3$, each type of data could be represented by a
single eigenvector. The authors of [21] report that according to an extensive
set of experiments (on a dataset of grayscale satellite images), when
$3<S_{gray}\leq h^{2}$, there is no significant change in the change detection
performance with respect to that of $S_{gray}=3$. For our application, we also
choose $S_{gray}=3$, and for the RGB features, we choose $S_{RGB}=9$. Based on
experiments that we have conducted on our domain, changing these parameters
leads to no significant change in the results. According to [21], change
detection is achieved by partitioning the feature vector space into two
clusters. We find that in the case of PCBs images, higher resolution
clustering is needed to successfully separate the changes from the other
pixels. The reason for that is the high resolution usually PCBs images are
taken in, and the bigger intricacy PCBs images have compared to satellite
images. The number of classes, $n$, is a parameter of our application, usually
is set around 10 and might differ from system to system to achieve optimized
results.

###  III-C Post-processing Stage: Classes analysis using an MSE heuristic

The previous stage clusters all the pixels into classes
$\\{C_{1},C_{2},...C_{n}\\}$. In the post-processing stage, we characterize
each class according to a certain heuristic, in order to determine which
classes represent a ’change’ in a more meaningful fashion. We use the Mean
Squared Error (MSE) measure, defined as follow for each class $C_{i}$:

$MSE(C_{i})=\frac{1}{3|C_{i}|}\sum_{(i,j)\in
C{i}}\sum_{c\in\\{R,G,B\\}}(I_{1}(i,j)_{c}-I_{2}(i,j)_{c})^{2}$

The classes are sorted according to the MSE scores, and each class receives a
color according to its index in that order, from blue for the class with the
minimum MSE score, to red for the class with the maximum MSE score. The
assignment of colors to the classes yields a meaningful heat map, which
indicates the intensity of the changes. However, in order to produce a map
with only the suspected changes tagged, some of the lower MSE score classes
have to be discarded. For this purpose, we suggest clustering the MSE scores
using the Density-Based Spatial Clustering of Applications with Noise (DBSCAN)
algorithm [25], followed by discarding the lowest MSE score class. The
advantage of DBSCAN over Kmeans, for example, is that DBSCAN does not require
to specify the number of clusters. Instead, it requires setting the density
parameter (eps), which according to our tests produces better results.
Although eps may be tuned and fixed based on some learning samples, it still
might be modified from application to application for optimized results.

##  IV Evaluation

###  IV-A Change-Detection-PCB (CD-PCB) Dataset

In order to evaluate the performance of ChangeChip, and as a proposed dataset
for future works, we present Change-Detection-PCB (CD-PCB) dataset, which is
composed of 20 pairs of PCB images with annotated changes between them. CD-PCB
includes images of PCBs with various types of challenging defects. Two of the
pairs were captured via radiographic technology, and the other 18 pairs were
synthesized and captured under an established environment that maintains the
same imaging conditions (although a better setting can be achieved [26]). To
our knowledge, such a database does not exist in this form. The imaging system
that was used to form CD-PCB, consists of a camera, tripod, 2 studio flashes,
and a fixed desk. A sketch of the system is presented in Fig. 4.

Figure 4: The capture system model used to construct CD-PCB

###  IV-B Evaluation Method

We use Precision and Recall to evaluate ChangeChip. Recall is the fraction of
true ’change’ pixels that were classified as such by the algorithm, out of all
the true ’change’ pixels. Precision is the fraction of true ’change’ pixels
among all pixels that were classified (correctly and incorrectly) as changes
by the algorithm. High Recall is vital for an algorithm to be considered as
providing high assurance. Since every defect on a PCB can be crucial to its
functioning, high assurance of defect detection is desired. On the other hand,
high Precision is also desired in order to focus on the true defects. High
rates of false positives detection point to the weakness of the algorithm to
characterize and extract defects correctly.

###  IV-C Results

In this section, we evaluate ChangeChip on CD-PCB. The parameters that were
used to conduct this experiment are $h=5$ (windows of size $5\times\ 5$) and
$n=16$ (number of classes to cluster the descriptors). For the optical images
we use, $eps=0.02$ (DBSCAN’s density parameter for clustering the $n$ MSE
values), and for the radiographic images we set $eps=0.05$, since the noises
are less concentrated around the lowest MSE score class. We note that this set
of parameters might change from application to application for optimized
results.
The recall of ChangeChip on CD-PCB is $0.87$ and the precision is $0.8$.
Although it seems that the precision is low, we note that the true changes are
tiny and therefore, in effect, the number of false positives is very small in
relation to the whole image. Fig. 5 presents couples of images from CD-PCB in
the first two columns, with ground truth changes in the third column, and the
result of ChangeChip on those images in the fourth column. The quantitative
results represent the performance of the algorithm on the pixel level, while
the qualitative results intuitively refer to the object level. It can be seen
that ChangeChip marks all the changes that were present in the ground truth,
without any significant noise except for couple number 14, which we discuss in
Chapter V. We conclude that ChangeChip excels qualitatively to highlight all
changes to the human expert that operates the system.

##  V Conclusions, Limitations and Future Work

In this work, we present ChangeChip, a comprehensive system for detecting
defects on PCBs by applying change detection. The system suggests a
comprehensive procedure to compare a golden PCB and an inspected PCB. We note
that ChangeChip is not limited to PCBs only and may be found suitable to other
domains as well. As a limitation of this work, we stress that the system is
targeted to find changes between the two images, while changes, in fact, might
include noise in addition to defects. To illustrate that point, refer to row
14 in table 5. In this example, certain changes are characterized as noise by
our experts, and thus are considered to be false positives. To cope with this
issue, one can specifically characterize the attributes of noise in his system
and discard pixels that appear to be noise in some probability. Moreover,
manual labeling of false positives over time can be used to build a classified
data-set of ’noise’ and ’defect’ pixels (or windows of pixels) in one’s
system, then use Supervised (e.g. Deep) Learning techniques to automatically
classify new pixels. Another direction for decreasing the effect of noise
might be by choosing more complicated descriptors that describe the pixels,
with higher tolerance to noise, and more sensitivity to defects. These
descriptors can include for example color, edges, and texture features. As
future work, we also suggest increasing the support of the system in
radiography images, by analyzing multiple images of adjacent layers, in order
to strengthen the assurance in a significant defect – if it appears as a
’change’ in multiple results. This future work will amplify the support of
this work to 3D scanning of multi-layered PCBs.

## Acknowledgments

We would like to thank Amir Ellenbogen and Gabi Dadush for a fruitful
collaboration with our team in this project. We also thank Ori Stern and Ami
Moskowitz for helping in the creation of CD-PCB. This work was supported by
the Lynn and William Frankel Center for Computer Science. Computational
support was provided by the NegevHPC project [27].

|  Reference PCB |  Target PCB |  Ground Truth |  ChangeChip Results
---|---|---|---|---
1 |  |  |  |
2 |  |  |  |
3 |  |  |  |
4 |  |  |  |
5 |  |  |  |
6 |  |  |  |
7 |  |  |  |
8 |  |  |  |
---|---|---|---|---
9 |  |  |  |
10 |  |  |  |
11 |  |  |  |
12 |  |  |  |
13 |  |  |  |
14 |  |  |  |
Figure 5: Example of pairs from CD-PCB, the ground truth changes and
ChangeChip results according to the parameters described in the Results
section. The red circles are for easy identification by the reader.

## References

  * [1] Kenneth Sundaraj.  Pcb inspection for missing or misaligned components using background subtraction.  WSEAS transactions on information science and applications, 6(5):778–787, 2009.
  * [2] Hasit Desai.  Defect free assembly of smt devices.  In Proceedings: Electrical Insulation Conference and Electrical Manufacturing and Coil Winding Conference, pages 677–682. IEEE, 1997.
  * [3] Paul HP Yeow and Rabindra Nath Sen.  Ergonomics improvements of the visual inspection process in a printed circuit assembly factory.  International journal of occupational safety and ergonomics, 10(4):369–385, 2004.
  * [4] Madhav Moganti, Fikret Ercal, Cihan H Dagli, and Shou Tsunekawa.  Automatic pcb inspection algorithms: a survey.  Computer vision and image understanding, 63(2):287–313, 1996.
  * [5] Jianjie Ma.  Defect detection and recognition of bare pcb based on computer vision.  In 2017 36th Chinese Control Conference (CCC), pages 11023–11028. IEEE, 2017.
  * [6] Abdel-Aziz IM Hassanin, Fathi E Abd El-Samie, and Ghada M El Banby.  A real-time approach for automatic defect detection from pcbs based on surf features and morphological operations.  Multimedia Tools and Applications, 78(24):34437–34457, 2019.
  * [7] Ulbert J Botero, Ronald Wilson, Hangwei Lu, Mir Tanjidur Rahman, Mukhil A Mallaiyan, Fatemeh Ganji, Navid Asadizanjani, Mark M Tehranipoor, Damon L Woodard, and Domenic Forte.  Hardware trust and assurance through reverse engineering: A survey and outlook from image analysis and machine learning perspectives.  arXiv preprint arXiv:2002.04210, 2020.
  * [8] Weibo Huang and Peng Wei.  A pcb dataset for defects detection and classification.  arXiv preprint arXiv:1901.08204, 2019.
  * [9] Mukhil Azhagan Mallaiyan Sathiaseelan, Olivia P Paradis, Shayan Taheri, and Navid Asadizanjani.  Why is deep learning challenging for printed circuit board (pcb) component recognition and how can we address it?  Cryptography, 5(1):9, 2021.
  * [10] Hangwei Lu, Dhwani Mehta, Olivia P Paradis, Navid Asadizanjani, Mark Tehranipoor, and Damon L Woodard.  Fics-pcb: A multi-modal image dataset for automated printed circuit board visual inspection.  IACR Cryptol. ePrint Arch., 2020:366, 2020.
  * [11] Ihar Volkau, Abdul Mujeeb, Dai Wenting, Erdt Marius, and Sourin Alexei.  Detection defect in printed circuit boards using unsupervised feature extraction upon transfer learning.  In 2019 International Conference on Cyberworlds (CW), pages 101–108. IEEE, 2019.
  * [12] Kevis-Kokitsi Maninis, Sergi Caelles, Jordi Pont-Tuset, and Luc Van Gool.  Deep extreme cut: From extreme points to object segmentation.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 616–625, 2018.
  * [13] Philippe Bolon Coltuc, Dinu and J-M. Chassery.  Exact histogram specification.  IEEE Transactions on Image Processing, 15.5:1143–1152, 2006.
  * [14] Huihui Wu, Guanglin Feng, Huiwen Li, and Xianrong Zeng.  Automated visual inspection of surface mounted chip components.  In 2010 IEEE International Conference on Mechatronics and Automation, pages 1789–1794. IEEE, 2010.
  * [15] EK Teoh, DP Mital, BW Lee, and LK Wee.  Automated visual inspection of surface mount pcbs.  In [Proceedings] IECON’90: 16th Annual Conference of IEEE Industrial Electronics Society, pages 576–580. IEEE, 1990.
  * [16] G Acciani, G Brunetti, E Chiarantoni, and G Fornarelli.  An automatic method to detect missing components in manufactured products.  In Proceedings. 2005 IEEE International Joint Conference on Neural Networks, 2005., volume 4, pages 2324–2329. IEEE, 2005.
  * [17] Ulbert J Botero, Fatemeh Ganji, Navid Asadizanjani, Damon L Woodard, and Domenic Forte.  Semi-supervised automated layer identification of x-ray tomography imaged pcbs.  In 2020 IEEE Physical Assurance and Inspection of Electronics (PAINE), pages 1–6. IEEE, 2020.
  * [18] Nathan T Jessurun, Olivia P Paradis, Mark Tehranipoor, and Navid Asadizanjani.  Shade: Automated refinement of pcb component estimates using detected shadows.  In 2020 IEEE Physical Assurance and Inspection of Electronics (PAINE), pages 1–6. IEEE, 2020.
  * [19] Eyal Weiss.  Electronic component solderability assessment algorithm by deep external visual inspection.  In 2020 IEEE Physical Assurance and Inspection of Electronics (PAINE), pages 1–6. IEEE, 2020.
  * [20] RD Taylor.  Determination of hole location on printed circuit boards.  Technical report, 1972.
  * [21] Turgay Celik.  Unsupervised change detection in satellite images using principal component analysis and $k$-means clustering.  IEEE Geoscience and Remote Sensing Letters, 6(4):772–776, 2009\.
  * [22] David G. Lowe.  Object recognition from local scale-invariant features.  In Proceedings of the International Conference on Computer Vision, volume 2, pages 1150–1157, 1999.
  * [23] David G. Lowe.  Distinctive image features from scale-invariant keypoints.  In International Journal of Computer Vision, volume 60 (2), pages 91–110, 2004.
  * [24] Konstantinos G Derpanis.  Overview of the ransac algorithm.  Image Rochester NY, 4(1):2–3, 2010.
  * [25] Kamran Khan, Saif Ur Rehman, Kamran Aziz, Simon Fong, and Sababady Sarasvady.  Dbscan: Past, present and future.  In The fifth international conference on the applications of digital information and web technologies (ICADIWT 2014), pages 232–238. IEEE, 2014.
  * [26] Gal Oren, Lior Amar, David Levy-Hevroni, and Guy Malamud.  The looking-glass system: A unidirectional network for secure data transfer using an optic medium.  In International Conference on Future Network Systems and Security, pages 81–97. Springer, 2016.
  * [27] NegevHPC Project.  www.negevhpc.com.  [Online].


# Race for Quantum Advantage using Random Circuit Sampling

Sangchul Oh  Department of Chemistry, Department of Physics and Astronomy, and
Purdue Quantum Science and Engineering Institute, Purdue University, West
Lafayette, IN, USA  Sabre Kais<EMAIL_ADDRESS>Department of Chemistry,
Department of Physics and Astronomy, and Purdue Quantum Science and
Engineering Institute, Purdue University, West Lafayette, IN, USA

###### Abstract

Random circuit sampling, the task to sample bit strings from a random unitary
operator, has been performed to demonstrate quantum advantage on the Sycamore
quantum processor with 53 qubits and on the Zuchongzhi quantum processor with
56 and 61 qubits. Recently, it has been claimed that classical computers using
tensor network simulation could catch on current noisy quantum processors for
random circuit sampling. While the linear cross entropy benchmark fidelity is
used to certify all these claims, it may not capture in detail statistical
properties of outputs. Here, we compare the bit strings sampled from classical
computers using tensor network simulation by Pan et al. [Phys. Rev. Lett. 129,
090502 (2022)] and by Kalachev et al. [arXiv:2112.15083 (2021)] and from the
Sycamore and Zuchongzhi quantum processors. It is shown that all Kalachev et
al.’s samples pass the NIST random number tests. The heat maps of bit strings
show that Pan et al.’s and Kalachev et al.’s samples are quite different from
the Sycamore or Zuzhongzhi samples. The analysis with the Marchenko-Pastur
distribution and the Wasssertein distances demonstrates that Kalachev et al.’s
samples are statistically close to the Sycamore samples than Pan et al.’s
while the three datasets have similar values of the linear cross entropy
fidelity. Our finding implies that further study is needed to certify or beat
the claims of quantum advantage for random circuit sampling.

##  I Introduction

The outperformance of quantum computers over classical computers for some
computational tasks, called quantum advantage or quantum supremacy, is one of
main subjects of quantum computing research [1, 2, 3]. Random quantum circuit
sampling [4, 5, 6] and Boson sampling [7, 8] are considered good computational
tasks for demonstrating quantum advantage with current noisy quantum
computers. In 2019, the Google team claimed the achievement of quantum
advantage for random circuit sampling [9]. The Google Sycamore quantum
processor with 53 qubits took 200 seconds to sample ten millions of bit
strings from a random quantum circuit while a supercomputer at that time was
estimated to take 10,000 years using the Schrödinger-Feynman algorithm to
simulate the random quantum circuits [9]. Subsequently, the USTC team
implemented random circuit sampling on the Zuchongzhi quantum processor with
56 and 61 qubits [10, 11]. Also, quantum advantage for Boson sampling with
optical qubits has been reported [12, 13].

On the opposite side, to close the quantum advantage gap, the tensor network
simulation of random quantum circuits has been performed on classical
computers [14, 15, 16, 17, 18, 19, 20]. For example, Kalachev et al. [17]
reported that a tensor-contraction algorithm on a graphic processing unit took
14 days to sample millions of bit-string from the same random circuit that was
implemented on the Sycamore processor with 53 qubits and 20 cycles. Similarly,
Pan et al. [19, 20] showed that graphic processing units using a tensor
network took 1 day to sample millions of bit-strings from the same random
quantum circuit performed on the Sycamore processor with 53 qubits and 20
cycles. Thus, it is claimed that classical computers could catch up the
current noisy quantum processors [21].

All of above arguments on quantum advantage for random circuit sampling on
quantum or classical processors are based on the two metrics: computational
time and the linear cross-entropy benchmark fidelity (linear XEB). The first
metric is just run time to accomplish a computational task and may be a
primary measure to compare the performance of computers. The second metric is
employed to ensure that the output of noisy quantum processors is close to an
ideal result. The linear XEB is zero if bit-strings are sampled from a
classical uniform probability distribution. While a non-zero value of the
linear XEB was used to support the claims of quantum advantage for random
circuit sampling, the liner XEB has some limitations in addition to the
intrinsic scalability issue [6]. The linear XEB could be spoofed without fully
simulating a random quantum circuit [22, 23]. Our previous works [24, 25]
demonstrated that the linear XEB cannot capture the statistical properties of
samples in the sense that the Sycamore bit-strings are statistically different
from the Zuchongzhi bit-strings while both processors show similar linear XEB
values as a function of the number of qubits or the number of cycles.

While the benchmark and certification of quantum processors are important to
ensure the faithful calculation with a quantum processor, they are not fully
explored [26]. To validate the claim of quantum advantage and to verify
quantum circuits, it would be better to compare directly the statistics of
outcomes generated by quantum or classical processors. Quantum advantage
regime means the verification with classical processors is very inefficient In
this paper, we compare bit strings sampled from the random circuits on the
Sycamore processor [27] and those calculated from tensor network simulation of
the Sycamore random circuit on classical computers by Kalachev et al. [28] and
Pan et al.’ [29]. The paper is organized as follows. In Sec. II, we will
perform the NIST random number tests on bit strings and plot the heat maps of
bit strings. Also, we will calculate the Marchenko-Pastur distances and the
Wasserstein distances between bit strings. These statistical analysis will
show the samples from classical simulation are quite different from the
samples from the Sycamore quantum processor. In Sec. III, we will summarize
our results and discuss the recent claim that quantum advantage for random
circuit sampling is faded by classical computers.

##  II Sampling from Random Circuits

###  II.1 Random Quantum Circuit Sampling

Let us start with a brief summary of random quantum circuit sampling. The task
of random quantum circuit sampling is to sample $n$-bit strings, $x\equiv
a_{1}a_{2}\cdots a_{n}$, with $a_{i}=0,1$ and $i=1,\dots,n$ from a probability
$p(x)=|\bra{x}U\ket{0^{\otimes n}}|^{2}$ defined by a random unitary operator
$U$. Here, $\ket{x}=\ket{a_{1}a_{2}\cdots a_{n}}$ stands for a computational
basis of $n$-qubits and $\ket{0^{\otimes n}}=\ket{00\cdots 0}$ for an initial
state. In an ideal situation, a random unitary operator $U$ is sampled
uniformly and randomly from the unitary group $\mathbb{U}(N)$ with $N=2^{n}$.
A unique probability measure that is invariant under group multiplication is
known the Haar measure. A set of random unitary matrices drawn from
$\mathbb{U}(N)$ according to the Haar measure is called the circular unitary
ensemble, which was introduced by Dyson [30].

A random unitary matrix with respect to the Haar measure can be drawn using
the Euler angle method [31, 32, 33] or the QR decomposition algorithm [34].
Hurwitz [31] showed that a unitary matrix can be constructed through the
parameterization of $N^{2}$ elements with Euler angles. A random unitary
matrix $U$ can be generated through the QR decomposition of a matrix $Z$ with
complex normal random entries, $Z=QR$ where $Q$ is a unitary matrix and $R$ is
a upper-triangle matrix. Then, $U=Q\Lambda$ is a Haar-measure random unitary
where $\Lambda$ is the diagonal matrix with diagonal entries
$\Lambda_{ii}=R_{ii}/|R_{ii}|$ and $R_{ii}$ are the diagonal entries of the
upper triangle matrix $R$ [34].

Emerson et al. [35] proposed the implementation of pseudo-random unitary
operators on a quantum computer by applying single-qubit gates for each qubits
and the simultaneous Ising interaction for all qubits. The Google and USTC
experiments implemented a pseudo-random unitary operator that is composed of
$m$ cycles of the sequence of $n$ single-qubit gates selected randomly from
$\\{\sqrt{X},\sqrt{Y},\sqrt{W}\\}$ and fixed two-qubit gates for according to
tiling patterns.

A random unitary operator $U$ transforms the initial state $\ket{0^{\otimes
n}}$ into a random quantum state

$\displaystyle\ket{\psi}=U\ket{0^{\otimes
n}}=\sum_{x=0}^{2^{n}-1}\sqrt{p_{x}}\,e^{i\theta_{x}}\ket{x}\,,$ (1)

where $p_{x}=p(x)$. It is known that the probability distribution $Pr$ of
$p_{x}$ with $x=0,\dots,2^{n}-1$ is given by $Pr(p)=(N-1)(1-p)^{N-2}$ with
$N=2^{n}$, the eigenvector distribution of the circular unitary ensemble [30].
For $N\gg 1$, it becomes $Pr(p)\approx N\exp(-Np)$ and is known as the
$\chi^{2}$ distribution with 2 degrees of freedom, or the exponential
distribution [36, 37, 38]. The output of random circuit sampling is a set of
$M$ bit-strings, ${\cal D}=\\{x_{1},x_{2},\dots,x_{M}\\}$, generated by $M$
measurements of $\ket{\psi}$ in computational basis. Note that a current
quantum processor could generate only a bit-strings $x$, but not an amplitude
$\sqrt{p_{x}}\,e^{i\theta_{x}}$, while a classical computer could calculate
both quantities. The phases $\theta_{x}$ are distributed uniformly in the
range $[-\pi,\pi]$.

Today quantum processors are noisy and imperfect so the output of quantum
computation could be deviated from an ideal solution. It is important to
verify whether quantum computation is performed faithfully [26]. Gilchrist et
al. [39] suggested a golden standard for measures of the distance between
ideal and noisy quantum processors. For sampling computation where the task is
to sample outcomes $x$ from an ideal probability distribution $p(x)$,
Gilchrist et al. proposed the the Kolmogorov distance
$D(p,q)\equiv\sum_{x}|p(x)-q(x)|/2$ and the Bhattacharya overlap
$F(p,q)\equiv\sum_{x}\sqrt{p(x)q(x)}$ where $p(x)$ is an ideal distribution
and $q(x)$ is a real probability distribution generated by a quantum
processor. On the other hand, Arute et al. introduced the linear XEB

$\displaystyle F_{\rm XEB}=\frac{2^{n}}{M}\sum_{i=1}^{M}p_{\rm
ideal}(x_{i})-1\,,$ (2)

where the set of bit strings, ${\cal D}=\\{x_{1},x_{2},\dots,x_{M}\\}$, is the
measurement outcomes of a noisy quantum processor and $p_{\rm
ideal}(x)=|\bra{x}U\ket{0}^{\otimes n}|^{2}$ is an ideal probability of a
quantum circuit $U$. For a classical uniform random distribution $p_{\rm
ideal}=1/2^{n}$, the linear XEB becomes zero. Basically, the linear XEB is
derived from the cross entropy between the ideal distribution $p_{\rm
ideal}(x)$ and the real distribution $p_{\rm real}(x)$ or the Kullback-Leibler
divergence. However, the Kullback-Leibler divergence is not a true metric
because it is not symmetric and does not satisfy the triangular inequality.

One of the main challenges in calculating the linear XEB or the Kolmogorov
distance is how to obtain the ideal probability distribution $p_{\rm
ideal}(x)$ or the real probability distribution $p_{\rm real}(x)$. For
sampling calculation, a quantum computer produces bit-strings rather than
$p(x)$. In quantum advantage regime, a classical computer is hard or
inefficient to calculate $p(x)$ from quantum simulation. Also, it is
impossible to estimate the empirical probability distribution from the
measurement data of a quantum processor since the number of measurements $M$
is much less than $2^{n}$. In the Google experiment with the Sycamore quantum
processor, the number of measurements, $M$, are order of $10^{7}$ but the
dimension of the Hilbert space of 53 qubits is $2^{53}\sim 9\times 10^{15}$.
Thus, to verify the performance of quantum computers in quantum advantage
regime, it would be better to compare directly two data sets, ${\cal
D}_{1}=\\{x_{1},x_{2},\dots,x_{M}\\}$ and ${\cal
D}_{2}=\\{y_{1},y_{2},\dots,y_{M}\\}$ without calculating the probability
distributions $p(x)$ or $q(x)$. For example, since classical uniform random
bit strings give rise to zero of the linear XEB, they could be one reference
point.

###  II.2 Statistics of Samples of Random Circuits

####  II.2.1 Datasets of random circuit sampling

We compare the 4 datasets of bit strings sampled from the Sycamore random
circuits. The first is the Google dataset used in quantum supremacy
experiments with the Sycamore quantum processor [27]. In order to compare with
other dataset, only sub-dataset for 53 qubits and cycle $m=12,14,16,18,20$ are
considered. Each data file is identified using five parameters: $n$ is the
number of qubits, $m$ the number of cycles $s$ seed for the pseudo-random
number generator, $e$number of elided gates, and $p$ sequence of coupler
activation patterns. For example, measurement-n53-m20-s0-e0-pABCDCDAB.txt
stand for the measurement data file for 53 qubits, 20 cycles and ABCDCDAB
activation pattern.

The second dataset is Kalachev et al.’s [17]. They produced one million bit-
strings for the Sycamore random circuit up to 20 cycles using the multi-tensor
contraction algorithms and a modified frugal rejection sampling method. To
produce $k$ bit-strings, they first calculated $2k$ random batches and then
selected $k$ bit-strings out of $2k$ batches using the frugal rejection
sampling. We obtained Kalachev et al.’s data from [28]. The data set is
composed of 5 text files of one million bit strings for the Sycamore random
circuit with $m=12,14,16,18$, and $20$ cycles (for example,
samples-m20-f0-002.txt). For cycle $m=12,14,16,18$, the target fidelity is
0.02 and for $m=20$, the target fidelity is 0.02. Also, one text file for
spoofing the linear XEB (spoofing-m20.txt) is provided.

The third dataset is generated by Pan et al [29] using the tensor network
simulation. However, only the spoofing dataset with about one million bit
strings, samples-metropolis.txt, is available. Finally, we prepare the file of
bit strings sampled from a classical uniform random distribution
(n50-classical.txt). Note that the Zuchongzhi data for $n=53$ qubit are not
available [10] so we compare only the Sycamore dataset with Kalachev et al.’s
data and Pan et al. data.

####  II.2.2 NIST random number tests.

If one take a close look at bit strings sampled from a random quantum circuit,
bit 0 and 1 seem to appear randomly. Note that we are focusing on randomness
of bits of a sample, not bit-strings. As discussed, the probability of finding
a bit string $x$ is given by $p(x)=|\bra{x}U\ket{0^{\otimes n}}|^{2}$ and is
not uniformly random. So some bit-strings are sampled more likely. While it is
not clear whether a random circuit produces bit 0 or 1 uniformly at random,
our previous studies [24, 25] showed that some of the Zuchongzhi samples pass
the NIST random number tests [40]. Recently, Shi et al. [41] demonstrated a
possibility of using the Boson sampling method to generate true random bits.
They showed that bit-strings generated from the Boson sampling on optical
qubits pass the NIST random number tests while more simple ways of generating
true quantum random numbers exist [42]. Here, we perform the NIST random
number tests [40, 43] on bit strings sampled form the Sycamore random circuit
using tensor network simulation: Kalachev et al.’s five samples and Panet
al.’s spoofing sample. As expected, Pan et al’s spoofing sample does not pass
the NIST random number tests because first several bits of 56-bit strings are
fixed to zeros or ones. We find that all five samples of Kalachev et al. pass
the NIST random number tests. The details of the NIST random number tests for
the six samples are shown in Supplemental Material [44]. However, we leave an
open question whether bits of bit-strings sampled from a random circuit should
be uniformly at random.

Figure 1:  (First row) For 2-dimensional arrays of the Sycamore quantum
processor with 53 qubits, the averages of finding outcome 1 of each qubit over
one million measurements are plotted for (a) the classical random bit strings,
(b) the Sycamore sample (measurement-n53-m20-s0-e0-pABCDCDAB.txt), (c)
Kalachev et al. sample (samples-m20-f0-002.txt), and (d) Pan sample (sample-
metropolis.txt). (Second row) The heat maps are plotted for (e) classical
random bits (f) the Sycamore sample, (f) Kalachev et al.’s sample, (g) Pan et
al.’s sample. The average of finding bit 1 over the sample is denoted by
$p_{1}$.

####  II.2.3 Heat maps of bit strings.

A heat map is a good way of visualizing the statistical properties of outcomes
of each qubits. The averages of outcomes at each qubit arranged in the
2-dimensional array of the Sycamore processor are plotted for the classical
random bits in Fig. 1 (a), the Sycamore sample in Fig. 1 (b), the Kalachev
sample in Fig. 1 (c), and the Pan et al.’s sample in Fig. 1 (d), Kalachev et
al.’s sample and Pan et al.’s sample, both of which were obtained using the
tensor network simulation, look like only that a random quantum circuit acts
on only some parts of 53 qubits. Especially, Pan et al’s spoofing sample is
quite different from the other three samples. As the animations in
Supplemental Material show [44], for Pan et al’ spoofing sample, some qubits
are turned always off (and turned on later).

To see the statistical difference among samples in detail, we slice a
$10^{6}\times 53$ bit array, ${\cal D}=(x_{1},x_{2},\dots,x_{M})^{t}$, into
$53\times 53$ square arrays $\\{A^{(i)}\\}$ with $i=1,\dots,L$ and
$L\equiv{\rm int}(M/53)$. The heat maps of the averages of $A^{(i)}$,
$(1/L)\sum_{i}A^{(i)}$, for the classical random bit sample, the Sycamore
sample, the Kalachev et al.’ sample, and the Pan et al.’s sample are plotted
in Figs. 1 (e), (f), (g), and (h), respectively. The heat map of the classical
random bits has no stripe patterns as expected. The heat map of Pan et al.
sample exhibits also the stripe patterns, but is quite different from those of
the Sycamore sample of Kalachev et al.’s. It is interesting to see the heat
map of Kalachev et al’s samples, Fig. 1 (f) shows stripe patterns similar to
that of the Sycamore sample, Fig. 1 (e) while the former passes the NIST
random number tests. One may expect that bit strings which pass the NIST
random number tests would do not show any patterns. However, Kalachev et al’s
samples here and Zuchongzhi samples in our previous work [25] may break a
common belief on random bit or random numbers.

We count the number of outcome $1$ of a sample composed of one million bit-
strings, denoted by $p_{1}$. As shown in Figs. 1 (e), (f), (g) and (h), the
averages of finding outcome 1 are estimated as $p_{1}=0.50001$ for the
classical random bit sample, $p_{1}=0.48360$ for the Sycamore sample,
$p_{1}=0.50004$ for Kalachev et al’s sample, and $p_{1}=0.545000$ for Pan et
al’s sample. Note that the averages of finding outcome 1 for Zuchongzhi
samples are almost 0.5. the Sycamore samples have $p_{1}$ less than 0.5
because of the readout errors [45].

####  II.2.4 Marchenko-Pastur distances and Wasserstein distances

One expects that two classical processors will produce the same output or
statistically similar outputs. Like that, two quantum processors implement
similar quantum gates will produce similar outcomes [46]. Suppose a unitary
operator $V$ of an ideal processor is approximate to a unitary operators $U$
of a real processor. The error between $U$ and $V$ is defined by
$E(U,V)\equiv\max_{\ket{\psi}}||(U-V)\ket{\psi}||$. If the difference in two
gates is small, then measurement statistics of $U\ket{\psi}$ is approximated
by $V\ket{\psi}$. This is written as the inequality [46]

$\displaystyle|p_{U}-p_{V}|\leq 2E(U,V)\,.$ (3)

If the error $E(U,V)$ is small, the difference in probabilities between
measurement outcomes, $p_{U}$ and $p_{V}$, is also small.

Figure 2:  (a) Marchenko-Pastur distribution of eigenvalues of
$(1/k)X^{t}\cdot X$, Eq. (5), for the Sycamore sample with $n=53$ and $m=20$
(red color) and for the Kalachev et al.’s sample with $n=53$ and $m=20$
(orange color). The inset shows the outlier peaks around $n/4$. (b) Marchenko-
Pastur distribution of eigenvalues of $(1/k)X^{t}\cdot X$ for Pan et al.’s
spoofing sample (samples-metropolis.txt) and Kalachev et al.’s spoofing sample
(spoofing-m20.txt). The inset shows the outliers. (c) Marchenko-Pastur
distances (outlier distances from $n/4$) as a function of the number of cycles
$m$ for the Sycamore samples and Kalachev et al.’s samples with $n=53$ qubits,
and for the Zuchongzhi samples with $n=56$.

If two processors, no matter processor are quantum or classical, implement the
similar unitary operators for random circuit sampling, statistical distances
between probabilities should be close to each other. In practice, the two
probability distributions, $p_{U}(x)$ and $p_{V}(y)$, may be unknown or
unavailable, but only outcomes ${\cal D}_{U}=\\{x_{1},x_{2},\dots,x_{M}\\}$
and ${\cal D}_{V}=\\{y_{1},y_{2},\dots,y_{M}\\}$ sampled from $p_{U}(x)$ and
$p_{V}(y)$ are available. Thus one needs to calculate the statistical
distances between two samples, ${\cal D}_{U}$ and ${\cal D}_{V}$. To this end,
we employ the Marchenko-Pastur distances and the Wasserstein distances between
samples.

The linear XEB, Eq. (2) is zero if bit strings are sampled from classical
uniform random bits and a sample generated from a random circuit would be
different. Random matrix theory, here the Marchenko-Pastur distribution of
random matrices [47], could be helpful in capturing this difference. An
$M\times n$ binary array ${\cal D}=(x_{1},x_{2},\dots,,x_{M})^{t}$ is sliced
into into the collection of $k\times n$ rectangular random matrices $X$. Here
we set $k=2n$. If all entries of $X$ are independent and identically binary
random variables $\\{0,1\\}$, then one has the mean $\mu_{X}=1/2$ and the
variance $\sigma_{X}^{2}=1/4$. We calculate the distribution of eigenvalues of
square matrices $(1/k)X^{t}\cdot X$ as shown in Fig. 2. The eigenvalue
distribution of $X$ is composed of the two parts: the bulk distribution is the
Marchenko-Pastur distribution corresponding to random noise and the outlier
represents the signal [24, 25]. The random matrices $X$ can be written as
$Y=2X-J$ where $J$ of all entries 1, and $Y$ has the zero mean $\mu_{Y}=0$ and
the variance $\sigma^{2}_{Y}=1$. The square random matrix $X$ can be expressed
as

$\displaystyle\frac{1}{k}X^{t}\cdot X=\frac{1}{4k}\left(Y^{t}\cdot
Y+X^{t}\cdot J+X^{t}\cdot Z+J^{t}\cdot J\right)\,,$ (4)

Here, the eigenvalue distribution of the first term in Eq. (4) is given by the
Marchenko-Pastur distribution [47]

$\displaystyle\rho(\lambda)=\frac{1}{2\pi\sigma^{2}\gamma}\frac{\sqrt{(\lambda_{+}-\lambda)(\lambda-\lambda_{-})}}{\lambda}\,,$
(5)

where $\gamma=n/k$ is the rectangular ratio and
$\lambda_{\pm}=\sigma^{2}(1\pm\sqrt{\gamma})^{2}$ are the upper and lower
bounds. For $k=2n$, one has $\gamma=1/2$ and $\lambda_{+}=1+\sqrt{1/2}$. The
eigenvalues of the last term $(1/4k)J^{t}\cdot J$ are 0 and $n/4$. This is the
location of the outlier. We call the Marchenko-Pastur distance the distance of
the peak position of the outlier of a sample from $n/4$.

Fig. 2 (a) plots the Marchenko-Pastur distribution of eigenvalues of the
ensemble of $(1/2n){X^{t}\cdot X}$ for the Sycamore sample
(measurement-53-m20-s0-e0-pABCDCDAB.txt) and the Kalachev et al.’s sample
(samples-m20-f0-002.txt). The majority of the eigenvalues follow the
Marchenko-Pastur distribution given by Eq. (5). The inset of Fig. 2 (a) shows
the outliers of the distribution of eigenvalues around $n/4$, here
$53/4=13.25$ for $n=53$ qubits. The distance of the peak of the outliers from
$n/4$ is called the Marchenko-Pastur distance which is used to measure how a
sample is deviated from classical random bit strings. Fig 2 (b) plots the
Marchenko-Pastur distribution for Pan et al.’s spoofing sample (samples-
metropolis.txt) and Kalachev et al. spoofing file (spoofing-m20.txt). As
expected, the eigenvalue distributions of both spoofing samples are quite
different from those of the Sycamore samples or those of Kalachev et al.’s
samples. There is the big peak at eigenvalue zero, $\lambda=0$. Also, the
outliers are far away from $n/4$.

Fig. 2 (c) plot the Marchenko-Pastur distance of samples from $n/4$ as a
function of the number of cycles $m=12,14,16,18,$ and $20$. As shown in Fig. 2
(b), the Marchenko-Pastur distances of Kalachev et al.’s samples are very
close to that of the classical random bit strings, approximately 0.36. The
closeness of Kalachev et al.’s samples to the classical random bit strings are
consistent with the fact that Kalachev et al.’s samples pass the NIST random
number tests discussed above. We could not calculate the Marchenko-Pastur
distribution for Pan et al.’s sample because

Figure 3:  Wasserstein-distances of the Sycamore samples qubits with $n=53$
and the Kalachev et al.’s samples with $n=53$ from the classical random bit
sample.

In order to support the analysis above using the Marchenko-Pastur distances
between samples, we calculate the Wasserstein distances for samples [48]. In
contrast to the Kullback-Leibler divergence, the Wasserstein distance is a
metric on probability distributions. It is symmetric and satisfies the
triangular inequality. The $\alpha$-Wasserstein distance between two
probability distribution $p(x)$ and $q(y)$ is given by [48]

$\displaystyle
W(p,q)=\left(\inf_{\pi\in\Gamma(p,q)}\int_{\mathbb{R}\times\mathbb{R}}|x-y|^{\alpha}\,d\pi(x,y)\right)^{1/\alpha}\,,$
(6)

where $\Gamma(p,q)$ is the set of joint distributions $\pi(x,y)$ whose
marginals are $p(x)$ and $q(y)$. If $p(x)$ is the empirical distribution of a
dataset $D_{1}=\\{x_{1},x_{2},\dots,x_{M}\\}$ and $q(y)$ is the empirical
distribution of a dataset $D_{2}=\\{y_{1},y_{2},\dots,y_{M}\\}$, then
$\alpha$-Wasserstein distance is written as

$\displaystyle
W_{\alpha}(p,q)=\left(\sum_{i=1}^{M}|x_{(i)}-y_{(i)}|^{\alpha}\right)^{1/\alpha}\,,$
(7)

where $x_{(k)}$ denotes the order statistic of rank $k$, i.e., $k$-th smallest
value in the dataset ${\cal D}_{1}$. Here, we calculate the first Wasserstein
distance for the datasets of random circuit sampling using the python optimal
transport library [49]. To this end, each bit string $x$ is converted to a
value in range $[0,1]$ by dividing it by $1/2^{n}$. Also, we consider only 1
million bit strings. Fig. 3 plots the Wasserstein distances of samples from
the classical random bit strings as a function of the number of cycles $m$. It
is shown that the Wasserstein distances between Kalachev et al.’s samples and
the classical random bit strings are very close. This result is consistent
with the analysis with the Marchenko-Pastur distance. The Wasserstein
distances of Pan et al’s spoofing sample from the classical random sample, the
Sycamore sample, Kalachev et al.’s spoofing sample are about 0.0234, 0.01390,
and 0.3275, respectively.

##  III Conclusion

In summary, we have investigated the statistical properties of bit-strings
sampled from random circuits on the Sycamore quantum processor and obtained
using the tensor network simulation on classical processors. We considered
Kalachev et al.’s samples [17, 18, 28] and Pan et al.’s spoofing sample [19,
20, 29], and compared them with the Sycamore dataset [27] and the classical
random bit strings. We found that Kalachev et al.’s samples pass the NIST
random number tests while they have stripe patterns of the heat maps. The heat
maps of the two spoofing samples were shown to be quite different from those
of the Sycamore samples. Since some parts of bit-strings of the spoofing
samples are fixed to zeros or ones, the outcomes of the spoofing samples look
like that a random circuit applied only a few qubits. So one can easily
distinguish them from the Sycamore samples while the linear XEB values of the
spoofing samples are greater than zeros. Using the Marchenko-Pastur
distribution of eigenvalues of the samples and the Wasserstein distances
between the samples show that Kalachev et al.’s sample is close to classical
random bit strings and the spoofing samples are far away from the classical
samples or from the Sycamore samples.

Google’s 2019 experiment on random circuit sampling [9] was a landmark in
quantum computing and has been much attraction as well as intense debates
[50]. The tensor network simulation on quantum circuits has been proven to be
effective compared the Schrödinger-Feynman simulation and could catch up with
current noisy quantum processors [17, 18, 19, 20, 21]. All claims are based on
two metrics: computational time and the linear XEB. While the linear XEB may
serve as a good measure to verify quantum advantage, its limitation has been
pointed out [22, 17, 20]. Other measures could be used to verify quantum
computation for sampling [39]. For example, the random matrix analysis or the
Wasserstein distances [24, 25] would tell us other aspects of outcomes of
random circuit sampling. As more results from quantum processors and advanced
simulation on classical processors become available, a comparative study could
deepen our understanding of quantum advantage. Our results raise a question
about the claim that quantum advantage for random circuit sampling is faded by
classical computers [20, 21].

###### Acknowledgements.

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, National Quantum Information Science Research Centers. We
also acknowledge the National Science Foundation under Award number 1955907.

## References

  * Feynman [1982] R. P. Feynman, International Journal of Theoretical Physics 21, 467 (1982).
  * Shor [1997] P. W. Shor, SIAM Journal on Computing 26, 1484 (1997).
  * Harrow  _et al._ [2009] A. W. Harrow, A. Hassidim, and S. Lloyd, Phys. Rev. Lett. 103, 150502 (2009).
  * Boixo  _et al._ [2018] S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush, N. Ding, Z. Jiang, M. J. Bremner, J. M. Martinis, and H. Neven, Nature Physics 14, 595 (2018).
  * Neill  _et al._ [2018] C. Neill, P. Roushan, K. Kechedzhi, S. Boixo, S. V. Isakov, V. Smelyanskiy, A. Megrant, B. Chiaro, A. Dunsworth, K. Arya, R. Barends, B. Burkett, Y. Chen, Z. Chen, A. Fowler, B. Foxen, M. Giustina, R. Graff, E. Jeffrey, T. Huang, J. Kelly, P. Klimov, E. Lucero, J. Mutus, M. Neeley, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, H. Neven, and J. M. Martinis, Science 360, 195 (2018).
  * Bouland  _et al._ [2019] A. Bouland, B. Fefferman, C. Nirkhe, and U. Vazirani, Nature Physics 15, 159 (2019).
  * Scheel [2004] S. Scheel, Permanents in linear optical networks, https://arxiv.org/abs/quant-ph/0406127 (2004).
  * Aaronson and Arkhipov [2011] S. Aaronson and A. Arkhipov, in  _Proceedings of the Forty-Third Annual ACM Symposium on Theory of Computing_, STOC ’11 (Association for Computing Machinery, New York, NY, USA, 2011) p. 333–342.
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth, E. Farhi, B. Foxen, A. Fowler, C. Gidney, M. Giustina, R. Graff, K. Guerin, S. Habegger, M. P. Harrigan, M. J. Hartmann, A. Ho, M. Hoffmann, T. Huang, T. S. Humble, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly, P. V. Klimov, S. Knysh, A. Korotkov, F. Kostritsa, D. Landhuis, M. Lindmark, E. Lucero, D. Lyakh, S. Mandrà, J. R. McClean, M. McEwen, A. Megrant, X. Mi, K. Michielsen, M. Mohseni, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, E. Ostby, A. Petukhov, J. C. Platt, C. Quintana, E. G. Rieffel, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, K. J. Sung, M. D. Trevithick, A. Vainsencher, B. Villalonga, T. White, Z. J. Yao, P. Yeh, A. Zalcman, H. Neven, and J. M. Martinis, Nature 574, 505 (2019).
  * Wu  _et al._ [2021] Y. Wu, W.-S. Bao, S. Cao, F. Chen, M.-C. Chen, X. Chen, T.-H. Chung, H. Deng, Y. Du, D. Fan, M. Gong, C. Guo, C. Guo, S. Guo, L. Han, L. Hong, H.-L. Huang, Y.-H. Huo, L. Li, N. Li, S. Li, Y. Li, F. Liang, C. Lin, J. Lin, H. Qian, D. Qiao, H. Rong, H. Su, L. Sun, L. Wang, S. Wang, D. Wu, Y. Xu, K. Yan, W. Yang, Y. Yang, Y. Ye, J. Yin, C. Ying, J. Yu, C. Zha, C. Zhang, H. Zhang, K. Zhang, Y. Zhang, H. Zhao, Y. Zhao, L. Zhou, Q. Zhu, C.-Y. Lu, C.-Z. Peng, X. Zhu, and J.-W. Pan, Phys. Rev. Lett. 127, 180501 (2021).
  * Zhu  _et al._ [2022] Q. Zhu, S. Cao, F. Chen, M.-C. Chen, X. Chen, T.-H. Chung, H. Deng, Y. Du, D. Fan, M. Gong, C. Guo, C. Guo, S. Guo, L. Han, L. Hong, H.-L. Huang, Y.-H. Huo, L. Li, N. Li, S. Li, Y. Li, F. Liang, C. Lin, J. Lin, H. Qian, D. Qiao, H. Rong, H. Su, L. Sun, L. Wang, S. Wang, D. Wu, Y. Wu, Y. Xu, K. Yan, W. Yang, Y. Yang, Y. Ye, J. Yin, C. Ying, J. Yu, C. Zha, C. Zhang, H. Zhang, K. Zhang, Y. Zhang, H. Zhao, Y. Zhao, L. Zhou, C.-Y. Lu, C.-Z. Peng, X. Zhu, and J.-W. Pan, Science Bulletin 67, 240 (2022).
  * Zhong  _et al._ [2020] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, P. Hu, X.-Y. Yang, W.-J. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N.-L. Liu, C.-Y. Lu, and J.-W. Pan, Science 370, 1460 (2020).
  * Zhong  _et al._ [2021] H.-S. Zhong, Y.-H. Deng, J. Qin, H. Wang, M.-C. Chen, L.-C. Peng, Y.-H. Luo, D. Wu, S.-Q. Gong, H. Su, Y. Hu, P. Hu, X.-Y. Yang, W.-J. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N.-L. Liu, J. J. Renema, C.-Y. Lu, and J.-W. Pan, Phys. Rev. Lett. 127, 180502 (2021).
  * Huang  _et al._ [2020] C. Huang, F. Zhang, M. Newman, J. Cai, X. Gao, Z. Tian, J. Wu, H. Xu, H. Yu, B. Yuan, M. Szegedy, Y. Shi, and J. Chen, Classical simulation of quantum supremacy circuits, https://arxiv.org/abs/2005.06787 (2020).
  * Gray and Kourtis [2021] J. Gray and S. Kourtis, Quantum 5, 410 (2021).
  * Guo  _et al._ [2021] C. Guo, Y. Zhao, and H.-L. Huang, Phys. Rev. Lett. 126, 070502 (2021).
  * Kalachev  _et al._ [2021a] G. Kalachev, P. Panteleev, P. Zhou, and M.-H. Yung, Classical sampling of random quantum circuits with bounded fidelity, https://arxiv.org/abs/2112.15083 (2021a).
  * Kalachev  _et al._ [2021b] G. Kalachev, P. Panteleev, and M.-H. Yung, Multi-tensor contraction for xeb verification of quantum circuits, https://arxiv.org/abs/2108.05665 (2021b).
  * Pan and Zhang [2022] F. Pan and P. Zhang, Phys. Rev. Lett. 128, 030501 (2022).
  * Pan  _et al._ [2022a] F. Pan, K. Chen, and P. Zhang, Phys. Rev. Lett. 129, 090502 (2022a).
  * Cho [2022] A. Cho, Science 377, 563 (2022).
  * Gao  _et al._ [2021] X. Gao, M. Kalinowski, C.-N. Chou, M. D. Lukin, B. Barak, and S. Choi, Limitations of linear cross-entropy as a measure for quantum advantage, https://arxiv.org/abs/2112.01657 (2021), arXiv:2112.01657 [quant-ph] .
  * Barak  _et al._ [2021] B. Barak, C.-N. Chou, and X. Gao, in  _12th Innovations in Theoretical Computer Science Conference (ITCS 2021)_, Leibniz International Proceedings in Informatics (LIPIcs), Vol. 185, edited by J. R. Lee (Schloss Dagstuhl–Leibniz-Zentrum für Informatik, Dagstuhl, Germany, 2021) pp. 30:1–30:20.
  * Oh and Kais [2022a] S. Oh and S. Kais, J. Phys. Chem. Lett. 13, 7469 (2022a).
  * Oh and Kais [2022b] S. Oh and S. Kais, Phys. Rev. A 106, 032433 (2022b).
  * Eisert  _et al._ [2020] J. Eisert, D. Hangleiter, N. Walk, I. Roth, D. Markham, R. Parekh, U. Chabaud, and E. Kashefi, Nature Reviews Physics 2, 382 (2020).
  * Martinis et al. [2022] J. M. Martinis et al., Quantum supremacy using a programmable superconducting processor, Dryad, Dataset, https://doi.org/10.5061/dryad.k6t1rj8 (2022).
  * Kalachev  _et al._ [2022] G. Kalachev, P. Panteleev, and M.-H. Yung, Huawei-hiq/supremacy, https://gitee.com/Huawei-HiQ/supremacy/tree/master (2022), accessed:2022-11-07.
  * Pan  _et al._ [2022b] F. Pan, K. Chen, and P. Zhang, solve_sycamore, https://github.com/Fanerst/solve_sycamore (2022b), accessed:2022-11-07.
  * Dyson [1962] F. J. Dyson, Journal of Mathematical Physics 3, 1199 (1962).
  * Hurwitz [1897] A. Hurwitz, Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse 1897, 71 (1897).
  * Zyczkowski and Kuś [1994] K. Zyczkowski and M. Kuś, Journal of Physics A: Mathematical and General 27, 4235 (1994).
  * Diaconis and Forrester [2017] P. Diaconis and P. J. Forrester, Random Matrices Theory and Application 06, 1730001 (2017), http://arxiv.org/pdf/1512.09229.
  * Mezzadri [2007] F. Mezzadri, Notices of the American Mathematical Society 54, 592 (2007).
  * Emerson  _et al._ [2003] J. Emerson, Y. S. Weinstein, M. Saraceno, S. Lloyd, and D. G. Cory, Science 302, 2098 (2003).
  * Brody  _et al._ [1981] T. A. Brody, J. Flores, J. B. French, P. A. Mello, A. Pandey, and S. S. M. Wong, Rev. Mod. Phys. 53, 385 (1981).
  * Kuś  _et al._ [1988] M. Kuś, J. Mostowski, and F. Haake, Journal of Physics A: Mathematical and General 21, L1073 (1988).
  * Haake  _et al._ [2010] F. Haake, S. Gnutzmann, and M. Kuś,  _Quantum Signatures of Chaos_ (Springer-Verlag Berlin Heidelberg, 2010).
  * Gilchrist  _et al._ [2005] A. Gilchrist, N. K. Langford, and M. A. Nielsen, Phys. Rev. A 71, 062310 (2005).
  * Bassham  _et al._ [2010] L. Bassham, A. Rukhin, J. Soto, J. Nechvatal, M. Smid, S. Leigh, M. Levenson, M. Vangel, N. Heckert, and D. Banks, A statistical test suite for random and pseudorandom number generators for cryptographic applications (2010).
  * Shi  _et al._ [2022] J. Shi, T. Zhao, Y. Wang, C. Yu, Y. Lu, R. Shi, S. Zhang, and J. Wu, An unbiased quantum random number generator based on boson sampling, https://arxiv.org/abs/2206.02292 (2022).
  * Herrero-Collantes and Garcia-Escartin [2017] M. Herrero-Collantes and J. C. Garcia-Escartin, Rev. Mod. Phys. 89, 015004 (2017).
  * Ang [2019] S. K. Ang, Randomness_testsuite, https://github.com/stevenang/randomness_testsuite#readme (2019).
  * [44] Supplemental material contains the results of the nist random number tests for six samples. also, it includes the 4 animation files to illustrate outcomes of bits for the 2-dimensional array of 53 qubits of the sycamore processor.
  * Rinott  _et al._ [2020] Y. Rinott, T. Shoham, and G. Kalai, Statistical aspects of the quantum supremacy demonstration, https://arxiv.org/abs/2008.05177 (2020).
  * Nielsen and Chuang [2011] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ , 10th ed. (Cambridge University Press, USA, 2011) pp. 194–195.
  * Marchenko and Pastur [1967] V. A. Marchenko and L. A. Pastur, Mat. Sb. (N.S.) 72, 507 (1967).
  * Villani [2008] C. Villani,  _Optimal Transport: Old and New_, Grundlehren der mathematischen Wissenschaften (Springer Berlin Heidelberg, 2008).
  * Flamary  _et al._ [2021] R. Flamary, N. Courty, A. Gramfort, M. Z. Alaya, A. Boisbunon, S. Chambon, L. Chapel, A. Corenflos, K. Fatras, N. Fournier, L. Gautheron, N. T. Gayraud, H. Janati, A. Rakotomamonjy, I. Redko, A. Rolet, A. Schutz, V. Seguy, D. J. Sutherland, R. Tavenard, A. Tong, and T. Vayer, Journal of Machine Learning Research 22, 1 (2021).
  * Kalai  _et al._ [2022] G. Kalai, Y. Rinott, and T. Shoham, Google’s 2019 ”quantum supremacy” claims: Data, documentation, and discussion, https://arxiv.org/abs/2210.12753 (2022).