# Remote Sensing Image Translation via
Style-Based Recalibration Module and Improved Style Discriminator

Tiange Zhang, Feng Gao, Junyu Dong, Qian Du This work was supported in part by
the National Key Research and Development Program of China under Grant
2018AAA0100602, in part by the National Natural Science Foundation of China
under Grant U1706218, and in part by the Key Research and Development Program
of Shandong Province under Grant 2019GHY112048. _(Corresponding auther: Junyu
Dong and Feng Gao)_ Tiange Zhang, Feng Gao, and Junyu Dong are with the School
of Information Science and Engineering, Ocean University of China, Qingdao
266100, China (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>dongjunyu@ouc.edu.cn;). Qian Du is with the Department of Electrical and
Computer Engineering, Mississippi State University, Starkville, MS 39762 USA
(e-mail: du@ece.msstate.edu).

###### Abstract

Existing remote sensing change detection methods are heavily affected by
seasonal variation. Since vegetation colors are different between winter and
summer, such variations are inclined to be falsely detected as changes. In
this letter, we proposed an image translation method to solve the problem. A
style-based recalibration module is introduced to capture seasonal features
effectively. Then, a new style discriminator is designed to improve the
translation performance. The discriminator can not only produce a decision for
the fake or real sample, but also return a style vector according to the
channel-wise correlations. Extensive experiments are conducted on season-
varying dataset. The experimental results show that the proposed method can
effectively perform image translation, thereby consistently improving the
season-varying image change detection performance. Our codes and data are
available at `https://github.com/summitgao/RSIT_SRM_ISD`

###### Index Terms:

Change detection, remote sensing, image-to-image translation, GAN.

##  I Introduction

Nowadays, with the rapid advancement of the Earth observation program, an
ever-growing number of optical remote sensing images are available. These
images are widely applied in land cover change detection [1], environmental
analysis [2], disaster monitoring [3], etc. Among these applications, change
detection is one of the most important techniques to monitor the changed
surfaces on the Earth. In most cases, researches are mainly interested in
changes caused by human activities and natural disasters. However, such
changes may be difficult to be detected due to some interfering factors. For
instance, vegetation colors are different between summer and winter. These
changes caused by seasonal variation may yield false alarms.

Many change detection methods are affected by seasonal variation. As
illustrated in Fig. 1, the change map generated by the PCAKM [1] contains many
noisy regions. These regions are mainly vegetations whose appearances are
different between winter and summer. Such variations are falsely detected as
changes. Therefore, if a function can be learnt to map remote sensing images
captured in summer to images captured in winter, the false alarms in change
detection caused by seasonal variation can be alleviated effectively.

Figure 1: Illustration of the problem in change detection caused by season
variation. The change map generated by PCAKM [1] is polluted with noisy
regions. Specifically, the vegetation colors are different between summer and
winter. The change detection method considers them as significant changes.

The image-to-image translation [4] task aims to transform one representation
of a scene into another. So image translation is natural to be adopted in a
change detection framework to solve the aforementioned problem caused by
seasonal variation. Many researchers have leveraged generative adversarial
networks [5] into image translation. Isola et al. [4] used a conditional
generative adversarial network (cGAN) for supervised pairwise image
translation. Zhu et al. [6] proposed an unsupervised CycleGAN by imposing a
cycle consistency loss, which is capable of learning the distinctive semantic
difference between image sets from two domains and translating them
correspondingly without direct pair-wise supervision.

Several attempts have been made to translate remote sensing images between
different domains. Enomoto et al. [7] proposed a GAN-based image translation
method for environmental monitoring. Niu et al. [8] proposed an image
translation method between optical and SAR data based on conditional GAN for
land cover classification. Saha et al. [9] employed the CycleGAN to learn the
transcoding between SAR and optical images, where the CycleGAN provides useful
semantic features for building change detection. The above-mentioned methods
are designed for the image translation task between SAR and optical images.
However, image translation focusing on seasonal variation is seldom explored.

Although the tasks of SAR-to-optical translation and seasonal translation look
similar, the latter is more challenging, since it focuses on transferring the
seasonal features while maintaining other more important features. However,
such seasonal features are difficult to be captured during the translation
procedure. To tackle this issue, we consider the seasonal features of a remote
sensing image as its style features, and introduce a style-based recalibration
module to better capture the style features. Besides, a new style
discriminator is designed to improve the translation performance during the
training stage.

In this letter, the main contributions of our work can be summarized in
twofold.

  * •

We present an image translation method for season-varying multitemporal remote
sensing images, which can effectively improve the change detection
performance. To the best of our knowledge, this is the first attempt to
incorporate image translation in the season-varying image change detection
task.

  * •

To effectively capture the seasonal features, a style-based recalibration
module and a new style discriminator are designed to improve the translation
performance during the training phase.

  * •

Extensive experiments on the season-varying dataset demonstrate the
superiority of the proposed method over the state-of-the-art methods. As a
byproduct, we released the codes and data to facilitate other researchers.

##  II Methodology

The architecture of the proposed method is shown in Fig. 2, where CycleGAN [6]
is used as the backbone since it can be trained in an unsupervised way by
utilizing two collections of unpaired data without labels. The model has two
generators: $G_{XY}$ represents the mapping from domain $X$ (summer) to domain
$Y$ (winter), while $G_{YX}$ denotes the mapping from domain $Y$ to domain
$X$. Given the image collections $\left\\{x_{i}\right\\}_{i=1}^{N_{s}}$ where
$x_{i}\in X$ and $\left\\{y_{j}\right\\}_{j=1}^{N_{w}}$ where $y_{j}\in Y$.
$N_{s}$ denotes the number of images in the summer domain, and $N_{w}$ denotes
the number of images in the winter domain. The discriminator $D_{X}$ is
trained to distinguish between images $\left\\{x\right\\}$ from real domain
$X$ and translated images $\left\\{G_{YX}(y)\right\\}$. The discriminator
$D_{Y}$ is trained to distinguish between images $\\{y\\}$ from real domain
$Y$ and translated images $\left\\{G_{XY}(x)\right\\}$.

Figure 2: The framework of our proposed network. Domain $X$ represents images
captured in summer while domain $Y$ represents images captured in winter.
Figure 3: The architecture of the basic components from our image-to-image
translation network. (a) represents the Generator, which takes an image as the
input and outputs the translated image. (b) represents the new Style
Discriminator, which is designed to figure out whether the input is a real
image or a translated one. It outputs a scalar of $[0,1]$ as the decision as
well as a style vector which encode the style of the input to compute the
style loss between real and translated images. (c) illustrates the
SRMConvBlock utilized in both Generator and Style Discriminator.

###  II-A The Generator with Style-Based Recalibration Module

As shown in Fig. 3(a), the generator follows an Encoder-Transformation-Decoder
structure. The encoder contains three convolutional layers to encode the image
to the feature space, while the decoder utilizes two fractionally-strided
convolutions and one convolutional layer to reconstruct the image from the
feature space. The essential transformation module is comprised of 9 residual
blocks which manipulate the feature maps to conduct the translation process.

In order to represent style features effectively in feature space, we
introduce the style-based recalibration module (SRM)[10] into the residual
block. The SRM is comprised of two parts: _style pooling_ and _style
integration_. The style pooling operator extracts style features from each
channel by summarizing feature responses across spatial dimensions. It is
followed by the style integration operator, which produces example-specific
style weights by utilizing style features via channel-wise operation. Finally,
the style weights recalibrate the feature maps to either emphasize or suppress
their information. Note that SRMConvBlock represents the residual block which
integrates SRM. Details can be found in Fig. 3 (c).

In _style pooling_ , the average mean and standard deviation of the feature
map are selected as style features. Denote the input feature map as
$F\in\mathbb{R}^{C\times H\times W}$, where $C$ is the total number of
channels, and $H$, $W$ denote spatial dimensions. The style feature
$T\in\mathbb{R}^{C\times 2}$ can be computed in each channel by:

$\displaystyle\mu_{c}=\frac{1}{HW}\sum\limits_{i=1}^{H},\sum\limits_{j=1}^{W}f_{cij}$
(1)

$\displaystyle\sigma_{c}=\sqrt{\frac{1}{HW}\sum\limits_{i=1}^{H}\sum\limits_{j=1}^{W}(f_{cij}-\mu_{c})^{2}},$
(2)

$\displaystyle t_{c}=[\mu_{c}~{},~{}\sigma_{c}],$ (3)

where the style vector $t_{c}$ serves as a summary descriptor of the style
information for $c^{th}$ channel. To be specific, $\mu_{c}$ and $\sigma_{c}$
are concatenated to form $t_{c}$.

In _style integration_ , the style features are converted into channel-wise
style weights. Ref. [10] applied a learnable channel-wise fully connected
layer to achieve style integration. However, it is found that the fully
connected layer has high computational burden. Therefore, in our
implementation, a 1-d convolutional layer with fewer parameters but a
competitive capability is used to integrate the style. Given
$T\in\mathbb{R}^{C\times 2}$ from style pooling as the input, the style
integration weight $G\in\mathbb{R}^{C\times 1}$ performs as the following in
each channel:

$g_{c}=\sigma(w*t_{c}),$ (4)

where $g_{c}$ denotes the weight for $c^{th}$ channel, $\sigma(\cdot)$ denotes
the sigmoid function, $w$ denotes the convolution kernel, and * denotes
convolution. Finally, the input feature $F\in\mathbb{R}^{C\times H\times W}$
is recalibrated by the channel-wise style integration weights
$G\in\mathbb{R}^{C\times 1}$ to produce the output
$\hat{F}\in\mathbb{R}^{C\times H\times W}$ as:

$\hat{F}=G\cdot F.$ (5)

Here, channel-wise multiplication is employed between the style integration
weight $G$ and the input feature map $F$.

###  II-B Style Discriminator

In order to perform better image translation, we propose a new style
discriminator, which not only produces a decision for the fake or real sample,
but also returns a style vector according to the channel-wise correlations. By
constraining the L1 loss between two style vectors of the translated image and
the real target domain image, the style discriminator can distinguish the
style information better.

The proposed style discriminator is illustrated in Fig. 3 (b). Given an input
image, it is encoded to the feature space by adopting several convolutional
layers and SRM layers. Then, the encoded feature $M$ is fed into two different
modules.

On the one hand, after convolution, $M$ is mapped into a set of overlapping
image patches. Then, the decision results of each patch will be generated, and
the results are averaged to a final scalar of [0, 1]. On the other hand, we
consider the correlations of different channels as another kind of style
representation for evaluating the style distance between translated images and
target images. This correlation style representation can explore more style
information across channels because in _style pooling_ , the style features
are computed in every separate channel. First, $M\in\mathbb{R}^{C\times
H\times W}$ is provided to max and average pooling layers. After each max
pooling and average pooling, the aggregated feature maps are combined and fed
into the next two pooling layers. Through four pooling and fusion operations,
the feature map $M$ is transformed to $V\in\mathbb{R}^{C\times 1\times 1}$.
Since each channel from total number of $C$ is now represent by a single
value, $V\in\mathbb{R}^{C\times 1\times 1}$ is reshaped to
$V\in\mathbb{R}^{1\times C}$. Therefore, the correlation between different
channels can be described as a matrix $V^{T}V\in\mathbb{R}^{C\times C}$ by
element-wise multiplication between the transpose of $V$ and itself. Only
keeping the upper triangular elements of $V^{T}V$ while setting the other
elements as 0, a style vector $v$ with the length of $C\times C$ is generated
after flattening the correlation matrix. Then the L1 loss can be applied to
measure the distance between two style vectors that represent translated and
target images.

###  II-C Training Overflow

As mentioned before, our goal is to learn the mapping functions $G_{XY}:X\to
Y$ and $G_{YX}:Y\to X$. The objective function contains three types of terms:
1) _adversarial loss_ $\mathcal{L}_{\textrm{GAN}}$ aims to match the
distribution of translated images to the data distribution in the target
domain; 2) _style loss_ $\mathcal{L}_{S}$ calculates the differences between
the translated image and the target image with two style vectors which are
generated by the style discriminator $D_{X}$ and $D_{Y}$; 3) _cycle
consistency loss_ $\mathcal{L}_{cyc}$ aims to prevent the mappings $G_{XY}$
and $G_{YX}$ from contradicting each other.

The overall objective function can be summarized as:

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle\mathcal{L}_{\textrm{GAN}}(G_{XY},D_{Y},X,Y)+\mathcal{L}_{S}(G_{XY}(X),Y)$
(6) $\displaystyle+$
$\displaystyle\mathcal{L}_{\textrm{GAN}}(G_{YX},D_{X},Y,X)+\mathcal{L}_{S}(G_{YX}(Y),X)$
$\displaystyle+$ $\displaystyle\lambda\mathcal{L}_{cyc}(G_{XY},G_{YX})$

where $\lambda$ controls the relative importance of the cycle consistency
loss, and $\lambda=10$ in the following experiments. For both
$\mathcal{L}_{\textrm{GAN}}$, the negative log-likelihood objective is
replaced with a least-squares loss, which is more stable during training and
can generate higher quality samples.

Specifically, the generator $G_{XY}$ aims to minimize the following objective
function:

$\displaystyle\mathcal{L}_{G_{XY}}$ $\displaystyle=$
$\displaystyle\mathcal{L}_{\textrm{GAN}}+\mathcal{L}_{cyc}+\mathcal{L}_{id},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(7) $\displaystyle\mathcal{L}_{\textrm{GAN}}$ $\displaystyle=$
$\displaystyle\mathbb{E}_{x\sim
p_{\textrm{data}}(x)}[(D_{Y}(G_{XY}(x))-1)^{2}],~{}~{}~{}~{}$ (8)
$\displaystyle\mathcal{L}_{cyc}$ $\displaystyle=$
$\displaystyle\mathbb{E}_{x\sim
p_{\textrm{data}}(x)}[\|G_{YX}(G_{XY}(x))-x\|_{1}],$ (9)
$\displaystyle\mathcal{L}_{id}$ $\displaystyle=$
$\displaystyle\mathbb{E}_{y\sim
p_{\textrm{data}}(y)}[\|G_{XY}(y)-y\|_{1}],~{}~{}~{}~{}~{}$ (10)

where an additional loss $\mathcal{L}_{id}$ is used to regularize the
generator to approximate an identity mapping when real samples of the target
domain are provided as the input to the generator. The generator $G_{YX}$
operates in the same way.

For the style discriminator $D_{Y}$, it aims to minimize:

$\displaystyle\mathcal{L}_{D_{Y}}=\mathcal{L}_{\textrm{GAN}}+\mathcal{L}_{S},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(11) $\displaystyle\mathcal{L}_{\textrm{GAN}}=\mathbb{E}_{y\sim
p_{\textrm{data}}(y)}[(D_{Y}(y)-1)^{2}]~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle+\mathbb{E}_{x\sim
p_{\textrm{data}}(x)}[(D_{Y}(G_{XY}(x)))^{2}],$ (12)
$\displaystyle\mathcal{L}_{S}=\mathbb{E}_{x\sim p_{\textrm{data}}(x),y\sim
p_{\textrm{data}}(y)}[\|(D_{Y}(G_{XY}(x))-D_{Y}(y)\|_{1}],$ (13)

where $\mathcal{L}_{\textrm{GAN}}$ minimizes the output decision whereas
$\mathcal{L}_{S}$ regularizes the style vectors. The style discriminator
$D_{X}$ is trained in the same manner.

##  III Experimental Results and Discussions

###  III-A Datasets and Training Details

Two sets of high resolution remote sensing images are collected for
experiments from Google Earth. They are acquired at the same area from three
cities in China, including Beijing, Tianjin, and Qingdao. Images in summer are
captured in December 2016, and images in winter are captured in June 2019. All
images are randomly cropped into 256$\times$256 pixels. We generate 869 images
as the training set, and 217 images as the testing set. The proposed method is
trained from scratch with a learning rate of 0.0002. We keep the same learning
rate for the first 100 epochs and linearly decay the rate to zero over the
next 100 epochs.

###  III-B Image-to-Image Translation Evaluation

In order to evaluate our translation model, we employ three criteria:
Inception Score (IS) [11], Fréchet Inception Distance (FID) [12], and Kernel
Inception Distance (KID) [13]. IS measures the quality and diversity of the
translated images. FID compares features of real and generated images that
extracted by a layer from a pre-trained Inception Network. KID computes the
squared maximum mean discrepancy between the Inception feature representations
of real target domain images and translated images.

TABLE I: Quantitative Evaluation on Translation Method | Winter to summer
---|---
IS $\uparrow$ |  FID $\downarrow$ |  KID $\downarrow$
MUNIT[14] | 3.61 | 210.09 | 7.72
U-GAT-IT[15] | 3.81 | 146.58 | 2.51
CycleGAN[6] | 3.71 | 130.76 | 2.41
SECycleGAN[16] | 3.77 | 125.08 | 1.49
Ours | 3.72 | 119.75 | 1.20
Method | Summer to winter
IS $\uparrow$ | FID $\downarrow$ | KID $\downarrow$
MUNIT[14] | 3.38 | 216.41 | 8.48
U-GAT-IT[15] | 3.88 | 183.26 | 5.12
CycleGAN[6] | 3.99 | 129.02 | 1.04
SECycleGAN[16] | 3.90 | 125.89 | 0.89
Ours | 4.18 | 124.19 | 0.82

Four state-of-the-art unsupervised image-to-image translation methods
including MUNIT [14], U-GAT-IT [15], CycleGAN [6], and SECycleGAN [16] are
selected as the baselines to evaluate the proposed method. The quantitative
results are presented in Table I. It shows that our proposed image-to-image
translation model achieves the lowest FID and KID scores, and it is quite
competitive on IS score with the CycleGAN, SECycleGAN and U-GAT-IT. These
quantitative evaluation results indicate that our image translation model has
learnt both mapping functions $G_{XY}:X\to Y$ and $G_{YX}:Y\to X$ well after
training.

Besides quantitative evaluation, we also conduct two computational efficiency
analysis on our proposed method, and the results are illustrated in Table II
and Table III. First, Table II presents the translation time of the proposed
method tested on images of 256$\times$256 and 512$\times$512 pixels. It
demonstrates that our method is capable of performing fast image translation.
Second, Table III shows the amount of parameters and GFLOPs of our image
translation model compared with others. It demonstrates that the introduced
style-based recalibration module and the new style discriminator can improve
the translation performance only with less computational burden added.

TABLE II: Computational Efficiency Analysis I Mappings | Translation time per image (s)
---|---
256$\times$256 |  512$\times$512
Translation (Winter to Summer) | 0.039 | 0.188
Translation (Summer to Winter) | 0.040 | 0.185
TABLE III: Computational Efficiency Analysis II Method |  Params |  GFLOPs
---|---|---
MUNIT[14] | 15.026 M | 77.320
U-GAT-IT[15] | 10.587 M | 52.506
CycleGAN[6] | 11.378 M | 56.832
SECycleGAN[16] | 11.378 M | 56.832
Ours | 12.634 M | 56.843

###  III-C Performance Analysis for Season-Varying Image Change Detection

We now consider the use of the proposed method to assist in season-varying
image change detection. In our approach, we first apply the method to
translate multitemporal images to the same domain. Then, two well-known
algorithms, PCAKM[1] and GETNET[17], are used to obtain the change detection
results. To verify the change detection performance, false alarms (FA), missed
alarms (MA), overall errors (OE), and percentage correct classification (PCC)
are adopted as the evaluation criteria. One typical region are selected for
discussion to verify the effectiveness of the proposed image translation
method.

Figure 4: Visualized change detection results of the selected region. TABLE IV: Change Detection Results of the selected region Methods | FA | MA | OE | PCC(%)
---|---|---|---|---
PCAKM | 21518 | 131 | 21649 | 66.96
GETNET | 17663 | 144 | 17807 | 72.83
(Winter to Summer) + PCAKM | 390 | 146 | 536 | 99.18
(Summer to Winter) + PCAKM | 2682 | 26 | 2708 | 95.87
(Winter to Summer) + GETNET | 149 | 119 | 268 | 99.59
(Summer to Winter) + GETNET | 1903 | 53 | 1956 | 97.01

As can be observed from Fig. 4, new structures are built in this region. At
first, we generate change maps directly by using PCAKM [1] and GETNET [17]. We
can see that if PCAKM or GETNET is directly employed to generate the change
map, many false alarms are generated since the vegetation colors are different
between summer and winter. Then, Fig. 4 also demonstrates that the results
generated through image translation are less noisy than the original methods.
It is evident that the false alarms are greatly reduced, which means that the
proposed image translation method improves the change detection performance to
some extent. The corresponding quantitative results are illustrated in Table
IV. By employing the seasonal translation, the PCC values are improved
significantly. Moreover, the change maps generated by winter-to-summer
translation achieve 99.18% and 99.59% for both methods, rather close to the
ground truth. All these compared quantitative results reveal that the proposed
scheme can produce much better change detection results generally, not only on
one change detection algorithm.

##  IV Conclusion

In this letter, we present a remote sensing image translation method, which
can capture seasonal features to improve the detection performance on real
changes of interest. A style-based recalibration module and a new style
discriminator are designed to improve the translation during the training
phase. The encouraging experimental results verify the effectiveness of the
proposed method and indicate the necessity of integrating the proposed
translation scheme in the season-varying image change detection task. Although
our method achieves good results in change detection, we believe that there is
still room for improvement, such as integrating the image translation and
change detection operators into an end-to-end network.

## References

  * [1] T. Celik, “Unsupervised change detection in satellite images using principal component analysis and k-means clustering,” _IEEE Geosci. Remote Sens. Lett._ , vol. 6, no. 4, pp. 772–776, Apr. 2009.
  * [2] L. Bruzzone and D. F. Prieto, “Automatic analysis of the difference image for unsupervised change detection,” _IEEE Trans. Geosci. Remote Sens._ , vol. 38, no. 3, pp. 1171–1182, Mar. 2000.
  * [3] F. Kahraman, M. Imamoglu, and H. F. Ates, “Disaster damage assessment of buildings using adaptive self-similarity descriptor,” _IEEE Geosci. Remote Sens. Lett._ , vol. 13, no. 8, pp. 1188–1192, 2016\.
  * [4] P. Isola, J.-Y. Zhu, T. Zhou, and A. A. Efros, “Image-to-image translation with conditional adversarial networks,” _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , July 2017.
  * [5] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, et al. “Generative adversarial nets,” _Proceedings of the International Conference on Neural Information Processing Systems (NeurIPS)_ , 2014, pp. 2672–2680.
  * [6] J.-Y. Zhu, T. Park, P. Isola, and A. A. Efros, “Unpaired image-to-image translation using cycle-consistent adversarial networks,” _Proceedings of the IEEE International Conference on Computer Vision (ICCV)_ , Oct. 2017.
  * [7] K. Enomoto, K. Sakurada, W. Wang, et al. “Image translation between SAR and optical imagery with generative adversarial nets,” _Proceedings of the IEEE International Geoscience and Remote Sensing Symposium (IGARSS)_ , 2018, pp. 1752–1755.
  * [8] X. Niu, D. Yang, H. Yang, and Y. Dou, “Image translation between high-resolution remote sensing optical and SAR data using conditional GAN,” _Proceedings on the Pacific Rim Conference on Multimedia (PCM)_ , 2018, pp. 245–255.
  * [9] S. Saha, F. Bovolo, and L. Bruzzone, “Building change detection in VHR SAR images via unsupervised deep transcoding,” _IEEE Trans. Geosci. Remote Sens._ , pp. 1–13, 2020.
  * [10] H. Lee, H.-E. Kim, and H. Nam, “SRM: A style-based recalibration module for convolutional neural networks,” _Proceedings of the IEEE International Conference on Computer Vision (ICCV)_ , Oct. 2019.
  * [11] T. Salimans, I. Goodfellow, W. Zaremba, V. Cheung, A. Radford, X. Chen, and X. Chen, “Improved techniques for training GANs,” _Proceedings of the International Conference on Neural Information Processing Systems (NeurIPS)_ , 2016.
  * [12] M. Heusel, H. Ramsauer, T. Unterthiner, et al. “GANs trained by a two time-scale update rule converge to a local nash equilibrium,” _Proceedings of the International Conference on Neural Information Processing Systems (NeurIPS)_ , 2017.
  * [13] M. Binkowski, D. J. Sutherland, M. Arbel, and A. Gretton, “Demystifying MMD GANs,” _Proceedings of the International Conference on Learning Representations (ICLR)_ , 2018.
  * [14] X. Huang, M.-Y. Liu, S. Belongie, and J. Kautz, “Multimodal Unsupervised Image-to-Image Translation,” _Proceedings of the European Conference on Computer Vision (ECCV)_ , 2018.
  * [15] J. Kim, M. Kim, H. Kang, and K. Lee, “U-GAT-IT: Unsupervised Generative Attentional Networks with Adaptive Layer-Instance Normalization for Image-to-Image Translation,” _Proceedings of the International Conference on Learning Representations (ICLR)_ , 2020.
  * [16] Q. Wang, J.-Y. Gao, W. Lin, and Y. Yuan, “Pixel-wise crowd understanding via synthetic data,” _International Journal of Computer Vision_ , vol. 129, no.1, pp. 225–245, Jan. 2021.
  * [17] Q. Wang, Z.-H. Yuan, Q. Du and X.-L. Li, “GETNET: A general end-to-end two-dimensional CNN framework for hyperspectral image change detection,” _IEEE Trans. Geosci. Remote Sens._ , vol. 57, no.1, pp. 3–13, Jan. 2019.


# Evaluating the Convergence Limit of Quantum Neural Tangent Kernel

Trong Duong  Korea Advanced Institute of Science and Technology (KAIST)

(Date: July 2022)

###### Abstract.

Quantum variational algorithms have been one of major applications of quantum
computing with current quantum devices. There are recent attempts to establish
the foundation for these algorithms. A possible approach is to characterize
the training dynamics with quantum neural tangent kernel. In this work, we
construct the kernel for two models, Quantun Ensemble and Quantum Neural
Network, and show the convergence of these models in the limit of infinitely
many qubits. We also show applications of the kernel limit in regression
tasks.

This work was conducted during my internship at the Centre for Quantum
Technologies (CQT), Singapore, under the valuable guidance of Patrick
Rebentrost.

##  1\. Introduction

The study of overparametrized models in machine learning has drawn attentions
thanks to the increasing interest in deep learning. These models, which have
many more parameters than training samples, have empirically good learning
capability and performance on the test data. This unexpected generalization
may be explained by a training algorithm such as gradient descent somehow
inducing implicit regularization of the model. In particular, prior work
showed that in some circumstances, even complex nonlinear models might behave
as if in the kernel-based training with regularization on the corresponding
Reproducing Kernel Hilbert Space (RKHS) [8]. In contrast, other work showed
that under different conditions some models can exhibit different non-Hilbert
norm implicit regularization [7, 6, 15, 19].

On the other hand, quantum variational models is a common class of algorithms
in the noisy intermediate-scale quantum (NISQ) era. Essentially, the model
uses the natural quantum objective function
$f(x,\theta))=\bra{\psi}H\ket{\psi}$ for some Hermitian observable $H$, where
$\ket{\psi}=U(x,\theta)\ket{0}$ is prepared by applying unitary identified by
the input and model’s parameters on the initial state. The unitary matrix
contain information about parameters and input data. There are prior work
devoted to establish the notion and convergence of quantum neural tangent
kernel [18, 12], prove the fair constant-ness of the tangent kernel at
initialization and during training [11, 1], and the effect of noisy
measurements on the kernel [10].

We have so far a certain number of quantum variational models that could
entail a converging tangent kernel. The limit of the kernel at the limit of
infinite dimension is particularly useful. For example, there are efficient
analytical form of the tangent kernel for classical fully connected neural
networks and convolutional neural networks at the infinite limit [8, 2], which
serve as an effective benchmark for models of the same kind. From my humble
experience, there has not yet an analytical form for quantum neural tangent
kernel. This is because quantum neural networks are usually realized by some
specific architecture the convergence limit of the kernel depends on. Hence to
utilize the benefits of the kernel in high dimensional space, one still have
to carry out operations of large unitary matrices. The computation quickly
becomes inefficient as the dimension of the Hilbert space scales exponentially
to the number of qubits.

We attempt to provide a quantum model where one can efficiently compute the
analytical neural tangent kernel at large dimension without perform large
matrix operations.

##  2\. Preliminaries

###  2.1. Neural Tangent Kernel

It has been shown that training a neural network with gradient descent is
highly related to “neural tangent kernel”. In particular, when the neural
network can be approximated by a linearized model, one can describe the
dynamics of the model by a simple first-order differential equation. The
situation happens in the so-called “kernel regime”, which usually occurs in an
overparametrized model in which every parameter stays almost unchanged [8]. We
first specify the relations between training a neural network and the neural
tangent kernel and conditions for the kernel regime to arise.

#### Neural Tangent Kernel

Let the loss function to be optimize is
$L(\theta)=\frac{1}{2}\sum_{i}|f(\mathbf{x}_{i};\theta)-y_{i}|^{2}=\frac{1}{2}\|\mathbf{f}(\mathbf{X},\theta)-\mathbf{y}\|^{2}$.
Parameters are updated at every step according to gradient descent
$\delta\theta=\theta(t+1)-\theta(t)=-\eta\nabla L(\theta(t))$. With
infinitesimal learning rate $\eta$, the update rule is equivalent to the
gradient flow

(1) $\displaystyle\dot{\theta}=\frac{\partial\theta}{\partial t}$
$\displaystyle=\nabla L(\theta(t))$ (2)
$\displaystyle=-\nabla\mathbf{f}(\theta)(\mathbf{f}(\theta)-\mathbf{y})$

where $\mathbf{f}(\theta)\equiv f(\mathbf{X};\theta)$. We can induce the
dynamics of the model in the function space

(3) $\displaystyle\dot{\mathbf{f}}(\theta)$
$\displaystyle=\nabla\mathbf{f}(\theta)^{T}\dot{\theta}$ (4)
$\displaystyle=-\nabla\mathbf{f}(\theta)^{T}\nabla\mathbf{f}(\theta)(\mathbf{f}(\theta)-\mathbf{y})$
(5) $\displaystyle=-\mathbf{K}(\theta)(\mathbf{f}(\theta)-\mathbf{y})$

The quantity
$\mathbf{K}(\theta)=\nabla\mathbf{f}(\theta)^{T}\nabla\mathbf{f}(\theta)$ is
called the “neural tangent kernel” (NTK). In some special cases when this
kernel stay almost constant during training
$\mathbf{K}(\theta(t))\approx\mathbf{K}(\theta_{0})$ the training dynamics can
be described by linear first-order equation
$\dot{\mathbf{f}}(\theta)=-\mathbf{K}(\theta_{0})(\mathbf{f}(\theta)-\mathbf{y})$,
which turns out to have a simple solution

(6)
$\displaystyle\mathbf{f}(\theta(t))=\mathbf{y}+e^{-\mathbf{K}(\theta_{0})t}(\mathbf{f}(\theta_{0}))-\mathbf{y})$

The situation when the kernel is constant is referred to as the “kernel
regime” or “lazy regime”, usually characterized by constant gradient
$\nabla\mathbf{f}(\theta(t))\approx\nabla\mathbf{f}(\theta_{0})\Rightarrow\mathbf{K}(\theta(t))\approx\mathbf{K}(\theta_{0})$.
This happens when the parameters do change much during training, i.e. they
stay close to the initialization. Expand the output function near the
initialization

(7) $\displaystyle f(\mathbf{x};\theta)=f(\mathbf{x};\theta_{0})+\nabla
f(\mathbf{x};\theta_{0})^{T}(\theta-\theta_{0})$

Notice that the output function $f(\mathbf{x};\theta)$ is linear in the
parameters but non-linear on the input data with data feature
$\Phi(\mathbf{x})=\nabla f(\mathbf{x};\theta_{0})^{T}$.

It was proved that the linearization that assumes fairly constant gradient
during training occurs when every hidden layer of the neural network has
infinite width and parameters initialized with normal distributions [8], and
when the output function is scaled by some big factor $f\mapsto\alpha f$ [3].
Moreover, for infinitely wide neural networks, the kernel is fairly constant
over different initializations $\theta_{0}$ and there is a recursive formula
to calculate the kernel $\mathbf{K}_{\infty}(\theta_{0})$ the kernel converges
to at the infinite limit. The proof of convergence for classical neural
network heavily relies on the fact that the neural network behaves like a
Gaussian process at the infinite-width limit. In this project we attempt to
find such a convergence for the neural tangent kernel of a quantum neural
network (QNN).

###  2.2. Reproducing kernel Hilbert space

To understand the importance of kernel to (quantum) machine learning, we have
to study Reproducing Kernel Hilbert Space (RKHS), the central concept of
kernel theory. The type of kernel that is particularly useful is positive
definite and symmetric (PDS) kernel, which is almost always implied in the
literature about kernel methods. We provide a brief introduction to the kernel
theory. Readers can find a comprehensive discussion about the kernel method in
[13].

######  Definition 1.

A kernel is a real-valued bivariate function
$k:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ defined over the input
space $\mathcal{X}$. The kernel is positive-definite and symmetric (PDS) if

  1. (1)

for all $m\in\mathbb{N},c_{i}\in\mathbb{R},x_{i}\in\mathcal{X},i\in[m]$,

(8) $\sum_{i,j=1}^{m}c_{i}c_{j}k(x_{i},x_{j})\geq 0.$

  2. (2)

$k(x,x^{\prime})=k(x^{\prime},x)$ for all $x,x^{\prime}\in\mathcal{X}$.

For example, the tangent kernel, or any kernel defined by an inner product, is
PDS because $\sum_{i,j}c_{i}c_{j}\nabla f(x_{i})^{T}\nabla
f(x_{j})=\|\sum_{i}c_{i}\nabla f(x_{i})\|^{2}$ and the symmetricity is
straightforward. A key property of PDS kernels is that there exist a Hilbert
space, i.e, a complete vector space equipped with inner product, called
Reproducing Kernel Hilbert Space. This RKHS should be distinguished from the
quantum Hilbert space $\mathcal{H}$, hence we denote it by $\mathcal{G}$.

######  Theorem 1 (Reproducing Kernel Hilbert Space).

Given a positive-definite and symmetric kernel function $k$ over the input
space $\mathcal{X}$, there exists a Hilbert space $\mathcal{G}$ and a map
$\phi:\mathcal{X}\rightarrow\mathcal{G}$ such that

  1. (1)

$\mathcal{G}=\operatorname{span}\\{g_{x}:x\in\mathcal{X}\\}$ and shall include
the limit of every Cauchy sequence, where
$g_{x}:\mathcal{X}\rightarrow\mathbb{R}$ is a functional given by
$g_{x}(\cdot)=k(x,\cdot)$.

  2. (2)

$\phi(x)\in\mathcal{G}$ defined by $\phi(x)=g_{x}=k(x,\cdot)$

  3. (3)

$k(x,x^{\prime})\equiv\expectationvalue{k(x,\cdot),k(x^{\prime},\cdot)}_{\mathcal{G}}=\expectationvalue{\phi(x),\phi(x^{\prime})}_{\mathcal{G}}$
for all $x,x^{\prime}\in\mathcal{X}$.

  4. (4)

(Reproducing property) For all $h\in\mathcal{G}$,

(9)
$h(x)=\expectationvalue{h,\phi(x)}_{\mathcal{G}}=\expectationvalue{h,k(x,\cdot)}_{\mathcal{G}}$

The reproducing property comes from the construction of $\mathcal{G}$ that we
can write $h=\sum_{i=1}a_{i}k(x_{i},\cdot),x_{i}\in\mathcal{X}$. Then
$\phi(x)$ “reproduces” the values $h(x)$ for every $h\in\mathcal{G}$,

(10)
$h(x)=\sum_{i}a_{i}k(x_{i},x)=\expectationvalue{\sum_{i}a_{i}k(x_{i},\cdot),k(x,\cdot)}_{\mathcal{G}}=\expectationvalue{h,\phi(x)}_{\mathcal{G}}.$

The application of kernel and RKHS in machine learning comes from the famous
Representer Theorem, which establishes the fact that we only need to search
through a linear combination of finite kernel functionals for a supervised
learning problem.

######  Theorem 2 (Representer Theorem).

Given a kernel $k$, its induced RKHS $\mathcal{G}$ and feature map $\phi$.
Define an objective function $F:\mathcal{G}\rightarrow\mathcal{R}$ as a
regularized empirical loss function wrt. some dataset
$\\{(x_{i},y_{i}),i\in[m]\\}$

(11) $F(h)=G(\|h\|_{\mathcal{G}})+L((h(x_{1}),y_{1}),\dots,(h(x_{m}),y_{m}))$

such that $G(\cdot)$ is non-decreasing. If
$\operatorname*{arg\,min}_{h\in\mathcal{G}}F(h)$ has a solution, then there
exists a solution of the form
$h^{\prime}=\sum_{i=1}^{m}a_{i}k(x_{i},\cdot)=\sum_{i=1}^{m}a_{i}\phi(x_{i})$.
Furthermore, if $G$ is strictly increasing, all solutions must have the given
form.

The problem of minimizing $F(h)$ with mean square loss function, where $h$
admits the form $\sum_{i=1}^{m}a_{i}\phi(x_{i})$ can be solved by ordinary
least square with the kernel trick that essentially utilizes the simple
evaluation of $k(x,x^{\prime})$ in many cases. Moreover, it was shown that
variational quantum models, a huge class of quantum machine learning models,
are equivalent to functions in the corresponding RKHS [16]. Hence finding an
optimal model can be solved by finding an optimal solution on the RKHS.

##  3\. Quantum Ensemble Model

We consider a model that, given a collection
$\\{(a_{n},U_{n},W_{n}):a_{n}\in\mathbb{R},U_{n},W_{n}\in\mathcal{U}(d)\\},n\in[N]$,
where $\mathcal{U}(d)$ is the unitary group of dimension $d$. For a qubit
system of $n$ qubits, the corresponding quantum Hilbert space has the
dimension of $d=2^{n}$. Define a model

(12)
$f(x)=\frac{\sqrt{d}}{\sqrt{N}}\sum_{n=1}^{N}a_{n}\bra{0}U_{n}^{\dagger}S^{\dagger}(x)W_{n}^{\dagger}HW_{n}S(x)U_{n}\ket{0}.$

This type of model, which looks like an ensemble of natural quantum objective
functions, is non-conventional. However it would shed light on the concept of
tangent kernel in the context of quantum machine learning. With a training
dataset $\\{(x_{p},y_{p}):y_{p}\in\mathbb{R}\\},p\in[P]$, the model can be
optimized with respect to $a_{n}$, which is an ordinary least square problem.

(13) $\displaystyle\min_{a}\frac{1}{2}\sum_{p=1}^{P}|f(x_{p})-y_{p}|^{2}$ (14)
$\displaystyle=$ $\displaystyle\min_{a}\frac{1}{2}||Fa-y||^{2},$

where F is a matrix whose elements are given by

(15)
$F_{pn}=\frac{\sqrt{d}}{\sqrt{N}}\bra{0}U_{n}^{\dagger}S^{\dagger}(x_{p})W_{n}^{\dagger}HW_{n}S(x_{p})U_{n}\ket{0},$

i.e. every element is the expectation of the observable $H$ evaluated at the
state $\ket{\psi_{pn}}=W_{n}S(x_{p})U_{n}\ket{0}$. The solution to this
problem is

(16) $\hat{a}=(FF^{T})^{-1}Fy$

For the model to have good performance, one has to has a sufficiently large
$N$, and the unitaries $U_{n},W_{n}$ should be diverse. A general guide to
choose $N$ and the unitaries seems to be lacking up to my knowledge.

On the other hand, we can observe that the gradient $\nabla_{a}f(x_{p})$ is
also the $p$-th row of $F$. From the definition of tangent kernel, the kernel
matrix is

(17) $\displaystyle K(x_{p},x_{p^{\prime}})$
$\displaystyle=\nabla_{a}f(x_{p})^{T}\nabla_{x}f(x_{p^{\prime}})$
$\displaystyle=\frac{d}{N}\sum_{n=1}^{N}\bra{\psi_{pn}}H\ket{\psi_{pn}}\bra{\psi_{p^{\prime}n}}H\ket{\psi_{p^{\prime}n}}.$

Although this kernel matrix is independent to $a$, it depends on the choice of
unitaries. Suppose we have no prior knowledge about selected unitaries and
assume them to be sampled independently from $\mathcal{U}(d)$ with the Haar
measure. In that case one can derive the expectation value of a polynomial of
entries of random unitaries with respect to the unitary Haar measure [14, 4].
Here we use a symbolic integration program using tensor network [5] to
simplify the computation.

$\displaystyle\mathbb{E}K(x_{p},x_{p^{\prime}})=$ $\displaystyle
d\times\mathbb{E}_{\mathcal{U},\mathcal{W}}\left[\bra{0}U^{\dagger}S^{\dagger}(x_{p})W^{\dagger}HWS(x_{p})U\ket{0}\bra{0}U^{\dagger}S^{\dagger}(x_{p^{\prime}})W^{\dagger}HWS(x_{p^{\prime}})U\ket{0}\right]$
$\displaystyle=$ $\displaystyle
d\times\frac{d^{2}\Tr(H)^{2}+d\Tr(H)^{2}+d\Tr(H^{2})+\Tr(H^{2})s_{pp^{\prime}}}{(d^{2}-1)^{2}}$
$\displaystyle+$ $\displaystyle
d\times\frac{d^{2}\Tr(H^{2})+d\Tr(H)^{2}+d\Tr(H^{2})+\Tr(H)^{2}s_{pp^{\prime}}}{d^{2}(d^{2}-1)^{2}}$
$\displaystyle-$ $\displaystyle
d\times\frac{d^{2}\Tr(H)^{2}+d^{2}\Tr(H^{2})+2d\Tr(H)^{2}+2d\Tr(H^{2})+(\Tr(H)^{2}+\Tr(H^{2}))s_{pp^{\prime}}}{d(d^{2}-1)^{2}}$
(18) $\displaystyle=$
$\displaystyle\frac{d^{3}+d^{2}-d-s_{pp^{\prime}}}{d(d^{3}+d^{2}-d-1)}\Tr(H)^{2}+\frac{s_{pp^{\prime}}-1}{d^{3}+d^{2}-d-1}\Tr(H^{2})$

where $s_{pp^{\prime}}=|\Tr(S(x_{p^{\prime}})S^{\dagger}(x_{p}))|^{2}$. A
class of common observable operators in practice is Pauli measurement, i.e.
$H\in\\{I,X,Y,Z\\}^{n}$, all of which satistfy $\Tr(H)=0$ and $\Tr(H^{2})=d$.
In the limit $d\rightarrow\infty$, the limit of
$\mathbb{E}K(x_{p},x_{p^{\prime}})$ aligns with the intuition that kernel is a
similarity measure.

(19) $\displaystyle\mathbb{E}K(x_{p},x_{p^{\prime}})$
$\displaystyle=\frac{d(s_{pp^{\prime}}-1)}{d^{3}+d^{2}-d-1}$ (20)
$\displaystyle\rightarrow\left\\{\begin{matrix}1,&\text{ if
}s_{pp^{\prime}}\approx d^{2}\\\ 0,&\text{ o.w.}\end{matrix}\right.$

The condition for the positive limit to hold is $S(x_{p})\approx
S(x_{p^{\prime}})$. For example, when the encoding block is
$S(x)=\left(e^{-ixZ/2}\right)^{\otimes n}$, the kernel looks like a partial
Fourier sum as shown in Appendix B

(21) $\displaystyle
s_{pp^{\prime}}=\binom{2n}{n}+2\sum_{k=1}^{n}\binom{2n}{n-k}\cos(k(x_{p}-x_{p^{\prime}}))$

In most of the cases, there are no such efficient analytical forms to compute
$s_{pp^{\prime}}$. Fortunately, one can estimate it with the Hadamard test.
Considering a qubit system with a probe qubit initialized at
$\outerproduct{0}{0}$ and $n$ qubits at the completely mixed state
$\frac{1}{d}I_{d},d=2^{n}$. The mixed state can be prepared by applying the
Pauli-$X$ gate on each qubit uniformly at random. One can prepare the system
into the state

(22)
$\displaystyle\rho_{n+1}=\frac{1}{2d}\left(\outerproduct{0}{0}I_{d}+\outerproduct{0}{1}U^{\dagger}+\outerproduct{1}{0}U+\outerproduct{1}{1}I_{d}\right)=\frac{1}{2d}\begin{pmatrix}I_{d}&U^{\dagger}\\\
U&I_{d}\end{pmatrix}$

by applying a Hadamard gate on the probe qubit, followed by
$U=S(x_{p^{\prime}})S^{\dagger}(x_{p})$ in the ancillary system controlled on
the probe qubit. Measuring the probe qubit gives sufficient information to
compute $\Tr(U)$ since $\expectationvalue{X}=\frac{1}{d}\real\Tr(U)$ and
$\expectationvalue{Y}=-\frac{1}{d}\imaginary\Tr(U)$. A simple Chernoff bound
can show that one needs $O(1/\varepsilon^{2})$ samples to estimate either the
real or imaginary part within the absolute error $\varepsilon$, independent to
the size of the system. Hence if restrained to use only multiple parallel
local encoding blocks for the encoding unitary
$S(x)=\left(e^{-ix\sigma}\right)^{\otimes n}$, we can evaluate the kernel
matrix efficiently regardless of the dimension. This provides a quick
alternative technique for regression and classification problems.

(a) Standard techniques (b) Kernel method

Figure 1. Optimization methods for quantum ensemble model (12) with $n=4$
qubits and $N=2000$ terms. The function to learn is
$f(x)=\cos(x)+3\sin(2x)-2\cos(3x)$ defined on $[0,2\pi]$.

Although the kernel in the ensemble model might have the form of a partial
Fourier sum (21), the corresponding kernel function fails be universal in the
sense that it can approximate any function in the RKHS within arbitrary
precision. This is because we have no control over the coefficients of the
Fourier sum even though we can have as many Fourier frequencies as desired.

##  4\. Quantum Neural Network Model

Define a quantum neural network to have the form $WS(x)U\ket{0}$, where $U$
and $W$ might be realized by sequences of parametrized and fixed quantum
gates. Let $Y(x)=\bra{0}U^{\dagger}S(x)^{\dagger}W^{\dagger}HWS(x)U\ket{0}$ be
the natural objective of the QNN with a Hermitian observable $H$. One can
assume the observable and the encoding block are diagonal, i.e.
$H=\text{diag}(h_{1},\dots,h_{d}),h_{j}\in\mathbb{R}$ and
$S(x)=\text{diag}(e^{i\lambda_{1}},\dots,e^{i\lambda_{d}})$. The function
value can be written as a partial Fourier sum

(23) $\displaystyle Y$
$\displaystyle=\sum_{ijk=1}^{d}(U^{\dagger})_{1i}(S^{\dagger})_{ii}(W^{\dagger})_{ij}H_{jj}W_{jk}S_{kk}U_{k1}$
$\displaystyle=\sum_{ijk}h_{j}e^{i(\lambda_{k}-\lambda_{i})}\bar{u}_{i1}u_{k1}\bar{w}_{ji}w_{jk}$
$\displaystyle=\sum_{ik}e^{i(\lambda_{k}-\lambda_{i})}\bar{u}_{i1}u_{k1}\sum_{j}h_{j}\bar{w}_{ji}w_{jk}$
$\displaystyle=\sum_{i}|u_{i1}|^{2}\sum_{j}h_{j}|w_{ji}|^{2}+\sum_{i\neq
k}e^{i(\lambda_{k}-\lambda_{i})}\bar{u}_{i1}u_{k1}\sum_{j}h_{j}\bar{w}_{ji}w_{jk}$

Denote (1) $\alpha_{i}=|u_{i1}|^{2}$, (2) $a_{i}=\sum_{j}h_{j}|w_{ji}|^{2}$,
(3) $\beta_{ik}=\bar{u}_{i1}u_{k1}$, (4)
$b_{ik}=h_{j}\bar{w}_{ji}w_{jk},i,k\in[d],i\neq k$. We show in Appendix C that

(24) $\displaystyle\alpha_{i}$ $\displaystyle\sim\operatorname{Beta}(1,d-1)$
$\displaystyle a_{i}$
$\displaystyle\sim\mathcal{N}\left(1,\frac{(d-1)}{(d+1)}\Tr(H^{2})\right)$
$\displaystyle\real\beta_{ik},\imaginary\beta_{ik}$
$\displaystyle\sim\operatorname{Laplace}(0,1/2d)$ $\displaystyle\real
b_{ik},\imaginary b_{ik}$
$\displaystyle\sim\mathcal{N}\left(0,\frac{1}{2d^{2}}\Tr(H^{2})\right)$

We continue to write

(25) $\displaystyle Y$ $\displaystyle=\sum_{i}\alpha_{i}a_{i}+\sum_{i\neq
k}e^{i(\lambda_{k}-\lambda_{i})}\beta_{ik}b_{ik}$
$\displaystyle=\sum_{i}\alpha_{i}a_{i}+2\sum_{i<k}\real\left\\{e^{i(\lambda_{k}-\lambda_{i})\beta_{ik}b_{ik}}\right\\}$
$\displaystyle=\sum_{i}\alpha_{i}a_{i}+2\sum_{i<k}\cos(\lambda_{k}-\lambda_{i})(\real\beta_{ik}\real
b_{ik}-\imaginary\beta_{ik}\imaginary b_{ik})$
$\displaystyle\quad\quad\quad\quad\quad\quad-\sin(\lambda_{k}-\lambda_{i})(\real\beta_{ik}\imaginary
b_{ik}+\imaginary\beta_{ik}\real b_{ik})$

We again define real variables
$\gamma_{ik}=\real\beta_{ik},\delta_{ik}=\imaginary\beta_{ik}$ and
$c_{ik}=\real b_{ik},d_{ik}=\real b_{ik}$. Then

(26)
$Y=\sum_{i}\alpha_{i}a_{i}+2\sum_{i<k}\cos(\lambda_{k}-\lambda_{i})(\gamma_{ik}c_{ik}-\delta_{ik}d_{ik})-\sin(\lambda_{k}-\lambda_{i})(\gamma_{ik}d_{ik}+\delta_{ik}c_{ik})$

Since $U$ and $W$ each has $d^{2}$ real degree of freedom, we have the full
control of $U$ via $d^{2}$ real variables
$\alpha_{i},\gamma_{ik},\delta_{ik}$. This suggests us to use
$\alpha_{i},\gamma_{ik},\delta_{ik}$ as independent parameters of $U$. It is
similar for $W$ and $b_{i},c_{ik},d_{ik}$. The gradient of the function with
respect to those variables is given by

(27) $\displaystyle\partial_{\alpha_{i}}Y$ $\displaystyle=a_{i}$
$\displaystyle\partial_{a_{i}}Y$ $\displaystyle=\alpha_{i}$
$\displaystyle\partial_{\gamma_{ik}}Y$
$\displaystyle=2\cos(\lambda_{k}-\lambda_{i})c_{ik}-2\sin(\lambda_{k}-\lambda_{i})d_{ik}$
$\displaystyle\partial_{\delta_{ik}}Y$
$\displaystyle=-2\cos(\lambda_{k}-\lambda_{i})d_{ik}-2\sin(\lambda_{k}-\lambda_{i})c_{ik}$
$\displaystyle\partial_{c_{ik}}Y$
$\displaystyle=2\cos(\lambda_{k}-\lambda_{i})\gamma_{ik}-2\sin(\lambda_{k}-\lambda_{i})\delta_{ik}$
$\displaystyle\partial_{d_{ik}}Y$
$\displaystyle=-2\cos(\lambda_{k}-\lambda_{i})\delta_{ik}-2\sin(\lambda_{k}-\lambda_{i})\gamma_{ik}$

The tangent kernel function $K(x,x^{\prime})=\nabla Y^{T}(x)\nabla
Y(x^{\prime})$ is equal to

(28) $\displaystyle K(x,x^{\prime})$
$\displaystyle=\sum_{i}a_{i}^{2}+\sum_{i}\alpha_{i}^{2}$
$\displaystyle+8\sum_{i<k}\cos(\lambda_{k}-\lambda_{i})\cos(\lambda_{k}^{\prime}-\lambda_{i}^{\prime})(c^{2}_{ik}+\gamma^{2}_{ik})$
$\displaystyle+8\sum_{i<k}\sin(\lambda_{k}-\lambda_{i})\sin(\lambda_{k}^{\prime}-\lambda_{i}^{\prime})(d^{2}_{ik}+\delta^{2}_{ik})$

This form of kernel function might be shown to be universal [17] given as many
qubits as needed. It is simple to verify the following expected value:
$\mathbb{E}[\alpha_{i}^{2}]=\frac{2}{d(d+1)}$,
$\mathbb{E}[a_{i}^{2}]=1+\frac{d-1}{d+1}\Tr(H^{2})$,
$\mathbb{E}[\gamma_{i}^{2}]=\mathbb{E}[\delta_{i}^{2}]=\frac{1}{2d^{2}}$, and
$\mathbb{E}[c_{i}^{2}]=\mathbb{E}[d_{i}^{2}]=\frac{1}{2d^{2}}\Tr(H^{2})$.

The analytical limit of the tangent kernel at the infinite limit is

(29) $\displaystyle\mathbb{E}[K(x,x^{\prime})]$
$\displaystyle=\frac{2}{d+1}+\left(d+\frac{d(d-1)}{d+1}\Tr(H^{2})\right)$
$\displaystyle+\frac{4}{d^{2}}(1+\Tr(H^{2}))\sum_{i<k}\cos(\lambda_{k}-\lambda_{i})\cos(\lambda_{k}^{\prime}-\lambda_{i}^{\prime})$
$\displaystyle+\frac{4}{d^{2}}(1+\Tr(H^{2}))\sum_{i<k}\sin(\lambda_{k}-\lambda_{i})\sin(\lambda_{k}^{\prime}-\lambda_{i}^{\prime})$
$\displaystyle=d+\frac{2+d(d-1)\Tr(H^{2})}{d+1}$
$\displaystyle+\frac{4}{d^{2}}(1+\Tr(H^{2}))\sum_{i<k}\cos(\lambda_{k}-\lambda_{i}-\lambda_{k}^{\prime}+\lambda_{i}^{\prime}).$

One might proceed to use this expected kernel for regression with kernel
method. For robustness, it is recommended to normalize $Y$ by
$\frac{1}{\sqrt{d}}$ to avoid $\Theta(d)$ terms in the expression of
$\mathbb{E}[K(x,x^{\prime})]$.

##  5\. Conclusion

In this manuscript we have show the convergence of quantum neural tangent
kernel at large dimension limit. For both Quantum Ensemble and Quantum Neural
Network models, the limit of the kernel function has a concise analytical form
from which efficient evaluation or estimation could be carried out. We expect
this could aid the finding of new quantum models as an effective benchmark for
the performance of supervised learning problems.

## References

  * [1] Erfan Abedi, Salman Beigi and Leila Taghavi  “Quantum lazy training”  In _Quantum_ 7 Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften, 2023, pp. 989
  * [2] Sanjeev Arora et al.  “On exact computation with an infinitely wide neural net”  In _Advances in Neural Information Processing Systems_ 32, 2019
  * [3] Lenaic Chizat, Edouard Oyallon and Francis Bach  “On lazy training in differentiable programming”  In _Advances in Neural Information Processing Systems_ 32, 2019
  * [4] Benoit Collins and Piotr Śniady  “Integration with respect to the Haar measure on unitary, orthogonal and symplectic group”  In _Communications in Mathematical Physics_ 264.3 Springer, 2006, pp. 773–795
  * [5] Motohisa Fukuda, Robert König and Ion Nechita  “RTNI—A symbolic integrator for Haar-random tensor networks”  In _Journal of Physics A: Mathematical and Theoretical_ 52.42 IOP Publishing, 2019, pp. 425303
  * [6] Suriya Gunasekar, Jason D Lee, Daniel Soudry and Nati Srebro  “Implicit bias of gradient descent on linear convolutional networks”  In _Advances in Neural Information Processing Systems_ 31, 2018
  * [7] Suriya Gunasekar et al.  “Implicit regularization in matrix factorization”  In _Advances in Neural Information Processing Systems_ 30, 2017
  * [8] Arthur Jacot, Franck Gabriel and Clément Hongler  “Neural tangent kernel: Convergence and generalization in neural networks”  In _Advances in neural information processing systems_ 31, 2018
  * [9] Olav Kallenberg and Olav Kallenberg  “Foundations of modern probability”  Springer, 1997
  * [10] Junyu Liu, Zexi Lin and Liang Jiang  “Laziness, barren plateau, and noise in machine learning”  In _arXiv preprint arXiv:2206.09313_ , 2022
  * [11] Junyu Liu et al.  “Analytic theory for the dynamics of wide quantum neural networks”  In _Physical Review Letters_ 130.15 APS, 2023, pp. 150601
  * [12] Junyu Liu et al.  “Representation learning via quantum neural tangent kernels”  In _PRX Quantum_ 3.3 APS, 2022, pp. 030323
  * [13] Mehryar Mohri, Afshin Rostamizadeh and Ameet Talwalkar  “Foundations of machine learning”  MIT press, 2018
  * [14] Zbigniew Puchała and Jarosław Adam Miszczak  “Symbolic integration with respect to the Haar measure on the unitary group”  In _arXiv preprint arXiv:1109.4244_ , 2011
  * [15] Pedro Savarese, Itay Evron, Daniel Soudry and Nathan Srebro  “How do infinite width bounded norm networks look in function space?”  In _Conference on Learning Theory_ , 2019, pp. 2667–2690  PMLR
  * [16] Maria Schuld  “Supervised quantum machine learning models are kernel methods”  In _arXiv preprint arXiv:2101.11020_ , 2021
  * [17] Maria Schuld, Ryan Sweke and Johannes Jakob Meyer  “Effect of data encoding on the expressive power of variational quantum-machine-learning models”  In _Physical Review A_ 103.3 APS, 2021, pp. 032430
  * [18] Norihito Shirai, Kenji Kubo, Kosuke Mitarai and Keisuke Fujii  “Quantum tangent kernel”  In _arXiv:2111.02951_ , 2021
  * [19] Blake Woodworth et al.  “Kernel and rich regimes in overparametrized models”  In _Conference on Learning Theory_ , 2020, pp. 3635–3673  PMLR
  * [20] Karol Życzkowski and Hans-Jürgen Sommers  “Average fidelity between random quantum states”  In _Physical Review A_ 71.3 APS, 2005, pp. 032313

##  Appendix A Probability distribution of expected value of Hermitian
observables

We want to determine the distribution of
$Y=Y(\ket{\psi})=\bra{\psi}H\ket{\psi}$ when $\ket{\psi}$ is sampled with
respect to the Haar measure of quantum states, i.e. finding pdf $p(y)$.

To begin with, we parametrize the complex projective space $CP^{N-1}$ almost
everywhere with the coordinates
$\xi=[\xi_{1},\dots,\xi_{N-1}]^{T}\in\mathbb{C}^{N-1}$. Express a random
$N$-dimensional normalized vector corresponding to the point $\xi$ using local
coordinates

(30)
$\ket{\phi(\xi,\bar{\xi})}=\frac{[1,\xi_{1},\dots,\xi_{N-1}]^{T}}{\sqrt{1+\xi^{\dagger}\xi}}$

The Fubini Study volume element is given by

(31)
$d\mu_{CP^{N-1}}(\xi,\bar{\xi})=\frac{(N-1)!}{\pi^{N-1}}\frac{\prod_{k=1}^{N-1}d\xi_{k}d\bar{\xi}_{k}}{(1+\xi^{\dagger}\xi)^{N}}.$

In general, assume $H$ can be diagonalized into
$\text{diag}(h_{0},\dots,h_{N-1})$ in which the eigenvalues are in descending
order. The pdf is given by:

$\displaystyle p(y)$
$\displaystyle=\int_{CP^{N-1}}\delta\left(\frac{h_{0}+\sum_{k=1}^{N-1}h_{k}|\xi_{k}|^{2}}{1+\xi^{\dagger}\xi}-y\right)\frac{(N-1)!}{\pi^{N-1}}d\mu_{CP^{N-1}}(\xi,\bar{\xi})$
(32)
$\displaystyle=\int_{CP^{N-1}}\delta\left(\frac{h_{0}+\sum_{k=1}^{N-1}h_{k}|\xi_{k}|^{2}}{1+\xi^{\dagger}\xi}-y\right)\frac{(N-1)!}{\pi^{N-1}}\frac{\prod_{k=1}^{N-1}d\xi_{k}d\bar{\xi}_{k}}{(1+\xi^{\dagger}\xi)^{N}}$

Since $|\xi_{k}|^{2}=\|[\Re{\xi_{k}},\Im{\xi_{j}}]\|$, we can replace
$\xi_{k}$ by two real variables that are its real and imaginary parts. Apply a
change of variable

(33)
$\real\xi_{k}=r_{k}\cos(\theta_{k}),\imaginary\xi_{k}=r_{j}\sin(\theta_{k}),$

with

(34)
$\left|\frac{\partial(\real\xi_{k},\imaginary\xi_{k})}{\partial(r_{k},\theta_{k})}\right|=r_{k}$

followed by another change of variable $z_{k}=r^{2}_{k},k=1,\dots,N-1$. The
pdf becomes

$\displaystyle=\frac{(N-1)!}{2^{N-1}}\int_{r_{1},\dots,r_{N-1}=0}^{\infty}\delta\left(\frac{h_{0}+\sum_{k=1}^{N-1}h_{k}r^{2}_{k}}{1+\sum_{k=1}^{N-1}r_{k}^{2}}-y\right)\frac{r_{1}\dots
r_{N-1}}{(1+\sum_{k=1}^{N-1}r_{k}^{2})^{N}}\prod_{k=1}^{N-1}dr_{k}$ (35)
$\displaystyle=(N-1)!\int_{z_{1},\dots,z_{N-1}=0}^{\infty}\delta\left(\frac{h_{0}+\sum_{k=1}^{N-1}h_{k}z_{k}}{1+\sum_{k=1}^{N-1}z_{k}}-y\right)\frac{\prod_{k=1}^{N-1}dz_{k}}{(1+\sum_{k=1}^{N-1}z_{k})^{N}}$

We want to get rid of the Dirac Delta function. Consider

(36)
$g(z_{N-1})=\frac{h_{0}+\sum_{k=1}^{N-2}h_{k}z_{k}+h_{N-1}z_{N-1}}{1+\sum_{k=1}^{N-2}z_{k}+z_{N-1}}$

Its derivative is

(37)
$g^{\prime}(z_{N-1})=\frac{h_{N-1}-h_{0}+\sum_{k=1}^{N-2}(h_{N-1}-h_{k})z_{k}}{(1+\sum_{k=1}^{N-2}z_{k}+z_{N-1})^{2}}$

Due to the order of $h_{k}$, the derivative $g^{\prime}(z_{N-1})$ is negative,
i.e. $|g^{\prime}(z_{N-1})|=-g^{\prime}(z_{N-1})$. Also, the peak of the Dirac
Delta function occurs at

(38) $\displaystyle z^{*}_{N-1}$
$\displaystyle=\frac{y-h_{0}+\sum_{k=1}^{N-2}(y-h_{k})z_{k}}{h_{N-1}-y}$ (39)
$\displaystyle 1+\sum_{k=1}^{N-2}z_{k}+z^{*}_{N-1}$
$\displaystyle=\frac{h_{N-1}-h_{0}+\sum_{k=1}^{N-2}(h_{N-1}-y+1-h_{k}}{h_{N-1}-y}$
(40) $\displaystyle\Rightarrow g^{\prime}(z^{*}_{N-1})$
$\displaystyle=\frac{(h_{N-1}-y)^{2}(h_{N-1}-h_{0}+\sum_{k=1}^{N-2}z_{k})}{\left(h_{N-1}-h_{0}+\sum_{k=1}^{N-2}(h_{N-1}-y+1-h_{k})z_{k}\right)^{2}}$

Therefore we rewrite

(41)
$\displaystyle\delta\left(\frac{h_{0}+\sum_{k=1}^{N-1}h_{k}z_{k}}{1+\sum_{k=1}^{N-1}z_{k}}-y\right)$
$\displaystyle=\frac{\delta(z_{N-1}-\frac{y-h_{0}+\sum_{k=1}^{N-2}(y-h_{k})z_{k}}{h_{N-1}-y})}{|g^{\prime}(z^{*}_{N-1})|}$
$\displaystyle=\delta\left(z_{N-1}-\frac{y-h_{0}+\sum_{k=1}^{N-2}(y-h_{k})z_{k}}{h_{N-1}-y}\right)$
$\displaystyle\times\frac{\left(h_{N-1}-h_{0}+\sum_{k=1}^{N-2}(h_{N-1}-h_{k}+1-y)z_{k}\right)^{2}}{(h_{N-1}-y)^{2}\left(h_{N-1}-h_{0}+\sum_{k=1}^{N-2}z_{k}\right)}.$

This allows us to eliminate the Dirac Delta function upon integrating over
$z_{N-1}$. In addition since $g(0)>g(\infty)$, the density $p(y)$ is positive
only when $g(0)\geq y>g(\infty)$, i.e.

(42) $\displaystyle h_{N-1}<y$
$\displaystyle\leq\frac{h_{0}+\sum_{k=1}^{N-2}h_{k}z_{k}}{1+\sum_{k=1}^{N-2}z_{k}}$

The first constraint is trivial. We use the second constraint to derive bounds
for $z_{k}$. It is equivalent to

(43) $h_{0}-y+\sum_{k=1}^{n-2}(h_{k}-y)z_{k}\geq 0$

Note that $p(y)=0$ for $y\geq h_{0}$. Assume $h_{m-1}>y>h_{m}$ for some
integer $1\leq m\leq N-1$. Rewrite the constraint as

(44) $\displaystyle\sum_{k=1}^{N-2}(y-h_{k})z_{k}\leq
h_{0}-y+\sum_{k=1}^{m-1}(h_{k}-y)z_{k}$

Notice $h_{k}-y>0,k=0,\dots,m-1$ and $y-h_{k}>0,k=m,\dots,N-2$. Thus the
bounds for $z_{1},\dots,z_{N-2}$ are given as follows

from to $\displaystyle z_{N-2}\quad$ $\displaystyle 0$
$\displaystyle\frac{h_{0}-y+\sum_{k=1}^{m-1}(h_{k}-y)z_{k}+\sum_{k=m}^{N-3}(h_{k}-y)z_{k}}{y-h_{N-2}}$
$\displaystyle z_{N-3}\quad$ $\displaystyle 0$
$\displaystyle\frac{h_{0}-y+\sum_{k=1}^{m-1}(h_{k}-y)z_{k}+\sum_{k=m}^{N-4}(h_{k}-y)z_{k}}{y-h_{N-3}}$
$\displaystyle\dots$ $\displaystyle z_{m}\quad$ $\displaystyle 0$
$\displaystyle\frac{h_{0}-y+\sum_{k=1}^{m-1}(h_{k}-y)z_{k}}{y-h_{m}}$
$\displaystyle z_{m-1},\dots,z_{2},z_{1}\quad$ $\displaystyle 0\quad$
$\displaystyle\infty$

The general procedure is follows. For each of non-trivial gaps
$y\in(y_{m-1},y_{m})$, we derive the bounds of variables $z_{1},\dots,z_{N-1}$
as above, which can be used to compute, say by a computer-algebra system
(CAS),

(45)
$p(y)=\frac{(N-1)!}{(h_{N-1}-y)^{2}}\int\frac{\left(h_{N-1}-h_{0}+\sum_{k=1}^{N-2}(h_{N-1}-h_{k}+1-y)z_{k}\right)^{2}}{h_{N-1}-h_{0}+\sum_{k=1}^{N-2}z_{k}}dz_{N-2}\dots
dz_{1}$

##  Appendix B Analytical form for Pauli-$Z$ encoding

Assume $S(x)=(e^{-ixZ})^{\otimes
n}=\operatorname{diag}(e^{-ix},e^{ix})^{\otimes n}$. Let
$\vec{v}(k)\in\\{0,1\\}^{n}$ be the binary representation of position index
$k,k\in\\{0,\dots,2^{n}-1\\}$ along the diagonal matrix. The number of “1” and
“0” in $v$ are $\expectationvalue{\vec{1},\vec{v}}$ and
$\expectationvalue{\vec{0},\vec{v}}$ respectively, where $\vec{1}$ is the all-
one vector. One can verify that the expansion of $S(x)$ is a diagonal matrix
whose $k$-th element is

(46) $\displaystyle
e^{ix\expectationvalue{\vec{1},\vec{v}(k)}-\expectationvalue{\vec{1},\vec{1}-\vec{v}(k)}}$
$\displaystyle=e^{ix\expectationvalue{\vec{1},2\vec{v}(k)-\vec{1}}}$
$\displaystyle\equiv e^{ixq_{k}}$

We can write

(47) $\displaystyle\Tr(S_{2}S_{1}^{\dagger})\Tr(S_{2}^{\dagger}S_{1})$
$\displaystyle=\sum_{j,k}e^{i(x_{2}-x_{1})(q_{j}-q_{k})}$
$\displaystyle=\sum_{j,k}\cos[(x_{2}-x_{1})(q_{j}-q_{k})]$

Observe that $q_{k}$ can take any value in the set $\\{-n,-n+2,\dots,n-2,n\\}$
with $\binom{n}{0},\binom{n}{1},\dots,\binom{n}{n-1},\binom{n}{n}$ occurrences
among $2^{n}$ values of $k$. Therefore the value $q_{j}-q_{k}$ is belongs the
set $\\{0,\pm 2,\pm 4,\dots,\pm 2n\\}$. We count the occurrences of those
values. There are
$\sum_{k=0}\binom{n}{k}\binom{n}{k}=\sum_{k=0}\binom{n-k}{k}\binom{n}{k}=\binom{2n}{n}$
combinations of $j$ and $k$ that leads to $q_{j}-q_{k}=0$. Similarly, the
occurrences of $\pm 2,\pm 4,\dots,\pm 2n$ are
$\binom{2n}{n-1},\binom{2n}{n-2},\binom{2n}{0}$ respectively, i.e.

(48)
$\left|\Tr(S_{2}S_{1}^{\dagger})\right|^{2}=\binom{2n}{n}+2\sum_{k=1}^{n}\binom{2n}{n-k}\cos(2k(x_{1}-x_{2}))$

##  Appendix C Distribution of coefficients

It is known that each entry from a Haar random unitary is
$\operatorname{Beta}(1,d-1)$ distributed, i.e.
$P(|u_{ij}|^{2}=x)=(d-1)(1-x)^{d-2}$ [20], which can also be derived following
the technique in Appendix A. The product $h_{j}|w_{ji}|^{2}$ therefore has the
mean of $\frac{h_{j}}{d}$ and the variance of
$\frac{h_{j}^{2}(d-1)}{d^{2}(d+1)}$. Moreover these random variables are
almost surely independent. Because
$\frac{1}{d}\sum_{j}\operatorname{Var}(h_{j}|w_{ji}|^{2})=\sum_{j}\frac{h_{j}^{2}(d-1)}{d(d+1)}$
converges to $0$ when $h_{j}\in O(1)$, the Lindeberg-Feller Central Limit
Theorem [9] (Theorem 6.13) says that the sum $\sum_{j}h_{j}|w_{ji}|^{2}$
converges in distribution to a normal random variable
$\mathcal{N}\left(1,\frac{(d-1)}{(d+1)}\sum_{j}h_{j}^{2}\right)$.

On the other hand, we expand $\bar{u}_{i1}u_{k1}$ with their real and
imaginary parts.

(49) $\displaystyle\bar{u}_{i1}u_{k1}$
$\displaystyle=(x_{i1}-iy_{i1})(x_{k1}+iy_{k1})$
$\displaystyle=(x_{i1}x_{k1}+y_{i1}y_{k1})+i(x_{i1}y_{k1}-x_{k1}y_{i1})$

We know that the real and imaginary part of any entry of a Haar random unitary
are identically independently distributed from
$\mathcal{N}\left(0,\frac{1}{2d}\right)$ with the corresponding characteristic
function $\varphi(t)=e^{-t^{2}/4d}$. Moreover, the components of $u_{i1}$ are
almost surely independent to the components of $u_{k1}$. The characteristic
function of $x_{i1}x_{k1}$ is therefore

(50) $\displaystyle\mathbb{E}[\exp(itx_{i1}x_{k1})]$
$\displaystyle=\mathbb{E}[\mathbb{E}[\exp(itx_{i1}x_{k1})|x_{k1}]]$
$\displaystyle=\mathbb{E}[\exp(-t^{2}x_{k1}^{2}/4d)]$
$\displaystyle=\int_{-\infty}^{\infty}\frac{\exp(-x_{k1}^{2}d)}{\sqrt{\pi/d}}\exp(-t^{2}x_{k1}^{2}/4d)dx_{k1}$

The last line is a Gaussian integral with
$\sigma^{2}=\frac{2d}{4d^{2}+t^{2}}$, hence evaluated to be
$1/\sqrt{1+(t/2d)^{2}}$. Thanks to the independence, the characteristic
function of $x_{i1}x_{k1}+y_{i1}y_{k1}$ is $1/(1+(t/2d)^{2})$, equal to the
characteristic function of the distribution $\operatorname{Laplace}(0,1/2d)$.
The same argument applies for $x_{i1}y_{k1}-x_{k1}y_{i1}$.

The previous argument also reveals that $\sum_{j}h_{j}\bar{w}_{ji}w_{jk}$ has
both real and imaginary parts being two linear combinations of (almost surely
independent) random variables with distribution from
$\operatorname{Laplace}(0,1/2d)$. In particular, both real and imaginary parts
have the form of $\sum_{j}q_{j}$, where
$q_{j}\sim\operatorname{Laplace}(0,|h_{j}|/2d)$ with
$\operatorname{Var}(q_{j})=h_{j}^{2}/2d^{2}$. Since
$\frac{1}{d}\sum_{j}\frac{h_{j}^{2}}{2d^{2}}=\frac{\sum_{j}h_{j}^{2}}{2d^{3}}\rightarrow
0$ when $h_{j}\in O(1)$, the Lindeberg-Feller Central Limit Theorem concludes
$\real\left\\{\sum_{j}h_{j}\bar{w}_{ji}w_{jk}\right\\}$ and
$\imaginary\left\\{\sum_{j}h_{j}\bar{w}_{ji}w_{jk}\right\\}\xrightarrow[]{\text{dist.}}\mathcal{N}\left(0,\frac{1}{2d^{2}}\sum_{j}h_{j}^{2}\right)$.