# Improved TokenPose with Sparsity

Anning Li
University of Electronic Science and Technology of China
<EMAIL_ADDRESS>

###### Abstract

Over the past few years, the vision transformer and its various forms have
gained significance in human pose estimation. By treating image patches as
tokens, transformers can capture global relationships wisely, estimate the
keypoint tokens by leveraging the visual tokens, and recognize the posture of
the human body. Nevertheless, global attention is computationally demanding,
which poses a challenge for scaling up transformer-based methods to high-
resolution features. In this paper, we introduce sparsity in both keypoint
token attention and visual token attention to improve human pose estimation.
Experimental results on the MPII dataset demonstrate that our model has a
higher level of accuracy and proved the feasibility of the method, achieving
new state-of-the-art results. The idea can also provide references for other
transformer-based models.

_K_ eywords Vision Transformer $\cdot$ Human Pose Estimation $\cdot$ Machine
Learning $\cdot$ Self-attention mechanism

##  1 Introduction

Human pose estimation Is the process of recovering anatomical keypoints of the
human body given an image or a video, which deeply relies on both visual cues
and keypoints constraint relationships. The most common approach for human
pose estimation is to use convolutional neural networks (CNNs) to extract
features from the input image [1, 2, 3, 4, 5, 6, 7]. However, they have some
limitations and challenges. Human body joints are highly correlated, and
constrained by strong kinetic and physical restrictions. [8]. Therefore, the
limited reception fields of convolutional neural networks (CNNs) can result in
poor capture of the long-range constraints that exist between joints[9].

The transformer-based Vision Transformer [10] and its variants have showcased
impressive performance on various visual tasks[11, 12]. In comparison with
CNN, the self-attention mechanism of transformers can capture long-range
dependencies[13] and learn the relationship among keypoints better. In single-
view 2D human pose estimation, TransPose [14] and TokenPose [9] have achieved
new state-of-the-art performances. Typically, the feature maps are flattened
into unidimensional token sequences and then fed into the transformer.
TokenPose divides the keypoints into two categories, visual tokens and
keypoint tokens. Visual tokens generally refer to the visual features or
representations used to encode the appearance and context information in an
image or video frame. Keypoint tokens, on the other hand, refer to the
specific body keypoints that are of interest in the pose estimation task,
represented as coordinates $(x,y)$ or $(x,y,z)$ in the image or 3D space.

However, The computation complexity is expensive due to dimensional
transformation, with the amount being quadratic to the length of input
sequences. Therefore, it is challenging to extend these methods to high-
resolution feature maps. What’s more, the attention maps are generally very
sparse and only attend to a small local region. Most of the redundant areas
can be mowed to speed up the training process. This also prompts our interest
to strengthen the connections between keypoint tokens, since we can generalize
mathematical matrices through certain joint contacts in human anatomy.

Figure 1: Attention map is usually very sparse and only attends to a small
local region.

In this paper, we propose a more efficient, yet accurate Transformer model
based on the TokenPose Transformer[9], named Improve TokenPose with Sparsity.
Specifically, we focus on two types of sparse tokens, human-joint keypoint
tokens, and attention-level visual tokens. On keypoint tokens, we refer to
human physiological structures to strengthen the connection of relevant joints
with a sparse joint mask. On visual tokens, we introduce sparsity at the
attention level by pruning attention connections. We reveal that pruning these
sparse tokens does not reduce the accuracy but can accelerate the entire
network.

Our main contributions to the Improved TokenPose with Sparsity are summarized
as follows:

  1. 1.

We propose an attention-based pruning strategy to dynamically upgrade a sparse
attention mask. By focusing on the more significant visual tokens, we can
simultaneously achieve both token and attention connections during training
while making use of their sparsities. Our method does not require a complex
token selector module or additional training loss or hyper-parameters.

  2. 2.

We propose a self-defined sparse joint mask based on human structure. With the
use of graph degree, we can pre-define a 0-1 mask to emphasize the connections
between adjacent joints and symmetrical joints.

  3. 3.

We conduct extensive experiments on MPII dataset, and prove the feasibility of
the method by experimental analysis.

##  2 Related Work

###  2.1 Vision Transformer

Vision transformer [15] was originally proposed for natural language
processing, which has been shown to be effective in capturing long-range
dependencies in sequences of tokens. The key idea of ViT is to apply the
Transformer architecture directly to image patches. Recently, they have shown
remarkable success in various applications such as classification [10, 11],
object detection[12, 16, 17], and semantic segmentation[18, 19, 20, 21]. Since
its introduction, ViT has achieved state-of-the-art performance on several
benchmark datasets for image recognition, including ImageNet[22]. ViT has also
inspired further research in the area of vision transformers, including ViT
variants with improved performance, such as DeiT [23], and methods for
efficient training of vision transformers, such as TNT[24]. However, their
extensive computation and memory usage during both training and inference
hinder their generalization. Existing compression algorithms typically focus
on efficient inference and rely on pre-trained dense models, neglecting time-
consuming training. Researchers have proposed several methods to improve the
efficiency of ViT. For example, Network pruning[25, 26, 27, 28], the vanilla
ViT [10], knowledge distillation [29, 11, 30], and quantization[31, 32] have
been applied to ViT, and other algorithms introduce CNN properties to
alleviate the burden of computing global attention[33, 34].

Additionally, some models slim the input tokens by aggregating neighboring
tokens[35, 36, 37, 38, 39, 40, 41]. Their core idea is mining tokens using
learnable attention weights or reducing and reorganizing image tokens based on
the classification token. For example the Token-to-tokens[35], the
TokenLearner,[37], the DynamicViT [38] and the EViT [40]. These models,
however, have only been designed for classification tasks, where the ultimate
prediction relies solely on the classification token, whereas in human pose
estimation, it is also important to take spatial clues and articular
association into consideration.

###  2.2 Human Pose Estimation

Human pose estimation is an important research area in computer vision that
involves estimating the joint locations and orientations of a person’s body
from images or videos. Previously, many CNN models were proposed in 2D human
pose estimation[2, 4, 5]. Despite its recognized effectiveness in capturing
local information, convolutional neural networks still struggle with modeling
long-range relationships. This is particularly challenging when it comes to
human pose estimation. Therefore, many works have been done to explore
efficient architecture design for human pose estimation. Such as HRNet[5], is
a multi-stage, high-resolution network designed specifically for human pose
estimation. It uses a high-resolution feature extraction stage, followed by a
set of parallel branches that combine information across different scales to
achieve state-of-the-art results. Lite-HRNet [42] is a lightweight version of
the HRNet with limited computational resources. It achieves this by reducing
the number of channels and layers in the original HRNet while retaining its
multi-resolution feature extraction and integration capabilities. Other works
are also designed for real-time pose estimation [43, 44]. However, their focus
is primarily on CNN-based networks, and none of their studies explore
transformer-based networks.

Recently, many works have used transformer networks for human pose estimation.
For example, TransPose utilizes transformers to illustrate the interdependence
of keypoint predictions[14], TokenPose incorporates extra keypoint tokens to
acquire knowledge of constraint relationships and appearance cues[9] and
employs a transformer module to capture the long-range dependencies between
body parts in human pose estimation for video object segmentation[45]. These
works highlight the growing use of transformer-based models in human pose
estimation. By leveraging self-attention mechanisms and capturing global
dependencies, transformer-based approaches have demonstrated promising results
and opened up new avenues for advancing the accuracy and robustness of pose
estimation systems.

Both CNNs and Transformer models are effective in human pose estimation tasks.
CNNs are well-suited for extracting local features and handling spatial
variability, while Transformers are better at capturing complex global
relationships between body parts. Both emphasize focus on global information,
efficiency, and robustness to noise.

##  3 Methodology

###  3.1 Overview

Figure 2: A schematic diagram of the Improved TokenPose.

Fig.2 is an overview of our Improved TokenPose with Sparsity. The input image
$\mathbf{I}$ first goes through a CNN backbone and outputs the feature map
$\mathbf{F}\in\mathbb{R}^{H\times W\times C}$. Then the feature map is split
into patches $\mathbf{F}_{p}\in\mathbb{R}^{N_{P}\times(C\cdot P_{h}\cdot
P_{w})}$ and flattened to 1D vectors $\mathbf{X}_{p}\in\mathbb{R}^{N_{v}\times
D}$. $(P_{h},P_{w})$ is the resolution of each image patch. $D$ is the
dimension of hidden embeddings. $N_{p}=\frac{H}{P_{h}}\cdot\frac{W}{P_{w}}$ is
the total number of patches. Then, by using linear projection, the flattened
vectors produce visual tokens $\mathbf{X}_{v}=\mathbf{X}_{p}+\mathbf{E}$.
$\mathbf{E}$ is the 2D positional encoding
$\mathbf{E}\in\mathbb{R}^{N_{p}\times D}$. $J$ specific types of tokens are
randomly initialized to represent individual keypoint tokens
$\mathbf{X}_{k}\in\mathbb{R}^{J\times D}$. The visual tokens and keypoint
tokens are then fed into the Transformer encoder. The input to the transformer
is
$\mathbf{X}^{0}=[\mathbf{X}_{k},\mathbf{X}_{v}]\in\mathbb{R}^{(N_{p}+J)\times
D}$.

The transformer has $L$ encoder layers in total. The encoder layer consists of
multi-headed self-attention. The multi-head attention in each layer captures
both appearance and anatomical constraint cues. At chosen
$L_{3}^{th},L_{6}^{th},L_{9}^{th}$ layers, we employ the masked multi-head
self-attention mechanism to make use of the sparsity in visual token attention
and better measure the complexity of images(see 3.2).

After the $L$ encoder layers, the output keypoint tokens are then put into the
$L_{1}(L_{1}\textless L)$ layers graph transformer with a self-defined sparse
joint mask $\mathbf{M_{J}}\in\mathbb{R}^{J\times J}$, which is initialized to
represent the connection of human joints. The graph transformer has the same
structure as the transformer, only this time we consider the joint mask as the
initialized attention mask. Through the graph transformer layers, we once
again emphasize the relevance of the human body structure. In the end, an MLP
head uses the keypoint tokens $\mathbf{X}_{k}^{L}\in\mathbb{R}^{J\times D}$
from the last Graph Transformer layer to predict the keypoints heatmaps
$\mathbf{H}\in\mathbb{R}^{J\times(H_{h}\cdot W_{h})}$.

Figure 3: Image map

Fig.3 is a schematic diagram of our proceed tokens. Based on the TokenPose
model[9], in order to improve the accuracy of its application in human pose
estimation, we propose improvements in two places mainly focusing on the upper
left keypoint token joint mask and lower right visual token attention mask.
Since the image size passed to the transformer encoder is
${(HW+J)\times(HW+J)}$, we can abstract a diagram like Fig.3 to help us
understand the tokens we are processing. By pruning the two areas ${HW\times
HW}$ and ${J\times J}$, and using the sparsity of the chosen tokens, we can
remove the relatively less important data and still get an accurate ultimate
result.

###  3.2 Masked Multi-head Self-attention(MMSA)

The Masked Multi-head Self-attention aims to match a query and a set of key-
value pairs to an output[15]. Given the inputs $\mathbf{X}$ and $\mathbf{M}$,
three linear projections are applied to transfer $\mathbf{X}$ into three
matrices of equal size, namely the query $\mathbf{Q}$, the key $\mathbf{K}$,
and the value $\mathbf{V}$. The Masked Multi-head Self-attention (MMSA)
operation is calculated by:

$\text{MMSA}(\mathbf{X})=\text{Softmax}(\frac{\mathbf{Q}\mathbf{K}^{T}\times\mathbf{M}}{\sqrt{D}})\mathbf{V},$
(1)

According to Fig.2, the tokens go through two different Transformer encoders
in turn. For the first Transformer encoder, in the process of MMSA, the inputs
are visual tokens $\mathbf{X}\in\mathbb{R}^{H\times W\times D}$ and attention
mask $\mathbf{M_{A}}\in\mathbb{R}^{J\times J}$. The attention mask
$\mathbf{M_{A}}\in\mathbb{R}^{HW\times HW}$ is initially set to an all-one
matrix, and then extracted from certain layers through attention pruning. In
this case, the attention mask generated from the $L_{3}^{th}$ layer will be
used in $L_{4}^{th}$ and $L_{5}^{th}$ layers, and will be dynamically pruned
and updated in the $L_{6}^{th}$ layer. Later in the second Transformer
encoder, the inputs for MMSA are keypoint tokens
$\mathbf{X}_{k}\in\mathbb{R}^{J\times D}$ and the self-defined joint mask
$\mathbf{M_{J}}\in\mathbb{R}^{J\times J}$. The calculation is the same as
Equation1.

###  3.3 Visual Token Attention Prune Module

We introduce sparsity at the attention level on visual tokens by pruning
attention connections. Specifically, for a given query token, we only
calculate its attention connections with a small subset of selected tokens.
The attention mask aims at determining the attention weights between different
input tokens. Since attention maps are typically sparse, it is unnecessary to
pursue intensive attention, and pruning it can selectively remove connections
that are deemed less important.

Prior versions of ViT have relied on predetermined sparse patterns, such as
neighborhood windows[46] or dilated sliding windows[47], which can be limited
for objects with varying shapes and sizes. Alternatively, some input-dependent
methods have used deformable attention to find corresponding tokens[48], but
this requires additional modules to learn the selected tokens. The Tri-level
ViT[49] proposed a relatively feasible scheme, but no concrete experiments
have been carried out on it. This paper intends to implement the method of
human posture detection and explore its influence on model accuracy and
efficiency. We employ the MMSA mechanism in the transformer to better capture
the complexity of images, as shown in Fig.2.

Figure 4: The extraction of sparse attention mask

As shown in Fig.4, given the image matrix $\mathbf{X}$ and $\mathbf{M}$, for
each row we dynamically select the top K values and form a new 0-1 attention
mask. After the calculation and normalization of the self-attention matrix, we
can simplify the attention mask. By dynamically choosing necessary tokens and
pruning the less relevant smaller numbers, we can remove the redundant data to
improve efficiency.

We initialize all elements in the mask to 1 and update the mask progressively.
The prediction modules take the current $\mathbf{M}\in\mathbb{R}^{N\times N}$
(When considering joint mask as input the mask size is
$\mathbf{M}\in\mathbb{R}^{N\times C}$) and the tokens
$\mathbf{X}\in\mathbb{R}^{N\times D}$ as input. $\mathbf{N}$ stands for
$\mathbb{H\times W}$. We compute the weighted sum over all values $v$ in the
sequence tokens, together with their respective query $q^{i}$ and key $k^{j}$
representations.

$\mathbf{[q,k,v]=zU_{qkv},\quad U_{qkv}\in\mathbb{R}^{D\times 3D_{h}}}$ (2)

The attention weights $A_{ij}$ are based on not only the pairwise similarity
between two elements of the sequence and their respective $q^{i}$ and $k^{j}$,
but also the Attention mask to progressively update the more important
attention weights. We will multiply $A$ with our input map $M$ to get a
simpler matrix $A^{\prime}$

$\mathbf{A^{\prime}=\mathbf{A}\times\mathbf{M},\quad A\in\mathbb{R}^{N\times
N},\;M\in\mathbb{R}^{N\times N}}$ (3)

Notice that we are only doing softmax when the value in the matrix $M_{ij}$ is
1. If $M_{ij}$ is 0 we will skip the calculation process, making the operation
easier.

$\text{softmax}\mathbf{(A^{\prime}_{ij})=\frac{exp(A^{\prime}_{ij})}{{\textstyle\sum_{j=0}^{J}exp(A^{\prime}_{ij})(\text{only
\, when}\;M_{ij}=1)}}}$ (4)

The superiority of the attention prune module is manifested in two aspects.
Firstly the attention mask is generally very sparse, there is no need to
consider all of them. Secondly, it is approved that attention is not the
bigger the better. Some researchers tended to throw away all the attention,
but this may result in losing information about the image. So in our method,
we are trying to make it shorter while still keeping the valid information, in
order to achieve substantial performance improvement. According to the
analysis of each operation in a Transformer block [50], the pruning operations
are performed on weight tensors of each attention head in the Multi-Self
Attention (MSA) module. Given an input sequence, $N\times D$, in which $N$ is
the input sequence length or the token number, and $D$ is the embedding
dimension of each token. In MSA, the input size was originally $N\times N$,
but by doing visual token attention pruning, we can reduce the input size to
$N\times D$, further simplifying the calculation.

###  3.4 Joint Attention Mask

Figure 5: The extraction of human joint mask

As shown in Fig.5. No matter what the posture of the human body is, the
connection between the joints is always there and has a certain pattern.
Therefore we can be able to summarize and pre-define a keypoint joint mask
$\mathbf{M_{J}}\in\mathbb{R}^{J\times J}$ based on the human graph to
emphasize the connection between limbs and strengthen the joint relationship
during training. By using the sparsity of the junction, for each human joint,
we only keep its adjacent nodes and symmetric nodes, set them to 1, and set
the other nodes to 0, to form a sparse joint mask. For instance, for joint
number $5$, the left shoulder, we only keep its neighboring left elbow joint
and its symmetric joint, the right shoulder.

A corresponding operation is performed on each joint of the human body
structure to form a ${J\times J}$ mask. It should be used as a constant value
and add to the whole process of training. The joint mask will later fed into
the Graph Transformer layer(as in Fig.2) to emphasize the invariability of
human joint connection.

###  3.5 Implementation Detail

On the Improved TokenPose with Sparsity, the number of encoder layers for
Transformer $L$ is set to $12$, the embedding size $D$ is set to $192$, and
the number of heads $H$ is set to $8$. And in Graph Transformer the number of
encoder layers $L^{\prime}$ is set to $8$. They take the shallow stem-net and
the HRNet-W32 as the CNN backbone.

##  4 Experiments and Validation

###  4.1 Settings

We evaluate our Improved TokenPose with Sparsity on monocular 2D human pose
estimation benchmarks, datasets MPII[51] contains about $25K$ images and $40K$
human instances with $16$ keypoints. The evaluation is based on the head-
normalized probability of correct keypoint (PCKh) score [51]. A keypoint is
correct if it falls within a predefined threshold to the ground-truth
location. We report the mean PCKh score by convention.

###  4.2 Ablation Studies

Method | Head | Sho | Elb | Wri | Hip | Kne | Ank | Mean |<EMAIL_ADDRESS>
---|---|---|---|---|---|---|---|---|---
TokenPose-S [9] | 96.3 | 93.9 | 86.1 | 79.2 | 86.2 | 79.6 | 74.5 | 85.8 | 30.9
TokenPose with joint mask | 96.3 | 94.2 | 87.1 | 80.0 | 85.8 | 80.8 | 75.8 | 86.3 | 31.17
Improved TokenPose( AKR = 0.2 ) | 95.8 | 93.3 | 85.4 | 78.3 | 85.4 | 79.4 | 74.2 | 85.2 | 29.98
Improved TokenPose( AKR = 0.5 ) | 96.1 | 94.3 | 87.3 | 80.3 | 86.7 | 81.1 | 76.2 | 86.7 | 31.23
Improved TokenPose( AKR = 0.6 ) | 96.1 | 94.0 | 87.0 | 80.6 | 86.9 | 81.4 | 75.9 | 86.8 | 31.28
Improved TokenPose( AKR = 0.7 ) | 96.0 | 94.0 | 86.7 | 79.8 | 85.8 | 80.0 | 73.9 | 85.9 | 31.24
Improved TokenPose( AKR = 0.8 ) | 96.3 | 94.5 | 87.0 | 80.1 | 86.8 | 81.2 | 75.5 | 86.6 | 31.11

Table 1:  Ablation studies on different Attention-Keep Ratio(AKR)

We first build our Improved TokenPose with Sparsity based on TokenPose-S for
comparison. We follow the same training recipes as TokenPose[9]. Table 1 shows
the progressive results of the human pose evaluation. Adding the joint mask
can bring a slight improvement. Making use of the sparsity in both visual
token attention and keypoint token attention can have an improvement of at
least one percentage point. As shown in the results, the accuracy is the best
when the attention-keep ratio is 0.6, which means focusing only 60 percent of
the visual tokens and throwing away the rest of them.

The results in Table 1 proved that the pruned visual token connections might
correspond to less critical relationships, and help with highlighting the more
influential relationships in the human pose estimation process. However, it’s
important to note that the pruning process may impact the overall accuracy of
the pose estimation model. There is definitely a trade-off. Removing certain
connections may discard valuable information, making the model lose some fine-
grained details of human poses. It still requires careful analysis and
experimentation to strike a balance between pruning for efficiency and
preserving accuracy.

###  4.3 Results

Method | Head | Sho | Elb | Wri | Hip | Kne | Ank | Mean
---|---|---|---|---|---|---|---|---
SimpleBaseline-Res50 [4] | 96.4 | 95.3 | 89.0 | 83.2 | 88.4 | 84.0 | 79.6 | 88.5
HRNet-W32. [5] | 96.9 | 96.0 | 90.6 | 85.8 | 88.7 | 86.6 | 82.6 | 90.1
TokenPose-S [9] | 96.3 | 93.9 | 86.1 | 79.2 | 86.2 | 79.6 | 74.5 | 85.8
Improved TokenPose-S | 96.1 | 94.1 | 87.0 | 80.6 | 87.0 | 81.4 | 75.9 | 86.8(+1.0)

Table 2: Results on the MPII validation set. The input size is $256\times
256$.

The experimental results on the MPII dataset are shown in Table 2. As we can
see our methods have some improvements compared with the original TokenPose
baseline. The idea of making use of attention-level visual tokens and keypoint
tokens sparsity is proven to be efficient and effective, and this idea can be
applied to other fields as well.

Figure 6: Sample heatmaps of our approach

###  4.4 Visualizations

We visualize the sample heatmaps of our approach in Fig.6 and Fig.8. We
visualize the procedure of attention pruning. Fig.7 shows how the attention
maps are pruned are dynamically select the important information during the
procedure. After the pruning, we will drop a few tokens that don’t make much
difference to the final result.

Figure 7: The process of attention pruning  Figure 8: Heatmaps with recognized
joints on MPII dataset

##  5 Conclusion

In this paper, we propose the Improved TokenPose with Sparsity for 2D human
pose estimation. Our superiority is reflected in two parts: the combination of
human body structure, and the use of attention pruning to generate a simpler
attention map. We reached improvements for applying the sparsity in keypoint
and visual tokens. Experiments on the MPII dataset show that the Improved
TokenPose with Sparsity achieves better accuracy compared with previous
TokenPose networks.

## References

  * [1] Alexander Toshev and Christian Szegedy.  Deeppose: Human pose estimation via deep neural networks.  In CVPR, 2014.
  * [2] Shih-En Wei, Varun Ramakrishna, Takeo Kanade, and Yaser Sheikh.  Convolutional pose machines.  In CVPR, 2016.
  * [3] Alejandro Newell, Kaiyu Yang, and Jia Deng.  Stacked hourglass networks for human pose estimation.  In ECCV, 2016.
  * [4] Bin Xiao, Haiping Wu, and Yichen Wei.  Simple baselines for human pose estimation and tracking.  In ECCV, 2018.
  * [5] Ke Sun, Bin Xiao, Dong Liu, and Jingdong Wang.  Deep high-resolution representation learning for human pose estimation.  In CVPR, 2019.
  * [6] Z. Wang, J. Yang, and C. Fowlkes.  The best of both worlds: Combining model-based and nonparametric approaches for 3d human body estimation.  In CVPR ABAW workshop, 2022.
  * [7] Zhe Wang, Daeyun Shin, and Charless C Fowlkes.  Predicting camera viewpoint improves cross-dataset generalization for 3d human pose estimation.  In ECCV 3DPW workshop, 2020.
  * [8] Jonathan J Tompson, Arjun Jain, Yann LeCun, and Christoph Bregler.  Joint training of a convolutional network and a graphical model for human pose estimation.  NIPS, 2014.
  * [9] Yanjie Li, Shoukui Zhang, Zhicheng Wang, Sen Yang, Wankou Yang, Shu-Tao Xia, and Erjin Zhou.  Tokenpose: Learning keypoint tokens for human pose estimation.  In ICCV, 2021.
  * [10] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  ICLR, 2021.
  * [11] Hugo Touvron, Matthieu Cord, Matthijs Douze, Francisco Massa, Alexandre Sablayrolles, and Hervé Jégou.  Training data-efficient image transformers & distillation through attention.  ICML, 2021.
  * [12] Nicolas Carion, Francisco Massa, Gabriel Synnaeve, Nicolas Usunier, Alexander Kirillov, and Sergey Zagoruyko.  End-to-end object detection with transformers.  In ECCV, 2020.
  * [13] Irwan Bello, Barret Zoph, Ashish Vaswani, Jonathon Shlens, and Quoc V. Le.  Attention augmented convolutional networks.  In Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV), October 2019.
  * [14] Sen Yang, Zhibin Quan, Mu Nie, and Wankou Yang.  Transpose: Keypoint localization via transformer.  In ICCV, 2021.
  * [15] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  NIPS, 2017.
  * [16] Xizhou Zhu, Weijie Su, Lewei Lu, Bin Li, Xiaogang Wang, and Jifeng Dai.  Deformable detr: Deformable transformers for end-to-end object detection.  ICLR, 2021.
  * [17] Yuxin Fang, Bencheng Liao, Xinggang Wang, Jiemin Fang, Jiyang Qi, Rui Wu, Jianwei Niu, and Wenyu Liu.  You only look at one sequence: Rethinking transformer in vision through object detection.  NeurIPS, 2021.
  * [18] Sixiao Zheng, Jiachen Lu, Hengshuang Zhao, Xiatian Zhu, Zekun Luo, Yabiao Wang, Yanwei Fu, Jianfeng Feng, Tao Xiang, Philip HS Torr, et al.  Rethinking semantic segmentation from a sequence-to-sequence perspective with transformers.  CVPR, 2021.
  * [19] Yuqing Wang, Zhaoliang Xu, Xinlong Wang, Chunhua Shen, Baoshan Cheng, Hao Shen, and Huaxia Xia.  End-to-end video instance segmentation with transformers.  CVPR, 2021.
  * [20] Xiangyi Yan, Hao Tang, Shanlin Sun, Haoyu Ma, Deying Kong, and Xiaohui Xie.  After-unet: Axial fusion transformer unet for medical image segmentation.  In WACV, 2022.
  * [21] Chenyu You, Ruihan Zhao, Fenglin Liu, Sandeep Chinchali, Ufuk Topcu, Lawrence Staib, and James S Duncan.  Class-aware generative adversarial transformers for medical image segmentation.  arXiv preprint arXiv:2201.10737, 2022.
  * [22] Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael S. Bernstein, Alexander C. Berg, and Li Fei-Fei.  Imagenet large scale visual recognition challenge.  International Journal of Computer Vision, 115:211–252, 2014.
  * [23] Hugo Touvron, Matthieu Cord, Matthijs Douze, Francisco Massa, Alexandre Sablayrolles, and Herve Jegou.  Training data-efficient image transformers & distillation through attention.  In Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, pages 10347–10357, 2021.
  * [24] Yuchen Dang, Zheyuan Hu, Miles Cranmer, Michael Eickenberg, and Shirley Ho.  Tnt: Vision transformer for turbulence simulations, 2022.
  * [25] Song Han, Huizi Mao, and William J Dally.  Deep compression: Compressing deep neural networks with pruning, trained quantization and huffman coding.  ICLR, 2016.
  * [26] Tianlong Chen, Yu Cheng, Zhe Gan, Lu Yuan, Lei Zhang, and Zhangyang Wang.  Chasing sparsity in vision transformers: An end-to-end exploration.  NeurIPS, 2021.
  * [27] Tianlong Chen, Zhenyu Zhang, Yu Cheng, Ahmed Awadallah, and Zhangyang Wang.  The principle of diversity: Training stronger vision transformers calls for reducing all levels of redundancy.  In CVPR, 2022.
  * [28] Shixing Yu, Tianlong Chen, Jiayi Shen, Huan Yuan, Jianchao Tan, Sen Yang, Ji Liu, and Zhangyang Wang.  Unified visual transformer compression.  In ICLR, 2022.
  * [29] Geoffrey Hinton, Oriol Vinyals, and Jeff Dean.  Distilling the knowledge in a neural network.  arXiv preprint arXiv:1503.02531, 2015.
  * [30] Xianing Chen, Qiong Cao, Yujie Zhong, Jing Zhang, Shenghua Gao, and Dacheng Tao.  Dearkd: Data-efficient early knowledge distillation for vision transformers.  In CVPR, 2022.
  * [31] Sheng Shen, Zhen Dong, Jiayu Ye, Linjian Ma, Zhewei Yao, Amir Gholami, Michael W Mahoney, and Kurt Keutzer.  Q-bert: Hessian based ultra low precision quantization of bert.  In AAAI, 2020.
  * [32] Mengshu Sun, Haoyu Ma, Guoliang Kang, Yifan Jiang, Tianlong Chen, Xiaolong Ma, Zhangyang Wang, and Yanzhi Wang.  Vaqf: Fully automatic software-hardware co-design framework for low-bit vision transformer.  arXiv preprint arXiv:2201.06618, 2022.
  * [33] Ze Liu, Yutong Lin, Yue Cao, Han Hu, Yixuan Wei, Zheng Zhang, Stephen Lin, and Baining Guo.  Swin transformer: Hierarchical vision transformer using shifted windows.  In ICCV, 2021.
  * [34] Chun-Fu Chen, Quanfu Fan, and Rameswar Panda.  Crossvit: Cross-attention multi-scale vision transformer for image classification.  ICCV, 2021.
  * [35] Li Yuan, Yunpeng Chen, Tao Wang, Weihao Yu, Yujun Shi, Zi-Hang Jiang, Francis EH Tay, Jiashi Feng, and Shuicheng Yan.  Tokens-to-token vit: Training vision transformers from scratch on imagenet.  In ICCV, 2021.
  * [36] Mathilde Caron, Hugo Touvron, Ishan Misra, Hervé Jégou, Julien Mairal, Piotr Bojanowski, and Armand Joulin.  Emerging properties in self-supervised vision transformers.  In ICCV, 2021.
  * [37] Michael Ryoo, AJ Piergiovanni, Anurag Arnab, Mostafa Dehghani, and Anelia Angelova.  Tokenlearner: Adaptive space-time tokenization for videos.  NeurIPS, 2021.
  * [38] Yongming Rao, Wenliang Zhao, Benlin Liu, Jiwen Lu, Jie Zhou, and Cho-Jui Hsieh.  Dynamicvit: Efficient vision transformers with dynamic token sparsification, 2021.
  * [39] Zhenglun Kong, Peiyan Dong, Xiaolong Ma, Xin Meng, Wei Niu, Mengshu Sun, Bin Ren, Minghai Qin, Hao Tang, and Yanzhi Wang.  Spvit: Enabling faster vision transformers via soft token pruning.  ECCV, 2022.
  * [40] Youwei Liang, Chongjian GE, Zhan Tong, Yibing Song, Jue Wang, and Pengtao Xie.  EVit: Expediting vision transformers via token reorganizations.  In ICLR, 2022.
  * [41] Lingchen Meng, Hengduo Li, Bor-Chun Chen, Shiyi Lan, Zuxuan Wu, Yu-Gang Jiang, and Ser-Nam Lim.  Adavit: Adaptive vision transformers for efficient image recognition.  In CVPR, 2022.
  * [42] Changqian Yu, Bin Xiao, Changxin Gao, Lu Yuan, Lei Zhang, Nong Sang, and Jingdong Wang.  Lite-hrnet: A lightweight high-resolution network.  In CVPR, 2021.
  * [43] Daniil Osokin.  Real-time 2d multi-person pose estimation on cpu: Lightweight openpose.  arXiv preprint arXiv:1811.12004, 2018.
  * [44] Christopher Neff, Aneri Sheth, Steven Furgurson, and Hamed Tabkhi.  Efficienthrnet: Efficient scaling for lightweight high-resolution multi-person pose estimation.  arXiv preprint arXiv:2007.08090, 2020.
  * [45] Minghao Yin, Zhuliang Yao, Yue Cao, Xiu Li, Zheng Zhang, Stephen Lin, and Han Hu.  Disentangled non-local neural networks, 2020.
  * [46] Zhisheng Lu, Juncheng Li, Hong Liu, Chaoyan Huang, Linlin Zhang, and Tieyong Zeng.  Transformer for single image super-resolution, Aug 2021.
  * [47] Iz Beltagy, MatthewE. Peters, and Arman Cohan.  Longformer: The long-document transformer, Apr 2020.
  * [48] Xizhou Zhu, Weijie Su, Lewei Lu, Bin Li, Xiaogang Wang, and Jifeng Dai.  Deformable detr: Deformable transformers for end-to-end object detection, Oct 2020.
  * [49] Zhenglun Kong, Haoyu Ma, Geng Yuan, Mengshu Sun, Yanyue Xie, Peiyan Dong, Xin Meng, Xuan Shen, Hao Tang, Minghai Qin, Tianlong Chen, Xiaolong Ma, Xiaohui Xie, Zhangyang Wang, and Yanzhi Wang.  Peeling the onion: Hierarchical reduction of data redundancy for efficient vision transformer training, Nov 2022.
  * [50] Zhenglun Kong, Peiyan Dong, Xiaolong Ma, Xin Meng, Mengshu Sun, Wei Niu, Xuan Shen, Geng Yuan, Bin Ren, Minghai Qin, Hao Tang, and Yanzhi Wang.  Spvit: Enabling faster vision transformers via soft token pruning, 2022\.
  * [51] Mykhaylo Andriluka, Leonid Pishchulin, Peter Gehler, and Bernt Schiele.  2d human pose estimation: New benchmark and state of the art analysis.  In CVPR, 2014.


# Quantum Computing for Power Flow Algorithms: Testing on real Quantum
Computers

Brynjar Sævarsson, Spyros Chatzivasileiadis, Hjörtur Jóhannsson, Jacob
Østergaard  Department of Wind and Energy Systems
Division of Electric Power and Energy
Technical University of Denmark
Kgs. Lyngby, Denmark
{brysa, spchatz, hjjo<EMAIL_ADDRESS>
###### Abstract

Quantum computing has the potential to solve many computational problems
exponentially faster than classical computers. The high shares of renewables
and the wide deployment of converter-interfaced resources require new tools
that shall drastically accelerate power system computations, including
optimization and security assessment, which can benefit from quantum
computing. To the best of our knowledge, this is the first paper that goes
beyond quantum computing simulations and performs an experimental application
of Quantum Computing for power systems on a real quantum computer. We use five
different quantum computers, apply the HHL quantum algorithm, and examine the
impact of current noisy quantum hardware on the accuracy and speed of an AC
power flow algorithm. We perform the same studies on a 3-bus and a 5-bus
system with real quantum computers to identify challenges and open research
questions related with the scalability of these algorithms.

###### Index Terms:

power grids, power system security, quantum computing, quantum power flow

##  I Introduction

The increasing penetration of distributed renewable energy sources (RES)
brings a number of challenges when it comes to ensuring a secure operation of
power systems. With traditional synchronous generation being replaced by
thousands of converter-interfaced generation sources, the complexity of the
power system greatly increases, its natural inertia decreases, and the
dynamics of the system change. This means traditional offline methods for
security assessment are expected to become insufficient as more detailed, and
computationally demanding, simulations are required to capture the faster
dynamics of RES. Additionally, the operating state of the system becomes
subject to the prevailing weather conditions, making it difficult to predict
hours ahead of time. The uncertainties introduced by RES can lead to
exponentially more contingency scenarios needed to be considered to ensure N-x
security. These challenges drive the need for new real-time and offline tools
for security analysis. Assessing large complex systems for millions of
contingency scenarios can be very computationally demanding and is currently
one of the major challenges that utilities are expected to face in the future.

Over the last few years there has been a great leap in the area of Quantum
Computing (QC), bringing us into the so called Noisy Intermediate-Scale
Quantum (NISQ) era [1] of quantum computing, where real quantum computers are
already containing over 100 qubits. According to roadmaps of vendors such as
IBM [2], there is ambitious development in building scalable quantum
computers, with the aim of creating large-scale and noise-free devices in the
near future. However, the current technology readiness level (TRL) of the
NISQ-era quantum computers is still rather low and only very basic
applications can be implemented. The use of quantum computations for power
system applications is very new. Quantum power flow algorithms have been
proposed to solve DC power flow [3] and AC power flow [4], and were tested
with simulated quantum computers. The use of QC for contingency assessment was
introduced in [5] and for performing EMT simulations in [6]. A common feature
of these papers is the use of the HHL quantum algorithm for solving linear
equations [7], as this promises exponential speedup compared to classical
methods.

To the best of our knowledge, this paper presents the first implementation of
an AC power flow on real Quantum Computers. The intention with this paper is
to explore the current capabilities of QC for power flow studies, which
challenges we need to address, and to investigate the foreseen future
capabilities and practical implications of QC for power systems. We test our
power flow algorithm for a simple 3 bus system on four of IBM’s publicly
available quantum computers to investigate the impact of noisy hardware. For a
5 bus system we use one of IBM’s larger quantum computers to investigate how
the method scales with increasing system size. This is important for future
implementations of full-scale power systems since the benefits of the
computational power of quantum computing will be more apparent for large power
systems. Larger systems exceed the current hardware capabilities of real
quantum computers and are not tested in this paper.

This paper is organized as follows. Section II provides an overview of quantum
computing and its potential for power systems. Section III describes the
quantum power flow method and its implementation. In Section IV, we describe
the simulation setup and the quantum hardware requirements. Section V provides
the results of the power flow on both real and simulated quantum computers. In
Section VI, we discuss potential issues with the scalability of the method.
Section VII concludes.

##  II The potential of quantum computing for power systems

The way quantum computers function is fundamentally different from classical
computers. Instead of using classical bits, which can only be in the states 0
and 1, quantum computers use quantum bits (qubits) which, in addition to being
in states 0 and 1, can form a linear combination of states, i.e. have a
certain probability to be measured either as 0 or 1. This state is referred to
as superposition. While in superposition, they can also be said to be in all
states simultaneously, i.e. both 0 and 1 or somewhere in-between, and only
when measured do they collapse to 1 or 0 with some probability. An important
feature of this is that the information contained in a quantum register grows
exponentially for every qubit which is added [8]. That is, one could say that
as a rule of thumb, instead of requiring $2^{n}$ bits for a computing process
with classical computing, quantum computing requires only $n$ qubits. For
example, if a process would require 1024 bits in a classical computer, a
quantum computer would aprroximately only need 10 qubits. Please bear in mind
that currently publicly available quantum computers do not exceed 7 qubits.
Currently existing Quantum Computers (but not publicly available), on the
other hand, go up to 100 qubits. That is a 30 orders of magnitude more
powerful computer, while we have increased the number of qubits by only 1
order of magnitude.

With qubits being able to be in superposition, i.e. 0 or 1 or somewhere in-
between, quantum states can also be used to evaluate functions for multiple
values simultaneously. In a power systems context, such a feature can become
extremely valuable for power flow analysis of multiple scenarios. Instead of
sequentially solving multiple power flows to assess a large number of
generation and demand scenarios, theory suggests that future QC infrastructure
can deliver the solution for many scenarios simultaneously, drastically
accelerating e.g. the security assessment. This, and other special features of
QC, allows them to solve certain problems exponentially faster than classical
computers. Theoretically, Quantum Computers can solve all the same problems as
classical computers; not all problems, though, allow for exponential speedups
and in some cases QC are even slower than classical computers. QC is therefore
seen as a supplement to classical computing for solving specific problems in a
hybrid quantum-classical computation.

In the next sections, we present the formulation of the quantum power flow
algorithm and the results we obtain after testing it on real quantum
computers. Considering that the quantum computers that were accessible at the
time of this paper only have 7 qubits, the application of our quantum power
flow algorithm on real quantum computers is limited to a 5-bus system.

##  III Quantum Power Flow method

The implementation of the quantum power flow method we follow in this paper is
based on the Fast Decoupled Load Flow (FDLF) method, which is a commonly used
adaptation of the Newton Raphson Power Flow (NRPF). It exploits inter-
dependence between $P-\theta$ and $Q-|V|$ to create two constant Jacobian sub-
matrices $B^{{}^{\prime}}$ and $B^{{}^{\prime\prime}}$. This replaces one of
the most computationally heavy parts of NRPF i.e. updating the Jacobian matrix
in every iteration. FDLF converges to the same solution as NRPF since the
mismatch functions are the same, but it usually requires more iterations than
NRPF [9]. On the other hand, each iteration of FDLF is much faster, meaning
that the overall computation time of FDLF can often be quite faster,
especially for large systems. In each iteration, we solve equations (1)-(3)
for $\Delta\theta$ and $\Delta V_{m}$ until the norm of $\Delta P$ and $\Delta
Q$ is less than the chosen tolerance $\xi$.

$\displaystyle\Delta
S=(S_{bus}-\bm{V}\circ\overline{(Ybus\cdot\bm{V})})\oslash V_{m}$ (1)
$\displaystyle\Delta P=\Re(\Delta S_{pv+pq})=B^{{}^{\prime}}\Delta\theta$ (2)
$\displaystyle\Delta Q=\Im(\Delta S_{pq})=B^{{}^{\prime\prime}}\Delta V_{m}$
(3)

In (1)-(3), $\bm{V}=V_{m}\angle\theta$ is a vector collecting the voltage
phasors for all buses, $\circ$ and $\oslash$ denote the element-wise product
and division respectively, and $pv$ and $pq$ are indices corresponding to the
set of PV and PQ buses respectively. The structure of FDLF makes it very
suitable for hybrid classical-quantum computing since we can easily replace
the numerical algorithm solving equations (2) and (3) with a method which is
capable of running on a quantum computer. HHL is probably the most popular
quantum algorithm for solving a set of linear equations at the moment, with a
runtime of $\mathcal{O}(log(N)\frac{s^{2}\kappa^{2}}{\epsilon})$ [7]; here, N
is the number of equations, $s$ is the sparsity and $\kappa$ the condition
number of the system matrix, and $\epsilon$ is the solution accuracy. In
theory, HHL can achieve an exponential speedup over the fastest classical
algorithm, the conjugate gradient method, which has a runtime of
$\mathcal{O}(Ns\kappa log(\frac{1}{\epsilon}))$.

A system of linear equations is usually given in the form:

$Ax=b,$ (4)

where, given a known matrix $A$ and vector $b$, we solve for $x$. For quantum
systems we use the Dirac (or bra-ket) notation, where we denote a vector with
a ket which has the form $\bm{\ket{v}}$ and represents the state of a quantum
system. We also denote a bra with the form $\bm{\bra{v}}$, where bra is the
conjugate transpose of ket. Quantum vectors are often very sparse, i.e. with
most amplitudes equal to zero, and this notation makes it much simpler to
represent values of interest. A single qubit system can be described as a
combination of the two states $\ket{0}$ and $\ket{1}$:

$\ket{q}=\alpha\ket{0}+\beta\ket{1},$ (5)

where the probability of it being measured as 0 is $|\alpha|^{2}$ and the
probability of it being measured as 1 is $|\beta|^{2}$. The combined
probability amplitudes must satisfy that $|\alpha|^{2}+|\beta|^{2}=1$. As also
mentioned earlier, this is one of the key strengths of quantum computing: as
qubits can be in both states at the same time, they can represent multiple
variables. For a single qubit, for example, we can associate one variable with
state 0 and one variable with state 1. The probability of the qubit to be
measured as 1 corresponds to the amplitude of that variable. Let us assume we
encode $\Delta\theta_{1}$ to the qubit state $\ket{0}$ and $\Delta\theta_{2}$
to the qubit state $\ket{1}$. The value of $\Delta\theta_{1}$ is then equal to
$C\cdot\alpha$, where $C$ is a scaling factor. As we will see in Section V and
(17), the scaling factor $C$ depends on how we have normalized our input to be
between 0 and 1.

Going back to the solution of a linear system with the HHL algorithm, we
rescale the system in (4) by normalizing $x$ and $b$, and we can then map them
to the quantum states $\ket{x}$ and $\ket{b}$.

$A\ket{x}=\ket{b}$ (6)

Portraying equation (2) in this form we get:

$B^{{}^{\prime}}\ket{\Delta\theta}=\ket{\Delta P}$ (7)

If we express the right hand side in the eigenbasis of the $B^{{}^{\prime}}$
matrix, we can write it as:

$\displaystyle\ket{\Delta P}=\sum_{j=0}^{N-1}\Delta
P_{j}\ket{u^{{}^{\prime}}_{j}}$ (8)

The HHL method requires matrix $A$ to be Hermitian, that is equal to its
conjugate transpose. For non-hermitian matrices, a hermitian matrix can be
constructed from A such that:

$\begin{bmatrix}0&A\\\ A^{\dagger}&0\end{bmatrix}\begin{bmatrix}0\\\
x\end{bmatrix}=\begin{bmatrix}b\\\ 0\end{bmatrix}$ (9)

However, $B^{{}^{\prime}}$ and $B^{{}^{\prime\prime}}$ of FDLF happen to be
hermitian. And since the matrices are Hermitian, they have a spectral
decomposition:

$\displaystyle
B^{{}^{\prime}}=\sum_{j=0}^{N-1}\lambda^{{}^{\prime}}_{j}\ket{u^{{}^{\prime}}_{j}}\bra{u^{{}^{\prime}}_{j}}$
(10) $\displaystyle
B^{{}^{\prime}-1}=\sum_{j=0}^{N-1}\lambda^{{}^{\prime}-1}_{j}\ket{u^{{}^{\prime}}_{j}}\bra{u^{{}^{\prime}}_{j}}$
(11)

Where $\lambda_{j}^{{}^{\prime}}$ and $u_{j}^{{}^{\prime}}$ are the $j_{th}$
eigenvalue and eigenvector of the $B^{{}^{\prime}}$ matrix. Putting (7), (8)
and (11) together we then get:

$\displaystyle\ket{\Delta\theta}=B^{{}^{\prime}-1}\ket{\Delta
P}=\sum_{j=0}^{N-1}\lambda^{{}^{\prime}-1}_{j}\Delta
P_{j}\ket{u^{{}^{\prime}}_{j}}\bra{u^{{}^{\prime}}_{j}}\ket{u^{{}^{\prime}}_{j}}$
(12)

As the scalar product $\bra{u^{{}^{\prime}}_{j}}\ket{u^{{}^{\prime}}_{j}}=1$,
(12) simplifies to:

$\displaystyle\ket{\Delta\theta}=\sum_{j=0}^{N-1}\lambda^{{}^{\prime}-1}_{j}\Delta
P_{j}\ket{u^{{}^{\prime}}_{j}}$ (13)

And following the same procedure for $\Delta V_{m}$ and
$B^{{}^{\prime\prime}}$ we get:

$\displaystyle\ket{\Delta V_{m}}=B^{{}^{\prime\prime}-1}\ket{\Delta
Q}=\sum_{j=0}^{N-1}\lambda^{{}^{\prime\prime}-1}_{j}\Delta
Q_{j}\ket{u^{{}^{\prime\prime}}_{j}}$ (14)

To solve (13) a quantum circuit is implemented as shown in Fig. 1. This
circuit is for solving a 2x2 system of equations, i.e., $N=2$. It consists of
three registers, $nb,nl,na$, which are all initialized as $\ket{0}$. $nb$ is
the data register where we load the data vector $b$ (here $b\equiv\Delta P$
onto $log_{2}(N)=1$ qubits):

$\ket{0}\mapsto\ket{\Delta P}$ (15)

Next, we apply quantum phase estimation to estimate the eigenvalues of the $B$
matrix. The $nl$ register consists of 3 qubits where we store the
approximation of the eigenvalues $\ket{\tilde{\lambda_{j}}}$. The more qubits
used in $nl$, the more accurately we can estimate the eigenvalues of $B$. The
last register, $na$, consists of a single auxiliary qubit. The output
auxiliary qubit indicates if we can obtain a binary estimation of the
eigenvalues and the data vector, i.e. if we can trust their numerical result
when we measure them, or not. To get the output of the auxiliary qubit, we
perform a rotation conditioned on $\ket{\tilde{\lambda_{j}}}$, and finally an
inverse phase estimation. When we measure the output of the auxiliary qubit as
1, then the registers $nb,nl$ are in the post measurement state; this means
that the phase estimation outputs a binary estimation of the eigenvalues of
the B matrix, and after an inverse phase estimation, the $nb$ register
contains the solution $\ket{\Delta\theta}$. If the auxiliary qubit is 0, then
the states are considered to contain no useful information and are, therefore,
discarded. For more information about the steps of the HHL algorithm, the
interested reader can refer to [7].

Figure 1: Quantum HHL circuit used for solving a Quantum Power Flow (QPF) on a
3-bus system. The blocks shown in the circuit are: Preparation of the data
vector (Load b), Quantum Phase Estimation (QPE), conditioned rotation (1/x)
and inverse QPE (QPE_dg). Finally, each qubit is measured and the result is
stored in a classical 5 bit register where “meas 0” is the least significant
bit and “meas 4” is the most significant bit.

The implemented Quantum Power Flow (QPF) algorithm is shown in Algorithm 1.

Algorithm 1 Quantum Power Flow

Input $Y_{bus},B^{{}^{\prime}},B^{{}^{\prime\prime}},V_{0},S_{bus}$

Update $\Delta P$, $\Delta Q$

while $\lVert\Delta P\rVert\geq\xi$ AND $\lVert\Delta Q\rVert\geq\xi$ do

$\Delta\theta=HHL(B^{{}^{\prime}},\Delta P)$

$\Delta V_{m}=HHL(B^{{}^{\prime\prime}},\Delta Q)$

Update $V_{m}\angle\theta$

Update $\Delta P$, $\Delta Q$

end while

return $V_{m}\angle\theta$

##  IV Simulations

The capabilities of both real and simulated Quantum Computers are still
limited. Some simplifications must therefore be made to run the power flow
algorithm on a real quantum computer. If each of the eigenvalues of the system
matrix cannot be exactly represented by the same number of classical bits as
there are qubits in the $nl$ register, the quantum phase estimation will only
give an approximation of the eigenvalues $\lambda_{j}$ of the $B^{\prime}$ and
$B^{\prime\prime}$ matrices; this introduces an error in the results. In order
to investigate the performance of real quantum computers, at this stage, we
select the parameters of the test systems so that each of the eigenvalues of
the $B$ matrices can be closely represented by 3 bits. This will correspond to
the need of $nl=3$ qubits in our quantum circuit. As our focus is to test a
simple system on a real QC, for the purpose of the investigations in this
paper we also consider only slack and PQ buses, and no phase-shifting
transformers. This means that $B^{{}^{\prime}}=B^{{}^{\prime\prime}}$, and we
can use the same HHL circuit for both (13) and (14) only with a different
vector $b$ loaded onto the $nb$ register (i.e., $\Delta\theta$ or $\Delta V$).
The two power systems considered for this paper are shown in Fig. 2 and the
parameters of the 3-bus system are shown in Table I.

132$G_{1}$$G_{2}$ (a) Three bus system

21345$G_{2}$$G_{1}$ (b) Five bus system

Figure 2: Systems for testing Quantum Power Flow (QPF). TABLE I: Three bus
system parameters

Bus | Type | $\mathbf{P_{MW}}$ | $\mathbf{Q_{MVAr}}$ | $\mathbf{V_{pu}}$ | $\bm{\theta^{\circ}}$
---|---|---|---|---|---
1 | Slack | 5 | 6.9 | 1.03 | 0
2 | PQ | 10 | 0 | 1.012 | 2.109
3 | PQ | -15 | -5 | 0.985 | -4.996

Line | From | To | $\mathbf{R_{pu}}$ | $\mathbf{X_{pu}}$
---|---|---|---|---
1 | 1 | 2 | 0 | 1j
2 | 1 | 3 | 0 | 1j
3 | 2 | 3 | 0 | 2j

The parameters in Table I give the $B$ matrices:

$\bm{B^{{}^{\prime}}}=\bm{B^{{}^{\prime\prime}}}=\begin{bmatrix}-1.5&0.5\\\
0.5&-1.5\\\ \end{bmatrix}$ (16)

The eigenvalues of matrices $\bm{B^{{}^{\prime}}}$ and
$\bm{B^{{}^{\prime\prime}}}$ are $\\{-1,-2\\}$.

###  IV-A Quantum hardware requirements

The number of qubits (circuit width) required for the HHL circuit depends on
the required accuracy of the phase estimation and the number of variables in
the data vector (plus the auxiliary qubit that we always need to indicate if
we have obtained a binary estimation or not). For example, the 3-bus system
has 2 variables for $\Delta P$ which are encoded in 1 qubit
($\\#qubits=log_{2}N$, where $N$ is the number of variables, as mentioned in
Section II) and we use a 2-bit estimation of the eigenvalues of
$B^{{}^{\prime}}$. Another qubit is needed for the phase estimation in case
the system matrix has negative eigenvalues. Including the auxiliary qubit,
this gives a total of 5 qubits required for the 3-bus system. The qubit for
negative eigenvalues could be eliminated if we can assume all eigenvalues have
the same sign. This is, however, not always the case for larger power systems
which could for example include capacitive branches.

We shall note here that the true benefits of Quantum Computing emerge in large
systems, since the number of required qubits have a sublinear increase with
larger system sizes (logarithmic to be precise). This means that we can
represent very large systems with only a few qubits. For example, applying the
$log_{2}(N)$ rule for the required number of qubits, we can represent the data
vector of a 500’000-bus system with only 20 qubits
($2^{20}=1$’$048$’$576\approx 2\times 500$’$000$ values).

Currently, however, the major issue with scalability is not as much the
circuit width, i.e. the number of qubits, but the circuit depth; that is, the
number of gates in the quantum circuit, as shown in Fig. 1 . The current
general implementation of the HHL algorithm uses methods which scale
exponentially with the number of qubits, and, so far, no algorithm exists for
preparing the HHL circuit with polynomially increasing resources for an
arbitrary matrix. Table II shows how the circuit sizes grow with the matrix
size.

TABLE II: Quantum Circuit Sizes | 3 bus | 5 bus | 9 bus | 17 bus
---|---|---|---|---
$B^{{}^{\prime}}$ matrix size | 2x2 | 4x4 | 8x8 | 16x16
Circuit Width | 5 | 7 | 9 | 11
Circuit Depth | 336 | 3528 | 76876 | 802737
CNOT Gates | 108 | 1181 | 28956 | 298594

##  V Simulation results

###  V-A 3-bus system on real Quantum Computers

The power flow application is implemented with IBM’s qiskit (0.34.1). The
3-bus system is tested on 4 of IBM’s open access quantum computers [10] listed
in Table III: ibmq_lima, ibmq_belem, ibmq_quito, and ibmq_bogota. The quantum
computers listed in Table III have different sizes, configurations and error
rates. The Quantum Volume (QV) is a measure of the device’s performance
regardless of the underlying technology or the number of qubits.

TABLE III: Quantum Computers used for testing Quantum Power Flow (values taken on 23.02.2022) QC | Qubits | QV | avg. CNOT error
---|---|---|---
ibmq_lima | 5 | 8 | 9.996e-3
ibmq_belem | 5 | 16 | 1.363e-2
ibmq_quito | 5 | 16 | 1.135e-2
ibmq_bogota | 5 | 32 | 1.206e-2
ibm_perth | 7 | 32 | 1.006e-2

The results of the 3-bus power flow on ibmq_quito are shown for each iteration
of the power flow calculation in Table IV. The rest of the quantum computers
we used, i.e. ibmq_lima, ibmq_belem, and ibmq_bogota give similar results. The
Quantum Power Flow converges to the same solution as the classical FDLF but,
due to the noisy hardware available at the moment, it requires a much larger
number of iterations. With the tolerance set to $10^{-5}$, the classical
method converges in 5 iterations while the Quantum Power Flow (QPF) takes
around 32-38 iterations.

The available quantum software (qiskit 0.34.1) does not allow us to implement
if-loops, for-loops, and while-loops within a real quantum algorithm at the
moment (March 2022). As a result, implementing a numerical solution algorithm
requires us to extract the value of the quantum data vector in every
iteration, perform the logical operations (e.g. is $\Delta P\geq\xi$, see
Algorithm 1) and plug $\Delta P$ back in the quantum circuit again for the
next iteration. Considering that the data vector is in a quantum state, we
need to ”probe” $\Delta P$ multiple times to measure the probability
distribution of the solution. This requires time of at least $\mathcal{O}(N)$
(see also the Discussion in Section VI). Obviously, this slows down our
quantum algorithm considerably at the moment. The next qiskit version, which
is expected around August 2022, is expected to allow for extended
functionality; this will unleash a wide range of opportunities, including if-
clauses, for-loops, and while-loops inside the quantum algorithm which could
eliminate the need for measuring the data vector in each iteration of the QPF.

Going back to the currently possible implementation of our QPF, where the data
vector $\Delta P$ needs to be measured in every iteration, Fig. 3 shows the
measurement from the first iteration with and without noise. The derivation of
(13) assumes perfect phase estimation. Based on the HHL circuit, that means
that qubits in the $nl$ register will be in the state $\ket{0}$ after the
inverse phase estimation, and the result state is conditioned on seeing the
auxiliary qubit in register $na$, as 1. In that case, the normalized values of
$\Delta\theta$ can be read from the $nb$ register, i.e. the states “10000” and
“10001” in the histogram in Fig. 3.

(a) Without noise

(b) With noise

Figure 3: Histogram showing measurements of $\ket{\Delta\theta}$ in the first
iteration. The values of $\Delta\theta$ are estimated from the probability
distribution of the data read in register $nb$ (see Fig. 1), when the
auxiliary bit $na=1$ i.e. from the states “10000” and “10001”

The quantum circuit in Fig. 1 is executed 1024 times to retrieve a probability
distribution of the quantum states as shown in Fig. 3. The probabilities shown
in the histogram sum up to 1 where the amplitudes of the quantum state
coefficients satisfy $\left|a\right|^{2}+\left|b\right|^{2}+...=1$. To get an
estimate of the amplitudes of $\Delta\theta$, as we also explained in Section
III, we must therefore take the square root of the probability and scale it
with the norm of the input $\Delta P$.

$\Delta\theta\approx\lVert\Delta P\rVert*\sqrt{probability}$ (17)

For the first iteration of the noise free simulation, this gives:

$\begin{bmatrix}\Delta\theta_{2}\\\ \Delta\theta_{3}\\\ \end{bmatrix}\approx
0.1803*\begin{bmatrix}\sqrt{0.045}\\\ \sqrt{0.229}\\\
\end{bmatrix}\approx\begin{bmatrix}0.03824\\\ 0.08628\\\ \end{bmatrix}$ (18)

Where $\lVert\Delta P\rVert=0.1803$. We can see that the result of the noise-
free quantum simulation is very close to the result for the first iteration of
the classical computing algorithm, which is $\Delta\theta=[0.0375,0.0875]$
rad. On the contrary, for noisy quantum computations, which are based on real
quantum computers, the output is $\Delta\theta=[0.0526,0.0748]$ rad after the
first iteration. This is the reason why noise-free quantum computing converges
almost equally fast as classical computing (same number of iterations, see
Fig. 4), while the algorithms implemented on actual quantum computers require
at the moment a significantly larger number of iterations.

TABLE IV: Iterations of Classical FDLF vs QPF on ibmq_quito

FDLF iteration | $V_{2}$ | $\theta_{2}$ | $V_{3}$ | $\theta_{3}$
---|---|---|---|---
1 | 1.01750000 | 2.14859173 | 0.99250000 | -5.01338071
2 | 1.01198104 | 2.08682333 | 0.98513520 | -4.94593291
3 | 1.01210909 | 2.10823443 | 0.98510531 | -4.99645096
4 | 1.01200181 | 2.10832465 | 0.98496362 | -4.99531212
5 | 1.01200181 | 2.10872362 | 0.98496362 | -4.99629051

QPF iteration | $V_{2}$ | $\theta_{2}$ | $V_{3}$ | $\theta_{3}$
---|---|---|---|---
1 | 1.00602829 | 2.05419579 | 0.99405120 | -1.77164967
2 | 1.00937372 | 3.12147980 | 0.99099606 | -2.73581045
3 | 1.01129425 | 2.50579482 | 0.04373437 | -3.30344432
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ |
34 | 1.01200263 | 2.10882674 | 0.98496970 | -4.99558292
35 | 1.01200492 | 2.10901696 | 0.98496697 | -4.99577486

A comparison between the convergence of the classical method, the 4 physical
QC devices and a simulated noise-free QC is shown in Fig. 4. It shows how a
simulated noise-free QC gives close to identical results with a classical
computer, while the real noisy QCs converge much slower. As the reader
observes, there is no notable difference between the implementations on the 4
real QCs.

Figure 4: Power flow convergence of the 3-bus system

###  V-B 5-bus system

We also test our QPF on a 5-bus system. Our goal is to investigate how current
noisy hardware impact the Quantum Power Flow solution as we scale to larger
systems. As shown in Table II, a QPF for a 5-bus system requires 7 qubits and
a much larger number of gates. We therefore use the 7 qubit machine
”ibm_perth”, shown in Table III, to test the 5 bus system.

Similar to the 3-bus system, the line parameters of the 5 bus system are
chosen so the eigenvalues of the $\mathbf{B}$ matrices can be closely
represented by 3 bits, as our primary focus is to examine the impact of the
noisy quantum hardware. The $\mathbf{B}$ matrices are presented in (19). Fully
acknowledging though that the characteristics of real power systems can vary
significantly, in our experiments we also perturb the system characteristics
(and the $\mathbf{B}$ matrices shown in (19)) to explore cases where
$\mathbf{B}$ matrices result to eigenvalues that cannot be accurately
represented by 3 bits (see Fig. 5).

$\bm{B^{{}^{\prime}}}=\bm{B^{{}^{\prime\prime}}}=\begin{bmatrix}-4&0.03&0&0\\\
0.03&-3&0.02&0\\\ 0&0.02&-1.55&0.5\\\ 0&0&0.5&-1.45\\\ \end{bmatrix}$ (19)

The eigenvalues of $\bm{B^{{}^{\prime}}}$ and $\bm{B^{{}^{\prime\prime}}}$ in
this case are $\\{-1,-2,-3,-4\\}$.

Figure 5 shows the iterations of the 5-bus power flow. As with the 3-bus
system, it converges to the correct solution and the noise-free simulation is
very close to the classical one. When we alter the line parameters, so that
the eigenvalues of the $\mathbf{B}$ matrices cannot be represented exactly by
3 bits, we observe that this introduces a small delay in the convergence (red
line in Fig. 5) and if the eigenvalues vary greatly then the QPF will not
converge at all without using more qubits. The major difference in convergence
speed emerges when we add the noise characteristics in the computation (green
line in Fig. 5).

As before, the presence of noise in the real hardware greatly impacts the
number of iterations required. The larger number of gates required for the
circuit compared to the 3-bus system also means greater impact from the noise,
as each gate adds noise to the result. The 5-bus system therefore requires
much larger number of iterations to converge on real hardware. Our simulations
show that unless we develop low-noise quantum hardware, it would be
challenging to perform computations for larger systems.

Figure 5: Power flow convergence of the 5-bus system

Based on the findings from our experiments, in the next section, we discuss
opportunities and challenges arising from the development of quantum
algorithms for power system applications.

##  VI Discussion

The result given by the HHL algorithm is encoded as a quantum state. If we
want to extract the full state vector from the quantum state we need time of
at least $\mathcal{O}(N)$. Therefore, if we do this in every iteration of our
power flow we lose any benefit of the quantum speedup [11]. For an HHL-based
AC QPF to achieve quantum advantage, the whole iteration process would need to
be performed by the QC, so that we only extract the final result. To do that,
we need some form of quantum memory in order to be able to perform logical
operations such as if-then-else clauses, for-loops and while-loops. The
quantum memory will allow us to store the quantum state between iterations so
we can update our power mismatches and perform the logical operations (e.g. is
$\Delta P\geq\xi$, see Algorithm 1) inside the QC algorithm, without requiring
us to extract e.g. the $\Delta P$ in every iteration, check if it is greater
than the defined tolerance $\xi$ and plug it back to the QC algorithm. Note
here that every time we have to extract a quantum state, e.g. determine the
value of $\Delta P$ outside the QC algorithm, we need to ”probe” $\Delta P$
multiple times to measure the probability distribution of the solution. This
requires time of at least $\mathcal{O}(N)$, which slows down our quantum
algorithm considerably at the moment. The functionality that will allow the
logical operations inside the QC algorithms is not yet available in current
real QC machines. The next release of qiskit, probably in August 2022, is
expected to include some form of looping capabilities which will greatly
improve what we are capable of doing with our quantum algorithms.

Even when we have some form of quantum memory, the process of extracting the
final results can still reduce the speedup of QC, as we need to probe the
output register several times to extract accurately the probability
distributions (see e.g. Fig. 3 for an example of the extracted histogram for
the $\Delta P$ output of the 3-bus system). Still, however, it is possible to
obtain some limited statistical information, such as the presence of very high
values, from the HHL algorithm without extracting the full solution vector,
i.e. without requiring us to probe the output register too many times. This
could be useful for contingency assessment where we have a large number of
scenarios and only want to identify the ones that are unstable or overloaded.
More research is necessary to identify the opportunities of such approaches.

A further discussion point that is often raised relates to the number of
qubits required for larger systems. Indeed, a larger number of qubits will be
required to implement a more general QPF method capable of accurately
representing both real and complex eigenvalues of the system matrices of real
power systems. However, the number of required quantum states scale only
logaritmically with the number of variables, i.e. $\\#qubits=\log_{2}(N)$,
where $N$ is the number of variables. For example, as mentioned in Section IV,
applying the $log_{2}(N)$ rule for the required number of qubits, we can
represent the data vector of a 500’000-bus system with only 20 qubits
($2^{20}=1$’$048$’$576\approx 2\times 500$’$000$ values). Same is the case for
the accuracy of the eigenvalues of the system matrices. Every qubit we add
increases the accuracy of the eigenvalues by one order of magnitude. With the
rapid development in quantum hardware, we expect that this will not be the
limiting factor for the QPF.

For the HHL algorithm to be efficient, the B matrix should be well-
conditioned. As B is Hermitian, the condition number is given by the ratio
between the largest and smallest eigenvalue
$\kappa=\frac{|\lambda_{max}|}{|\lambda_{min}|}$. If $\kappa$ grows
significantly with the size of B, then the exponential speedup of HHL is lost
[11]. As shown in Fig. 6, the condition number seems to grow exponentially for
larger test cases. This means that there is a need to design a preprocessing
procedure for the system matrices to reduce their condition number while
maintaining their characteristics, before plugging them in a Quantum Computer,
if we wish to maintain the quantum advantage. Further research is required on
how to achieve that in order to enable efficient quantum computations.

Figure 6: Condition number $\kappa$ of $B^{{}^{\prime}}$ for 6 IEEE test
systems (3 - 300 buses) and the 2736-bus Polish system.

The quantum power flow converges to the correct solution even with current
noisy hardware. However, the high number of iterations required could make it
challenging to achieve quantum advantage while quantum computers have a high
level of noise. We expect that applications of quantum computing will become
attractive when the noise levels decrease. A growing number of scientists and
engineers are currently working on reducing the noise level, as this will
enable a wide range of new applications. This is expected to happen in the
following years. Meanwhile, power system researchers can focus on inventing
the methods that can adapt existing algorithms and processes to take advantage
of the emerging quantum computing capabilities.

The implementation of QPF, as presented in this paper, aims to adapt the
classical FDLF method to run on a quantum computer. As we showed, there is a
number of challenges with current available hardware and the lack of quantum
memory. However, since quantum computers are fundamentally different from
classical computers, future implementations of QPF might even take a different
approach in order to better utilize the unique capabilities of QC of
evaluating multiple values simultaneously and possibly avoid the iterative
process altogether. We identify here two key opportunities. First, the ability
to drastically accelerate approaches that require some form of probability
distribution, e.g. Monte-Carlo methods. The fact that quantum states can
encode a probability distribution could for example significantly accelerate
the computation of the mean and value-at-risk [12]. This can help accelerate,
for example, N-1 computations and help assess whether a probability
distribution of generation profiles (due to e.g. uncertain wind power and
solar PV infeed) leads to N-1 violations. We plan to investigate this in our
future work. Another opportunity relates to the Power Flow computation itself,
which belongs to the wider family of numerical algorithms for non-linear
systems. Designing a method that can efficiently implement the equivalent of a
Newton-Raphson algorithm in a Quantum Computing environment can unleash
opportunities not only for power systems for a very large family of problems
in analysis, optimization, and control. Intensive research efforts are
required towards these directions to uncover the still unexplored potential of
quantum computing.

Last but not least, the new capabilities that QC brings along require us to
think differently. QC has been shown to solve problems that have so far been
impossible with classical computing; such an example is the creation and
observation of a new phase of matter, popularly known as time crystals [13].
Pursuing an out-of-the-box thinking, where instead of replicating existing
power system algorithms from the classical computing domain to the QC domain,
could uncover opportunities for power system computation, assessment, and
control that have so far been impossible with classical computing. Here, we
need to stress that as QC researchers say, we talk about “quantum advantage”
and not about “quantum supremacy”. Our goal is not to completely replace
classical computing and high performance computing (HPC). On the contrary, our
goal shall be to identify which are the processes that Quantum Computing can
do well and deliver new capabilities, by complementing the strengths offered
by HPC.

##  VII Conclusion

In this paper, we have successfully implemented and tested, for the first time
to our knowledge, a quantum AC power flow application on real Noisy-
Intermediate-Scale-Quantum-era (NISQ-era) quantum computers. We have shown
that current hardware is capable of performing a power flow for small test
systems, but scalability is currently a major issue.

There are still several challenges to be resolved before practical quantum
power flow applications can achieve quantum advantage. But this also creates
room for major opportunities. In our paper, we try to identify potential
directions for future research that can reap the benefits of quantum
computing. With further development and increasing capabilities of quantum
computers, quantum applications for power systems could become extremely
useful for future power system analysis, control, and optimization. An
increased focus by the power system researchers on this area can lead to the
development of the necessary quantum computing methods that can exploit these
benefits as soon as they become available.

##  VIII Acknowledgements

The authors would like to thank Jan Lillelund, CTO of IBM Denmark, and Gianni
del Bimbo, Head of Quantum Engineering, and Cristina Sanz-Fernandez, Quantum
Software Developer, from Multiverse Computing for the insightful discussions.

## References

  * [1] J. Preskill, “Quantum computing in the nisq era and beyond,” _Quantum_ , vol. 2, p. 79, Aug 2018. [Online]. Available: http://dx.doi.org/10.22331/q-2018-08-06-79
  * [2] “Ibm’s roadmap for scaling quantum technology,” https://research.ibm.com/blog/ibm-quantum-roadmap, 2020.
  * [3] R. Eskandarpour, K. Ghosh, A. Khodaei, and A. Paaso, “Experimental quantum computing to solve network dc power flow problem,” 2021.
  * [4] F. Feng, Y. Zhou, and P. Zhang, “Quantum power flow,” 2021.
  * [5] R. Eskandarpour, P. Gokhale, A. Khodaei, F. T. Chong, A. Passo, and S. Bahramirad, “Quantum computing for enhancing grid security,” _IEEE Transactions on Power Systems_ , vol. 35, no. 5, pp. 4135–4137, 2020.
  * [6] Y. Zhou, F. Feng, and P. Zhang, “Quantum electromagnetic transients program,” _IEEE Transactions on Power Systems_ , vol. 36, no. 4, pp. 3813–3816, 2021\.
  * [7] A. W. Harrow, A. Hassidim, and S. Lloyd, “Quantum algorithm for linear systems of equations,” _Phys. Rev. Lett._ , vol. 103, p. 150502, Oct 2009. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevLett.103.150502
  * [8] M. A. Nielsen and I. L. Chuang, _Quantum Computation and Quantum Information: 10th Anniversary Edition_ , 10th ed. USA: Cambridge University Press, 2011.
  * [9] B. Stott and O. Alsac, “Fast decoupled load flow,” _IEEE Transactions on Power Apparatus and Systems_ , vol. PAS-93, no. 3, pp. 859–869, 1974.
  * [10] “Ibm quantum.” https://quantum-computing.ibm.com/, 2022.
  * [11] S. Aaronson, “Read the fine print,” _Nature Physics_ , vol. 11, no. 4, pp. 291–293, Apr 2015. [Online]. Available: https://doi.org/10.1038/nphys3272
  * [12] C. Sanz-Fernandez, R. Hernandez, C. D. Marciniak, I. Pogorelov, T. Monz, F. Benfenati, S. Mugel, and R. Orus, “Quantum portfolio value forecasting,” 2021\. [Online]. Available: https://arxiv.org/abs/2111.14970
  * [13] X. Mi, M. Ippoliti, and C. e. a. Quintana, “Time-crystalline eigenstate order on a quantum processor,” _Nature_ , vol. 601, p. 531–536, 2022.


# Phase retrieval with physics informed zero-shot learning

Sanjeev Kumar

###### Abstract

Phase can be reliably estimated from a single diffracted intensity image, if a
faithful prior information about the object is available. Examples include
amplitude bounds, object support, sparsity in the spatial or a transform
domain, deep image prior and the prior learnt from the labelled datasets by a
deep neural network. Deep learning facilitates state of art reconstruction
quality but requires a large labelled dataset (ground truth-measurement pair
acquired in the same experimental conditions) for training. To alleviate this
data requirement problem, this letter proposes a zero-shot learning method.
The letter demonstrates that the object-prior learnt by a deep neural network
while being trained for a denoising task can also be utilized for the phase
retrieval, if the diffraction physics is effectively enforced on the network
output. The letter additionally demonstrates that the incorporation of total
variation in the proposed zero-shot framework facilitates the reconstruction
of similar quality in lesser time ( e.g. $\sim$8.5 fold, for a test reported
in this letter).

††journal: ol

Phase retrieval from a single diffracted intensity image is an ill-posed
inverse problem which can be solved using an iterative algorithm or a pre-
trained deep neural network, in coherent or partially coherent illuminations
[1, 2]. It enables imaging of transparent or nearly-transparent samples,
commonly referred as weakly scattering samples, which otherwise do not produce
sufficient contrast in the intensity images. Phase imaging is of utter
importance in the live-cell monitoring and histopathological microscopy of
unstained tissue samples.

As any other ill-posed inverse problem, the choice of object-prior plays an
integral role in determining the quality of solution obtained in a phase
retrieval problem. In the earlier iterative methods like the Gerchberg-
Saxton’s method [3] and the Fienup’s method [1], the amplitude constraints [4]
or the object supports (the presence of features only in parts of the field of
view) [5] were chosen as the priors. Sparsity in the total-variation (TV)
domain, or in a wavelet domain or in a domain learnt from data (dictionary
learning) were established as excellent hand-crafted object-priors in the
image restoration problems like degradation and denoising [6, 7]. Katkovnik et
al, 2011 used sparse modelling of both the amplitude and the phase for complex
object estimation [8]. Some other hand-crafted priors were used by Jolivet et
al, 2018 in their optimization based reconstruction framework, for in-line
holography [9]. They included bounds on the transmittance values, maximum and
minimum phase and spatial smoothness.

Rivenson et al, 2018 reported phase imaging in holography through end-to-end
learning with a deep convolutional neural network (deep CNN) [2, 10, 11]. This
network was trained with ground truth images obtained by multi-height phase
recovery method [12, 13] (8 multi-height holograms per image were used). Wang
et al, 2019 reported another deep network for the same task but they trained
their network with ground truth phase images reconstructed by the convolution
method [14]. In 2018, Ulyanov et al reported the use of untrained deep CNN as
a handcrafted prior and called it "deep image prior" (DIP) [15]. They
demonstrated that the natural structure of the deep CNN promotes the natural
signal and avoids the noise. While training, it inherently poses a higher
impedence to the values of network parameters which will output the noise than
the natural images. They demonstrated that network will represent noise only
if it is extensively trained to do so. In 2020, Wang et al demonstrated that
the deep image prior can also be used for reconstructing the phase images
[16]. They demonstrated that a deep CNN similarly poses a higher impedence to
the values of network parameters which will output the artifacts associated
with the missing phase information.

End-to-end learning frameworks have the burden of learning both the physics
i.e. the light transport and measurement model as well as the prior.
Alternatively, Physics informed (or physics aware) learning frameworks
integrate the partially or fully available physical or mathematical models
with the data-driven regression methods [17, 18]. The additional knowledge
from the former models can act as a regularizer to constrain the space of
admissible solutions to a smaller and manageable size [19]. The motivation is
to improve the generalization [20], reduce the big data requirement or
training time and achieve a better accuracy [18, 19]. In 1992, Psichogios et
al, proposed a hybrid neural network where a neural network estimates the
unknown parameters of interest, which are fed to a model predefined using the
known physical considerations [21]. Raissi et al, 2019 similarly applied a
neural network and their mathematical model in series, specifically in this
order to constraint the final output in a smaller admissible solution space
[19]. Chan et al, 2021 used this principle in coherent diffractive imaging.
They used their physics-aware optimization procedure as a refinement step on
the deep network’s output. They trained their network on simulated data,
obtained using a physics-informed data preparation pipeline to serve as a
representative of experimental data [22].

Figure 1: Recurrent approach of physics informed deep learning for phase
retrieval. Each iteration enforces the known physical laws and the deep prior
learnt from the data.

Unlike the previous physics informed machine learning methods, this letter
combines the knowledge of physical model and deep learning in a recurrent
manner. The letter uses a deep neural network originally trained for the
removal of additive white Gaussian noise (AWGN) to retrieve the phase
information from a diffracted intensity recording, without any fine-tuning of
network parameters. This transfer of knowledge has been achieved here by
enforcing the additional knowledge of diffraction physics in a recurrent
manner. The letter first presents a method labelled Phy-ZSN (short for physics
informed zero-shot network), which combines the classical error-reduction
algorithm and deep learning for obtaining more accurate phase estimations. The
method outperforms the phase retrieval with total-variation (TV) prior in
terms of image quality but at the cost of larger reconstruction time. To
alleviate this problem, this letter presents a second method labelled PhyTV-
ZSN, which combines both the handcrafted prior TV and the deep neural network
to obtain the improved reconstruction accuracy in a significantly less time.
Unlike other deep learning based phase reconstruction methods, the proposed
method does not need any labelled phase image data (ground truth).

The forward diffraction model can be represented as:

$I=\wp(|g|^{2}),\ \textnormal{where}\ g=f\otimes h(z,\lambda)$ (1)

$I$ is the recorded intensity image (diffraction pattern), $\wp$ is an
operator accounting for the noise of unknown statistics, $g$ is the complex
wave-field at the sensor plane, $f$ is the complex wave-field at the object
plane and $h$ is the depth and wavelength dependent impulse response function.
$z$ denotes the distance from the object to the sensor plane and $\lambda$
denotes the wavelength of the light. $|.|$ denotes the amplitude part of the
complex quantity (.) and $\otimes$ denotes the convolution operation. Under
the assumption of spatial and temporal coherence, impulse response function
$h(z,\lambda)$ is calculated using the scalar-diffraction equations [23]. The
problem statement is to estimate the complex object $f$ from the measurement
$I$. Additionally, the single scattering (first Born-approximation [24]) and
single depth approximations have been followed in this letter.

Zero-shot learning (ZSL) belongs to the family of transfer learning [25] where
a model performs the given task on samples from classes which were never
observed during the training phase. In ZSL, to facilitate this broader
generalization, some additional information about these unseen classes is
provided in the form of auxiliary information [26]. For example, a neural
network trained to identify horses can be used to identify zebras if the
auxiliary information is provided that zebras look like striped horses [25].
The proposed method uses a deep CNN trained for the additive white Gaussian
noise (AWGN) removal [27, 28], and transfers this learning to the phase
retrieval problem. The proposed framework operates in the recurrent manner
shown in the figure 1. The Physics fidelity block enforces the known
diffraction physics and measurement model (the auxiliary information in ZSL),
followed by passing the estimate through the pre-trained deep CNN for
enforcing the learnt object-prior. The process is repeated until the
convergence curve approaches a stationary region. The two blocks will be
discussed in details below.

Figure 2: Physics fidelity block based on the error-reduction algorithm.
Figure 3: PhyTV-ZSN: The object-prior block shows the cascaded TV and ZSN
priors. PhyTV-ZSN is identical to the Phy-ZSN except that the two TV units are
absent in the latter. Figure 4:  (a) The Poisson noise corrupted diffraction
pattern of the object "Cell" at $z$ = 1 mm. The reconstructed phase images
with (b) Angular spectrum wave propagation method. (c) Gerchberg-Saxton method
(amplitude constriant). (d) Error-reduction method with TV prior applied to
both amplitude and phase. (e) Phy-ZSN (proposed method 1). (f) PhyTV-ZSN
(proposed method 2). (g) Ground truth. Phase values were in the range 0 to
2$\pi$ but have been scaled to the range 0 to 255 for visualization. Figure 5:
Convergence pattern for the reconstructions: (a) Cell at $z$ = 1 mm,
reconstructed with Phy-ZSN and (b) with PhyTV-ZSN. (c) Lena at $z$ = 0.43 mm,
reconstructed with Phy-ZSN and (d) with PhyTV-ZSN. For PhyTV-ZSN (b and d),
the initial estimate is the reconstruction with TV prior.

In physics fidelity block, the error-reduction algorithm has been adapted. The
algorithm has been shown equivalent to the steepest-descent gradient search
method in the literature [1]. Rigorous mathematical analysis on the
convergence behaviour of this algorithm can be found in the reference [29].
The steps have been illustrated in figure 2. The object estimate is projected
to the sensor plane using the impulse response function $h(z,\lambda)$ via the
point-wise multiplication in the Fourier domain, amplitude is replaced by the
square-root of the measurement and finally backprojected to the object plane
via the point-wise division in the Fourier domain. Unlike most of impulse
response functions encountered in imaging, the $h(z,\lambda)$ in the frequency
domain $H(z,\lambda)$ is a phase only function. And hence the backprojection
using division operation in Fourier domain is possible (i.e. well-conditioned,
noise is not amplified by this division). $H(z,\lambda)$ is commonly referred
as the optical transfer function (OTF). The magnitude of $H(z,\lambda)$ is one
in the supported bandwidth and zero outside. Hence, the sampling rate must be
capped by this limit.

For free space (refractive index, $n=1$), the optical transfer function is
obtained by the following equation [30]:

$H(\textbf{v})=exp\Big{(}jk_{0}z\sqrt{1-(\lambda v_{x})^{2}-(\lambda
v_{y})^{2}}\Big{)};\ v^{2}_{x}+v^{2}_{y}<\frac{1}{\lambda^{2}}$ (2)
$H(\textbf{v})=0;\ \textnormal{otherwise}$ (3)

where, $\textbf{v}=(v_{x},v_{y})$ is the frequency coordinate vector,
$k_{0}=\frac{2\pi}{\lambda}$ is the wave-number.

Table 1: Quantitative assessment of reconstructions.

| Cell | Lena
---|---|---
Method | SSIM | PSNR | SSIM | PSNR
Angular spectrum Prop. | 0.6692 | 25.19 | 0.4740 | 13.59
Gerchberg-Saxton method | 0.4438 | 20.21 | 0.2990 | 15.00
Error-reduction with TV | 0.9161 | 29.24 | 0.8630 | 25.41
Phy-ZSN | 0.9629 | 34.45 | 0.8875 | 23.16
PhyTV-ZSN | 0.9601 | 34.84 | 0.9145 | 27.32

In learning based object-prior block in figure 1, a deep CNN based on the
U-Net and the ResNet architecture, referred as the DRUNet [27, 28] has been
used. The network was trained for denoising AWGN of randomly chosen variances
in the range [0,50]. The detailed steps for training this neural network can
be found in the reference [27]. This letter uses the terminology zero-shot
learning/network (ZSL/ZSN) because: (1) This trained network has not seen any
image from the class considered in this letter i.e. the estimates in a phase
retrieval process (figure 4b for instance is a backprojected diffraction
pattern without any object prior). (2) All the simulated diffraction patterns
used in this letter were corrupted with Poisson noise (not with AWGN, which
was used during the training).

The amplitude and phase values obtained from the physics fidelity block are
fed separately to this deep CNN. The inputs must be scaled in the range
between 0 to 255 because the network was trained with the examples of
intensity values in this specific range. The output must be scaled back to the
original range, by the applying the inverse of previous scaling.

As compared to the handcrafted prior TV, the deep neural network has a larger
computational complexity. To compare the reconstruction time, 10 iterations of
error reduction algorithm with TV prior take $\sim$1.3 seconds while the same
number of iterations of Phy-ZSN take $\sim$3.5 minutes. The testing was
performed on a system with 2.20 GHz Intel(R) Xeon(R) CPU and 13 GB RAM. To
alleviate this problem and for additional hand-crafted prior based
regularization control, another method which combines the two priors (TV and
deep network) in a cascaded manner has been proposed in this letter. This
method has been labelled as PhyTV-ZSN and has been shown in figure 3.

For faster reconstruction than Phy-ZSN, the following two-step implementation
with PhyTV-ZSN can be used:

Step 1: Disable the deep CNN part in PhyTV-ZSN until the iterative method
converges to a stationary point. This is a solution with TV prior.

Step 2: Enable the CNN part in the same. The iterative method will converge to
another stationary point. This is a solution with TV and deep prior.

10 iterations of PhyTV-ZSN take $\sim$3.67 minutes.

Figure 4 shows the phase reconstructions with different phase retrieval
methods, for comparison. The object is a HeLa cell phase image captured by an
Olympus IX83 microscope using a 20x/0.4 Ph2 Objective [31] and will be
referred as the "Cell" in the rest of this letter. The diffraction pattern was
simulated by the forward diffraction model described earlier, with the
following parameter values: pixel-pitch = 1.12 $\mu m$, light wavelength = 670
nm and object to sensor distance, $z$ = 1 mm. The simulated diffraction
pattern was corrupted by the Poisson noise. For figure 4d, TV prior
implementation was obtained using Chambolle’s method [6] and the TV prior was
applied on both the amplitude and phase separately. As the figure 4 and table
1 show, the proposed methods Phy-ZSN and PhyTV-ZSN give improved
reconstruction accuracies (TV is considered one of the best hand-crafted
object prior in the image processing community). The SSIM (structural
similarity index) and PSNR (peak signal-to-noise ratio) values are very low
for Gerchberg-Saxton method because of the non-suppression of the noise in the
reconstructed images. The convergence patterns have been shown in figure 5. In
figure 5 and table 1, the results corresponding to an another simulated
diffraction pattern from an object "Lena" at $z$ = 0.43 mm have also been
included. For Lena, when PhyTV-ZSN method was used, the total reconstruction
time reduced by a factor of 8.5, as compared to the Phy-ZSN. The convergence
rates are different for Lena and Cell because of the two different $z$ values,
i.e. because of the distance dependence of the impulse response function. The
Fienup’s hybrid input-output algorithm can also be used in the Physics
fidelity block (instead of error-reduction algorithm), former has been
reported to have better convergence speed and immunity to the stagnation
problem [1, 29].

In conclusion, this letter demonstrated that a deep neural network trained for
the image restoration task can also be used for solving the ill-posed phase
retrieval problem (from a single diffracted intensity image). This transfer of
knowledge was achieved here because the laws of diffraction physics were
enforced on the output of the deep neural network. The method did not require
any ground truth phase imaging data in any stage. The reconstruction quality
was better as compared to the reconstruction obtained with the total-variation
prior. Additionally, it was concluded that the joint use of handcrafted
priors, such as the TV prior and deep neural networks can facilitate a
faithful estimation in much lesser time, for a recurrent physics enforced deep
learning framework.

## Disclosures

The author declares no conflicts of interest.

## References

  * [1] J. R. Fienup, Applied optics 21, 2758 (1982).
  * [2] Y. Rivenson, Y. Zhang, H. Günaydın, D. Teng, and A. Ozcan, Light: Science & Applications 7, 17141 (2018).
  * [3] R. W. Gerchberg, Optik 35, 237 (1972).
  * [4] T. Latychevskaia and H.-W. Fink, Physical review letters 98, 233901 (2007).
  * [5] O. Mudanyali, D. Tseng, C. Oh, S. O. Isikman, I. Sencan, W. Bishara, C. Oztoprak, S. Seo, B. Khademhosseini, and A. Ozcan, Lab on a Chip 10, 1417 (2010).
  * [6] A. Chambolle, Journal of Mathematical imaging and vision 20, 89 (2004).
  * [7] J. M. Bioucas-Dias and M. A. Figueiredo, IEEE Transactions on Image processing 16, 2992 (2007).
  * [8] V. Katkovnik and J. Astola, JOSA A 29, 44 (2012).
  * [9] F. Jolivet, F. Momey, L. Denis, L. Méès, N. Faure, N. Grosjean, F. Pinston, J.-L. Marié, and C. Fournier, Optics express 26, 8923 (2018).
  * [10] Y. Rivenson, Y. Wu, and A. Ozcan, Light: Science & Applications 8, 1 (2019).
  * [11] Y. Wu, Y. Rivenson, Y. Zhang, Z. Wei, H. Günaydin, X. Lin, and A. Ozcan, Optica 5, 704 (2018).
  * [12] D. Misell, Journal of Physics D: Applied Physics 6, 2200 (1973).
  * [13] D. Paganin, A. Barty, P. McMahon, and K. A. Nugent, Journal of microscopy 214, 51 (2004).
  * [14] K. Wang, J. Dou, Q. Kemao, J. Di, and J. Zhao, Optics letters 44, 4765 (2019).
  * [15] D. Ulyanov, A. Vedaldi, and V. Lempitsky, “Deep image prior,” in _Proceedings of the IEEE conference on computer vision and pattern recognition,_ (2018), pp. 9446–9454.
  * [16] F. Wang, Y. Bian, H. Wang, M. Lyu, G. Pedrini, W. Osten, G. Barbastathis, and G. Situ, Light: Science & Applications 9, 1 (2020).
  * [17] G. Barbastathis, A. Ozcan, and G. Situ, Optica 6, 921 (2019).
  * [18] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, Nature Reviews Physics pp. 1–19 (2021).
  * [19] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Journal of Computational Physics 378, 686 (2019).
  * [20] D. Heaven, Nature 574, 163 (2019).
  * [21] D. C. Psichogios and L. H. Ungar, AIChE Journal 38, 1499 (1992).
  * [22] H. Chan, Y. S. Nashed, S. Kandel, S. O. Hruszkewycz, S. K. Sankaranarayanan, R. J. Harder, and M. J. Cherukara, Applied Physics Reviews 8, 021407 (2021).
  * [23] J. W. Goodman, _Introduction to Fourier optics_ (Roberts and Company Publishers, 2005).
  * [24] M. Born and E. Wolf, _Principles of optics: electromagnetic theory of propagation, interference and diffraction of light_ (Elsevier, 2013).
  * [25] W. Wang, V. W. Zheng, H. Yu, and C. Miao, ACM Transactions on Intelligent Systems and Technology (TIST) 10, 1 (2019).
  * [26] Y. Xian, C. H. Lampert, B. Schiele, and Z. Akata, IEEE transactions on pattern analysis and machine intelligence 41, 2251 (2018).
  * [27] K. Zhang, Y. Li, W. Zuo, L. Zhang, L. Van Gool, and R. Timofte, arXiv preprint arXiv:2008.13751 (2020).
  * [28] K. Zhang, W. Zuo, S. Gu, and L. Zhang, “Learning deep cnn denoiser prior for image restoration,” in _Proceedings of the IEEE conference on computer vision and pattern recognition,_ (2017), pp. 3929–3938.
  * [29] H. H. Bauschke, P. L. Combettes, and D. R. Luke, JOSA A 19, 1334 (2002).
  * [30] A. Ozcan and E. McLeod, Annual review of biomedical engineering 18, 77 (2016).
  * [31] O. Burri and R. Guiet, http://doi.org/10.5281/zenodo.3232478 (2019).

PhyZSN_OL