# Learning PAC-Bayes Priors
for Probabilistic Neural Networks

María Pérez-Ortiz∗, Omar Rivasplata∗, Benjamin Guedj∗†, Matthew Gleeson∗,
Jingyu Zhang∗, John Shawe-Taylor∗, Miroslaw Bober§ and Josef Kittler§
∗Centre for AI and Dept. of Computer Science, University College London (UK)

†Inria, Lille Nord-Europe research centre and Inria London Programme (France)

§Center for Vision, Speech and Signal Processing (CVSSP), University of Surrey
(UK)

###### Abstract

Recent works have investigated deep learning models trained by optimising PAC-
Bayes bounds, with priors that are learnt on subsets of the data. This
combination has been shown to lead not only to accurate classifiers, but also
to remarkably tight risk certificates, bearing promise towards self-certified
learning (_i.e._ use all the data to learn a predictor and certify its
quality). In this work, we empirically investigate the role of the prior. We
experiment on 6 datasets with different strategies and amounts of data to
learn data-dependent PAC-Bayes priors, and we compare them in terms of their
effect on test performance of the learnt predictors and tightness of their
risk certificate. We ask what is the optimal amount of data which should be
allocated for building the prior and show that the optimum may be dataset
dependent. We demonstrate that using a small percentage of the prior-building
data for validation of the prior leads to promising results. We include a
comparison of underparameterised and overparameterised models, along with an
empirical study of different training objectives and regularisation strategies
to learn the prior distribution.

##  1 Introduction

The majority of deep learning algorithms output weights of neural networks: in
recent years, a growing body of works has investigated algorithms which rather
output _probability distributions_ over the connection weights of a network,
with a number of advantages. This is evidenced _e.g._ by works inspired by
Bayesian learning (see e.g. [1, 2, 3] among many others) or by the frequentist
PAC-Bayes bounds (see e.g. [4, 5, 6, 7]). In both cases, a probability
distribution over neural network weights defines what can be called a
Probabilistic Neural Network (PNN).

Recently, PNNs learnt by optimising PAC-Bayes bounds have shown promising
results on performance guarantees, by delivering tight risk certificates
(generalisation bounds) for predictive models that are competitive compared to
standard empirical risk minimisation (see [7]). Importantly, the PNN paradigm,
coupled with PAC-Bayes bounds, is an example of self-certified learning, which
proposes to use all the available data for learning a predictor and providing
a reasonably tight numerical risk bound value that certifies the predictor’s
performance at the population level. In this case, the risk certificates can
be evaluated on a subset of the data used for training and thus do not require
a held-out test set, allowing efficient use of the available data. These
principled learning and certification strategies based on PAC-Bayes bounds
deserve further study to unfold their practical properties and limitations.

PAC-Bayes bounds (pioneered by [8, 9, 10]) are typically composed of two key
quantities: i) a term that measures the empirical performance of a so-called
‘posterior’ distribution, and ii) a term involving the divergence of the
posterior to a ‘prior’ distribution, which in most bounds is the Kullback-
Leibler (KL) divergence (we refer to [11] for a comprehensive presentation).
Thus, when using a PAC-Bayes bound as an optimisation objective, these two
terms must interact to balance fitness to data with fitness to the chosen
prior. A classical assumption underlying PAC-Bayes priors is that they must be
independent from the data on which the empirical term of the PAC-Bayes bound
is evaluated. This assumption is well-known in the PAC-Bayes literature (as
discussed by [12, 13]). Interestingly, however, the chosen prior greatly
impacts the bound via the KL term, which often amounts to the dominating
contribution to bound values (as pointed out by [14, 7]). This prominent role
of the KL term has implications both for optimisation and for risk
certification based on PAC-Bayes bounds: (i) the KL term effectively
constrains the posterior such that it cannot move too far from the prior; (ii)
the large values of the KL term when using data-independent priors or
uninformed priors suggests that these priors may not be able to give tight
risk certificates, therefore calling for data-dependent priors.

These considerations give clues on the importance of the prior distribution in
PAC-Bayes bounds: An arbitrarily chosen prior may mislead the optimisation,
since the posterior is constrained to the prior by the KL term, while a prior
representing a good solution for the problem at hand would give a better
starting point. Accordingly, some works on PAC-Bayes bounds for neural
networks have considered ways to connect PAC-Bayes priors to the data. In
particular, the recent works in [14] and [7] explored PAC-Bayes priors learnt
on a subset of the data, which does not overlap with the data used for
computing the empirical term in the PAC-Bayes bound. This way, these data-
dependent PAC-Bayes priors are in line with the classical assumption
underlying PAC-Bayes priors, while at the same time the priors have a more
sensible connection to the data when compared to arbitrarily chosen ones.

In this work, we investigate strategies for learning the PAC-Bayes prior, and
the role of the prior in obtaining PNNs with good performance by optimising
PAC-Bayes bounds.

##### Our Contributions.

  * •

We demonstrate consistency of the tightness of the risk certificates across
$6$ classification datasets.

  * •

We contribute supporting evidence on the crucial role of learning the PAC-
Bayes prior for achieving efficient neural network classifiers and tight risk
certificates.

  * •

We show that there is potential for tight risk certificates even for highly
over-parameterised architectures.

##  2 Related Literature

Bayesian Neural Networks (BNNs). BNNs are special cases of PNNs, in the sense
that BNN methods also output distributions over network weights. Works from
this front have shown very promising results, e.g. that BNNs provide an
intuitive approach to uncertainty quantification and principled
implementations of model pruning/distillation [1]. Furthermore, training
objectives inspired by the evidence lower bound (ELBO) and commonly used in
BNNs have been seen to act as an implicit regulariser. The resemblance between
the ELBO and PAC-Bayes bounds, which has been pointed out in the literature
([15, 16], among others), suggests that some of these properties may be shared
by methods inspired by PAC-Bayes bounds.

Bayes vs. PAC-Bayes. The ‘prior’ and ‘posterior’ distributions that appear in
PAC-Bayes bounds should not be confused with their Bayesian counterparts. In
PAC-Bayes bounds, what is called ‘prior’ is a reference distribution, and what
is called ‘posterior’ is an unrestricted distribution, in the sense that there
is no likelihood factor connecting them (cf. [11, 17]).

Data-dependent PAC-Bayes priors. In this line of thought, a precursor work is
that in [18]. The protocol for partitioning the data goes back to the works in
[19, 20, 21, 22].

##  3 Elements of Statistical Learning

Supervised classification algorithms receive training data
$S=((X_{1},Y_{1}),\ldots,(X_{n},Y_{n}))$ consisting of pairs that encode
inputs $X_{i}\in\mathcal{X}\subseteq\mathbb{R}^{d}$ and their labels
$Y_{i}\in\mathcal{Y}$. Classifiers $h_{w}:\mathcal{X}\to\mathcal{Y}$ are
mappings from input space $\mathcal{X}$ to label space $\mathcal{Y}$, and we
assume they are parametrised by ‘weight vectors’
$w\in\mathcal{W}\subseteq\mathbb{R}^{p}$. The quality of $h_{w}$ is given by
its risk $L(w)$, which by definition is the expected classification error on a
randomly chosen pair $(X,Y)$. However, $L(w)$ is an inaccessible measure of
quality, since the distribution that generates the data is unknown. An
accessible measure of quality is given by the empirical risk functional
$\hat{L}_{S}(w)=n^{-1}\sum_{i=1}^{n}\ell(h_{w}(X_{i}),Y_{i})$, defined in
terms of a loss function
$\ell:\mathbb{R}\times\mathcal{Y}\rightarrow[0,\infty)$. Indeed, the empirical
risk minimisation (ERM) paradigm aims to find $w\in\mathcal{W}$ that minimises
this functional for some choice of loss function (zero-one loss or a surrogate
loss).

The outcome of training a PNN is a distribution $Q$ over weight space and this
distribution depends on the sample $S$. Then, given a fresh input $X$, the
randomised classifier predicts its label by drawing a weight vector $W$ at
random from $Q$ and applying the predictor $h_{W}$ to $X$. For the sake of
simplicity, we identify the randomised predictor with the distribution $Q$
that defines it. The quality of this randomised predictor is measured by the
expected loss notions under the random draws of weights. Thus, the loss of $Q$
is given by $Q[L]=\int_{\mathcal{W}}L(w)Q(dw)$; and the empirical loss of $Q$
is given by $Q[\hat{L}_{S}]=\int_{\mathcal{W}}\hat{L}_{S}(w)Q(dw)$.

The PAC-Bayes-quadratic bound in [7] says that for any $\delta\in(0,1)$, with
probability of at least $1-\delta$ over size-$n$ i.i.d. random samples $S$,
simultaneously for all distributions $Q$ over $\mathcal{W}$ we have:

$Q[L]\leq\left(\sqrt{Q[\hat{L}_{S}]+\frac{\operatorname{KL}(Q\|Q^{0})+\log(\frac{2\sqrt{n}}{\delta})}{2n}}\right.+\left.\sqrt{\frac{\operatorname{KL}(Q\|Q^{0})+\log(\frac{2\sqrt{n}}{\delta})}{2n}}\right)^{2}.$

In this case the prior $Q^{0}$ must be chosen without any dependence on the
data $S$ on which the empirical term $\hat{L}_{S}$ is evaluated. In this work,
we use a partitioning scheme for the training data $S=S_{\mathrm{pri}}\cup
S_{\mathrm{cert}}$ such that the prior is trained on $S_{\mathrm{pri}}$, the
posterior is trained on the whole set $S$ and the risk certificate is
evaluated on $S_{\mathrm{cert}}$. See Figure 1.

Figure 1: Comparison of data partitions in traditional machine learning (where
the model is certified with a point estimate) vs self-certified learning
through PAC-Bayes bounds. Note that all $n$ training data are used by the
self-certified learning algorithm.

##  4 Learning and Certification Strategy

In a nutshell, the learning and certification strategy has three components:
(1) choose/learn a prior; (2) learn a posterior; and (3) evaluate the risk
certificate for the posterior.

###  4.1 Data-dependent PAC-Bayes priors

We experiment with Gaussian priors $Q^{0}$ with a diagonal covariance matrix
centered at i) random weights (uninformed data-free priors) and ii) learnt
weights (data-dependent priors) on a subset of the dataset which is
independent of the subset used to compute the risk certificate (see Figure 1).
In all cases, the posterior is initialised to the prior. Similar approaches
have been considered before in the PAC-Bayesian literature (we refer to [21,
7, 14].

####  4.1.1 Training objectives to learn the prior

We evaluate several training strategies to learn the PAC-Bayes prior:

  1. 1.

Deterministic training objectives: We compare learning the prior mean by non-
regularised ERM (referred to as $f_{\texttt{erm}}$) to ERM with dropout
($f_{\texttt{erm}}$ \+ dropout) and Mixup [23] ($f_{\texttt{mixup}}$), which
has been shown to act as an implicit regulariser [24]. Note that all of these
training objectives learn a deterministic neural network with a fixed setting
of weights. The prior distribution is then centered at the learnt weight
vector, and the scale of the distribution is set as a hyper-parameter.

  2. 2.

Probabilistic training objectives: We compare learning the prior distribution
via the Bayes by backprop ($f_{\texttt{bbb}}$) [1] objective inspired by the
evidence lower bound, to learning it via the objective inspired by the PAC-
Bayes-quadratic bound ($f_{\texttt{quad}}$) [7]. Since these objectives are
probabilistic, they learn the full prior distribution, and not only the mean,
as it is the case with the deterministic objectives. For both of these
objectives, we first need to set a pre-prior, for which we use uninformed
priors centered at zero or at a random weights, setting the scale as a hyper-
parameter. We then use this uninformed prior in the KL computation of the
training objective, in order to learn the data-dependent prior, which is later
used to train the posterior.

####  4.1.2 Validation of prior and posterior

The work in [7] highlighted the difficulty of avoiding overfitting of the
prior learnt by ERM and suggested the use of dropout. In this paper, we
propose to use the following strategy to monitor the performance of prior and
posterior and choose the best prior/posterior accordingly:

  1. 1.

Prior validation: We propose to select a small percentage of the prior-
building data to use as a validation set, with the objective of monitoring the
validation loss while learning the prior predictor. Then, at the end of the
prior-building stage we choose the one that performs the best in the
validation set.

  2. 2.

Posterior validation: In this case we compute the risk certificate at
different stages during posterior learning and choose the posterior with the
best risk certificate, thus not needing any additional data partition.

A summary of the partitions used at the different stages of learning and
certification is shown in Figure 1. Note that the prior validation set is
still being used to train the posterior, hence not discarding any data from
the learning process.

###  4.2 Posterior Optimisation & Certification

We now present the essential idea of training PNNs by minimizing a PAC-Bayes
upper bound on the risk. We use a recently proposed PAC-Bayes inspired
training objective [7], derived from Eq. (3) in the context of neural network
classifiers:

$f_{\mathrm{quad}}(Q)=\left(\sqrt{Q[\hat{L}^{\text{x-e}}_{S}]+\frac{\operatorname{KL}(Q\|Q^{0})+\log(\frac{2\sqrt{n}}{\delta})}{2n}}\right.+\left.\sqrt{\frac{\operatorname{KL}(Q\|Q^{0})+\log(\frac{2\sqrt{n}}{\delta})}{2n}}\right)^{2}.$

This objective is implemented using the cross-entropy loss, which is the
standard surrogate loss commonly used on these problems. Since the PAC-Bayes
bounds of Eq. (3) require the loss within [0,1], we construct a ‘bounded
cross-entropy’ loss by lower-bounding the network probabilities by a value
$p_{\mathrm{min}}>0$ (cf. [18, 7]). The empirical risk term is
$\hat{L}^{\text{x-e}}_{S}(w)$ is calculated with the ‘bounded’ version of the
cross-entropy loss.

Optimization of the objective in Eq. (4.2) entails minimizing over $Q$. By
choosing $Q$ in a parametric family of distributions, we can use the pathwise
gradient estimator [25, 26] as done by [1]. The details of the
reparametrisation strategy are outlined in [7]. Following [1], the
reparametrization we use is $W=\mu+\sigma\odot V$ with Gaussian distributions
for each coordinate of $V$. The optimization uses $\sigma=\log(1+\exp(\rho))$,
thus gradient updates are with respect to $\mu$ and $\rho$.

##### Evaluation of the Risk Certificates.

After optimising the posterior distribution over network weights through the
previously presented training objective, we compute a risk certificate on the
error of the stochastic predictor. To do so, we follow the procedure outlined
in [7], which was used before by [4] and goes back to the work of [27]. This
certification procedure uses the PAC-Bayes-kl bound. In particular, the
procedure is based on numerical inversion of the binary KL divergence, as done
by [4, 7].

##  5 Experiments

We conduct the following five experiments on the effect of PAC-Bayes priors on
PNN predictors:

  1. E1

We study the relationship between the learnt PAC-Bayes prior predictive
performance and the posterior risk certificate in a large set of experiments
on MNIST.

  2. E2

We test if the validation of prior and posterior leads to better performance,
even when the amount of data used for prior learning is effectively reduced.

  3. E3

We experiment with different trade-offs of the amount of data to learn the
prior and certify the predictor, extending the results in [7] to 5 additional
datasets, which demonstrates the role of data-dependent priors on the
tightness of the risk certificates.

  4. E4

We study the role of the number of parameters in the architecture in the risk
certificate and KL term.

  5. E5

Finally, we compare several training objectives and regularisation strategies
for learning the prior.

Table 1 shows the datasets used (all except for MNIST available at
OpenML.org), selected so as to represent a wide range of characteristics
(small vs large, low vs high dimensional and binary vs multiclass). For all
datasets except MNIST we select 20% of the data as test set (stratified with
class label). For MNIST, we use the standard data partitions.

Dataset | $n$ | $\\#f$ | $\\#c$
---|---|---|---
Spambase | 4601 | 58 | 2
Bioresponse | 3751 | 1777 | 2
Har | 10299 | 562 | 6
Electricity | 45312 | 9 | 2
Mammography | 11183 | 7 | 2
MNIST | 70000 | 784 | 10
Table 1: Datasets used: $n$ is the total number of data points, $\\#f$ the number of features and $\\#c$ the number of classes.  |  | Posterior | Prior
---|---|---|---
Dataset | Val. | R. Cert. | S. 01 err. |  KL/$n_{\texttt{cert}}$ | Post. mean 01 err. | Prior mean 01 err. | S. 01 err.
Spambase ($n_{\texttt{cert}}=1840$) | ✗ | $0.140$ | $0.082$ | $5e^{-5}$ | $0.077$ | $0.077$ | $0.083$
✓ | $0.127$ | $0.065$ | $3e^{-4}$ | $0.059$ | $0.056$ | $0.067$
Bioresponse ($n_{\texttt{cert}}=1500$) | ✗ | $0.318$ | $0.267$ | $9e^{-4}$ | $0.262$ | $0.261$ | $0.268$
✓ | $0.291$ | $0.257$ | $5e^{-5}$ | $0.250$ | $0.248$ | $0.258$
Har ($n_{\texttt{cert}}=4119$) | ✗ | $0.035$ | $0.021$ | $3e^{-4}$ | $0.017$ | $0.024$ | $0.026$
✓ | $0.037$ | $0.024$ | $8e^{-6}$ | $0.021$ | $0.020$ | $0.025$
Electricity ($n_{\texttt{cert}}=18124$) | ✗ | $0.223$ | $0.214$ | $6e^{-5}$ | $0.203$ | $0.221$ | $0.227$
✓ | $0.221$ | $0.212$ | $3e^{-5}$ | $0.203$ | $0.205$ | $0.213$
Mammography ($n_{\texttt{cert}}=4473$) | ✗ | $0.022$ | $0.015$ | $3e^{-6}$ | $0.015$ | $0.015$ | $0.016$
✓ | $0.023$ | $0.017$ | $3e^{-7}$ | $0.017$ | $0.017$ | $0.017$
MNIST ($n_{\texttt{cert}}=30000$) | ✗ | $0.034$ | $0.026$ | $1e^{-5}$ | $0.024$ | $0.025$ | $0.027$
✓ | $0.030$ | $0.027$ | $2e^{-5}$ | $0.027$ | $0.028$ | $0.031$
Table 2: Results with/without the validation set to choose the prior/posterior in an experiment using data-dependent priors. The stochastic 01 test set errors (S. 01 err.) for both prior and posterior are averaged over 10 samples. We do not include the standard deviation because they are low (around 0.001 in 6 cases and lower than 0.001 in the rest). The table also includes the 01 error of the posterior and prior mean. $n_{\texttt{cert}}$ indicates the number of data points used for the computation of the risk certificate.  Dataset | Spambase | Bioresponse | Har
---|---|---|---
| Posterior | Prior | Posterior | Prior | Posterior | Prior
Perc. of data for prior | R. Cert. | S. 01 err. | S. 01 err. | R. Cert. | S. 01 err. | S. 01 err. | R. Cert. | S. 01 err. | S. 01 err.
Data-free (0%) | $0.485$ | $0.406$ | $0.545$ | $0.531$ | $0.494$ | $0.505$ | $0.789$ | $0.758$ | $0.759$
Data-depend. (10%) | $0.213$ | $0.153$ | $0.182$ | $0.344$ | $0.315$ | $0.315$ | $0.072$ | $0.063$ | $0.077$
Data-depend. (25%) | $0.197$ | $0.133$ | $0.100$ | $0.314$ | $0.270$ | $0.266$ | $0.042$ | $0.036$ | $0.038$
Data-depend. (50%) | $\mathbf{0.127}$ | $0.065$ | $0.067$ | $\mathbf{0.291}$ | $0.258$ | $0.257$ | $0.037$ | $0.024$ | $0.025$
Data-depend. (75%) | $0.144$ | $0.056$ | $0.062$ | $0.305$ | $0.258$ | $0.258$ | $\mathbf{0.035}$ | $0.023$ | $0.024$
Data-depend. (90%) | $0.182$ | $0.062$ | $0.065$ | $0.345$ | $0.255$ | $0.255$ | $0.037$ | $0.018$ | $0.018$
ERM baseline (01 err.) | $0.045$ | $0.229$ | $0.015$
Dataset | Electricity | Mammography | MNIST
| Posterior | Prior | Posterior | Prior | Posterior | Prior
Perc. of data for prior | R. Cert. | S. 01 err. | S. 01 err. | R. Cert. | S. 01 err. | S. 01 err. | R. Cert. | S. 01 err. | S. 01 err.
Data-free (0%) | $0.441$ | $0.427$ | $0.427$ | $0.038$ | $0.023$ | $0.214$ | $0.321$ | $0.220$ | $0.891$
Data-depend. (10%) | $0.244$ | $0.234$ | $0.235$ | $0.024$ | $0.020$ | $0.019$ | $0.082$ | $0.075$ | $0.076$
Data-depend. (25%) | $0.234$ | $0.224$ | $0.232$ | $0.024$ | $0.018$ | $0.018$ | $0.047$ | $0.046$ | $0.047$
Data-depend. (50%) | $0.233$ | $0.212$ | $0.227$ | $\mathbf{0.023}$ | $0.017$ | $0.017$ | $0.030$ | $0.027$ | $0.031$
Data-depend. (75%) | $\mathbf{0.028}$ | $0.023$ | $0.208$ | $0.027$ | $0.018$ | $0.018$ | $0.028$ | $0.023$ | $0.029$
Data-depend. (90%) | $0.223$ | $0.201$ | $0.202$ | $0.031$ | $0.017$ | $0.017$ | $\mathbf{0.026}$ | $0.022$ | $0.026$
ERM baseline (01 err.) | $0.192$ | $0.017$ | $0.022$
Table 3: Results when using different percentages of the dataset for prior building and evaluation of the risk certificate. Best risk certificates for each dataset are highlighted in bold face. Dataset | Spambase | Bioresponse | Har | Electricity | Mammography
---|---|---|---|---|---
Prior objective | R. Cert. | S. 01 err. | R. Cert. | S. 01 err. | R. Cert. | S. 01 err. | R. Cert. | S. 01 err. | R. Cert. | S. 01 err.
$f_{\texttt{erm}}$ | $0.132$ | $0.088$ | $0.291$ | $0.270$ | $0.036$ | $0.024$ | $0.222$ | $0.215$ | $0.022$ | $0.017$
$f_{\texttt{erm}}$ \+ dropout | $0.126$ | $0.063$ | $0.290$ | $0.255$ | $0.033$ | $0.021$ | $0.219$ | $0.205$ | $0.022$ | $0.017$
$f_{\texttt{mixup}}$ | $0.166$ | $0.105$ | $0.299$ | $0.272$ | $0.038$ | $0.019$ | $0.241$ | $0.231$ | $0.028$ | $0.021$
$f_{\texttt{bbb}}$ (random prior) | $0.123$ | $0.068$ | $0.290$ | $0.249$ | $0.031$ | $0.020$ | $0.222$ | $0.209$ | $0.022$ | $0.017$
$f_{\texttt{bbb}}$ (zero prior) | $0.120$ | $0.060$ | $0.288$ | $0.254$ | $0.031$ | $0.020$ | $0.214$ | $0.204$ | $0.022$ | $0.017$
$f_{\texttt{quad}}$ (random prior) | $0.186$ | $0.129$ | $0.309$ | $0.272$ | $0.045$ | $0.033$ | $0.247$ | $0.237$ | $0.027$ | $0.018$
$f_{\texttt{quad}}$ (zero prior) | $0.243$ | $0.168$ | $0.309$ | $0.272$ | $0.045$ | $0.033$ | $0.247$ | $0.237$ | $0.027$ | $0.018$
Table 4: Results with different training objectives to learn the prior.
$f_{\texttt{erm}}$ and $f_{\texttt{mixup}}$ based objectives are deterministic
and only learn the mean of the prior distribution. $f_{\texttt{bbb}}$ and
$f_{\texttt{quad}}$ are probabilistic and learn the full prior distribution.

###  5.1 Experimental setup

In all experiments the models are compared under the same conditions, i.e.
weight initialisation and optimiser (vanilla SGD with momentum). The mean
parameters $\mu_{0}$ of the prior are initialised randomly from a truncated
centered Gaussian distribution with standard deviation set to
$1/\sqrt{n_{\mathrm{in}}}$, where $n_{\mathrm{in}}$ is the dimension of the
inputs to a particular layer, truncating at $\pm 2$ standard deviations. All
risk certificates are computed using the the PAC-Bayes-kl inequality, as
explained in Section 6 of [7], with $\delta=0.025$ and $\delta^{\prime}=0.01$
and $m=150000$ Monte Carlo model samples. We also report the average 01 error
of the stochastic predictor, where we randomly sample fresh model weights for
each test example $10$ times and compute the average 01 error. Input data was
standardised for all datasets.

For all experiments we performed a grid search over all hyper-parameters and
selected the run with the best risk certificate. We did a grid sweep over the
prior distribution scale hyper-parameter (i.e. standard deviation
$\sigma_{0}$) with values in $[0.1,0.05,0.04,0.03,0.02,0.01]$. For SGD with
momentum we performed a grid sweep over learning rate in
$[1\mathrm{e}^{-3},5\mathrm{e}^{-3},1\mathrm{e}^{-2}]$ and momentum in
$[0.95,0.99]$. We also performed a grid sweep over the learning rate and
momentum used for learning the prior (testing the same values as before). The
dropout rate used for learning the prior was selected from
$[0.01,0.05,0.1,0.2]$. The grid sweep is done for experiment E2. In the
subsequent experiments, we use the same best performing hyper-parameters from
E2. We experiment with fully connected neural networks (FCN) with 2/3 layers
(excluding the ‘input layer’) and 100 units per hidden layer (unless specified
otherwise). ReLU activations are used in each hidden layer. For learning the
prior we ran the training for 500 epochs (except for MNIST for which we ran
100). Posterior training was run for 100 epochs. We use a training batch size
of $250$. PyTorch code will be released anonymously at the repository
associated to this
project111https://anonymous.4open.science/r/pacbayespriors-F355/.

###  5.2 E1: The role of the PAC-Bayes prior in PNNs

Figure 2: Comparison of several metrics vs the risk certificate for more than 600 runs with different hyper-parameters on MNIST as done in [7]. We use a reduced subset of MNIST for these experiments (10% of training data).  | FCN Architecture | Posterior
---|---|---
Dataset | #Params | Neurons | H. Lay. | R. Cert. | S. 01 err. |  KL/$n_{\texttt{cert}}$
Spamb. | $1404$ | $10$ | $2$ | $0.130$ | $0.088$ | $1e^{-4}$
$32004$ | $100$ | $2$ | $0.132$ | $0.088$ | $9e^{-5}$
$792004$ | $600$ | $2$ | $0.120$ | $0.060$ | $4e^{-4}$
$1512004$ | $600$ | $3$ | $0.120$ | $0.060$ | $6e^{-4}$
Bioresp. | $7128$ | $3$ | $2$ | $0.521$ | $0.483$ | $3e^{-8}$
$35784$ | $10$ | $2$ | $0.287$ | $0.254$ | $9e^{-6}$
$375804$ | $100$ | $2$ | $0.291$ | $0.270$ | $4e^{-5}$
$395804$ | $100$ | $3$ | $0.290$ | $0.248$ | $2e^{-5}$
$2854804$ | $600$ | $2$ | $0.283$ | $0.257$ | $3e^{-5}$
Har | $4588$ | $4$ | $2$ | $0.405$ | $0.376$ | $7e^{-6}$
$11572$ | $10$ | $2$ | $0.040$ | $0.024$ | $3e^{-5}$
$133612$ | $100$ | $2$ | $0.036$ | $0.024$ | $9e^{-6}$
$153612$ | $100$ | $3$ | $0.037$ | $0.027$ | $5e^{-5}$
$1401612$ | $600$ | $2$ | $0.038$ | $0.023$ | $3e^{-5}$
Elect. | $424$ | $10$ | $2$ | $0.242$ | $0.230$ | $7e^{-6}$
$22204$ | $100$ | $2$ | $0.222$ | $0.215$ | $2e^{-5}$
$733204$ | $600$ | $2$ | $0.220$ | $0.209$ | $1e^{-5}$
$1453204$ | $600$ | $3$ | $0.211$ | $0.203$ | $6e^{-5}$
Mammo. | $384$ | $10$ | $2$ | $0.026$ | $0.018$ | $3e^{-7}$
$21804$ | $100$ | $2$ | $0.022$ | $0.017$ | $4e^{-7}$
$730804$ | $600$ | $2$ | $0.025$ | $0.018$ | $6e^{-5}$
$1450804$ | $600$ | $3$ | $0.027$ | $0.018$ | $9e^{-5}$
MNIST | $16100$ | $10$ | $2$ | $0.071$ | $0.068$ | $2e^{-4}$
$178820$ | $100$ | $2$ | $0.030$ | $0.027$ | $2e^{-4}$
$1672820$ | $600$ | $2$ | $0.028$ | $0.023$ | $9e^{-5}$
$2392820$ | $600$ | $3$ | $0.029$ | $0.024$ | $2e^{-4}$
Table 5: Results with architectures with different numbers of parameters (both
under and over parametrised with respect to the dataset size [28]).

The work in [7] showed that PAC-Bayes bounds can be used not only as training
objectives but also for model selection. Specifically, the authors compared
the risk certificates obtained to the test stochastic 01 error in a large-
scale hyper-parameter grid search experiment involving more than 600 runs on
MNIST, showing a strong relationship which motivated using the risk
certificates for model selection.

We begin by replicating this experiment to study the role of the prior on the
risk certificate. To do this, we perform the same grid search with the
architecture reported in [7] over 6 hyper-parameters: prior scale, dropout
rate, and the learning rate and momentum both for learning the prior and the
posterior. These experiments use a reduced MNIST training set (i.e. 10%),
while the test set is maintained. Figure 2 shows the relationship between
three metrics of interest and the risk certificate. We start by comparing the
test 01 loss of the stochastic predictor to the risk certificate as done in
[7] (right plot). We observe a similarly strong relationship. The middle plot
shows the relationship between the deterministic test 01 loss of the prior
mean learnt by ERM and the risk certificate. Interestingly, this plot shows
that when the prior is competitive, then the risk certificate is as well.
However, good posteriors and risk certificates can be achieved even if the
prior is not very competitive. Studying the hyper-parameters in this
experiment, we observe that both the standard deviation used in the prior and
the use of dropout when learning the prior have a big impact on the posterior
risk certificate, independently of the prior’s test set 01 error. Finally, the
right part of the plot shows the relationship between the KL term in the bound
and the risk certificate. This plot demonstrates the significant influence of
the KL term (and thus the prior) in the final risk certificate.

###  5.3 E2: Validation of prior and posterior

We now test the use of the validation strategy to choose the best fitted prior
and posterior. To do so, we use 50% of the training data to build the prior,
and 5% of those 50% for validation. The results of this experiment can be seen
in Table 2, where the column named _Validation_ indicates whether we use the
above mentioned strategy to select the best prior and posterior. We note that
even when discarding a small part of the prior set, this strategy improves the
risk certificate in most cases (4 out of 6 datasets). This seems to lead to
much smaller values of the KL term. However, the posterior stochastic 01 error
does not improve in all cases (worsens in 3 out of the 6 datasets and improves
in the rest). All the subsequent experiments in the present paper use this
validation strategy to choose the prior and posterior. Table 2 also includes a
comparison of the posterior and prior in terms of the stochastic 01 error and
the deterministic 01 error of the mean of the distribution. The results show a
consistent improvement of the stochastic error of the posterior over the
stochastic error of the prior, while often the error of the mean of the
distribution is kept similar or slightly worsened.

###  5.4 E3: Percentage of data devoted to prior learning

Table 3 shows the results obtained when using different percentages of the
dataset to learn the prior by ERM (the rest of the dataset being used to
certify the predictor). The first relevant conclusion is that data-free priors
(represented in the table as 0%) are far from competitive, thus justifying
that the tightness of the risk certificates are indeed due to the use of data-
dependent priors. The table also shows that the amount of data for the prior
that leads to the best risk certificates is dataset dependent, with 50% and
75% approaching the optimal trade-off. Interestingly, using more data for the
prior (e.g. 90%) improves the performance of the predictor but often harms the
risk certificates, as the number of data points used to evaluate the bound
decreases. This is in line with the findings of [14]. We note that the best
trade-off on the amount of data for prior building and evaluation of the risk
certificate may depend on the size of the dataset. For example, in the case of
electricity, which is the second largest dataset, using 90% of the data for
building the prior does not harm the risk certificate considerably (0.2%), as
we will still be using a reasonable amount to evaluate the bound. It is also
worth noting the tightness of the risk certificate in this case (2.2% risk
certificate vs 2.1% stochastic test set 01 error). The same is applicable to
MNIST, where the risk certificates improve when using 90% of the data for
prior building (since we will still be using 6000 examples for computing the
risk certificate). However, in the case of Bioresponse and Spambase (smallest
datasets) the risk certificates worsen by 6.4% and 4.6% respectively when
comparing 50% of the data for prior building to 90%. The table also includes
the results of the ERM baseline. For ERM, we used the same range for
optimising the learning rate, momentum and dropout rate. However, given that
in this case we do not have a risk certificate we need to set aside some data
for validation and hyper-parameter tuning. In all cases we set 5% of the data
for validation. Note that this predictor is deterministic, which usually shows
slightly better performance than the stochastic counterpart (as also discussed
by [18, 7]). The results on these datasets show a slightly larger gap between
PNNs and standard neural networks learnt by ERM when compared to the results
on MNIST and CIFAR-10 in [7]. Nonetheless, the risk certificates are tight.

###  5.5 E4: Differently parametrised architectures

The experiments in [7] also showed that the KL term decreased with the depth
of the architecture in CIFAR-10. This is a noteworthy observation and possibly
non intuitive initially, as there are many more parameters in the computation
of the KL for deep architectures and one would expect the KL term to increase
significantly with the number of parameters. However, the authors conjectured
that this may be because in a higher-dimensional weight space, the weight
updates have smaller euclidean norms (see the discussion in [29]), hence the
smaller KL. We further validate with our experiments that the KL term does not
increase with the number of parameters in the architecture. The results are
shown in Table 5, where we test multiple fully connected architectures
architectures with different amount of hidden layers and neurons. Our results
demonstrate that highly over parameterised models show a KL term of similar
factor than those architectures with less parameters. This could be because
over parameterised models have optimization properties qualitatively different
from those that are under parameterised [28] or that over parameterised priors
are so efficient that posteriors do not need to travel far. Either way, we
believe this to be a promising result, demonstrating that there is potential
for tight risk certificates even for highly over parameterised architectures.

###  5.6 E5: Training objectives for learning the prior

The results with different training objectives to learn the prior distribution
are gathered in Table 4. The first conclusion is that regularisation of the
prior using dropout improves the performance of $f_{\texttt{erm}}$
consistently. However, mixup regularisation seems to worsen the risk
certificates and the stochastic error, at least for tabular datasets and FCNs.
Note that this objective has an additional hyper-parameter for which we do the
same grid sweep as suggested in [23] for MNIST. The Bayesian inspired learning
objective $f_{\texttt{bbb}}$, which learns the full prior distribution,
improves the results over $f_{\texttt{erm}}$, generally both in terms of the
risk certificate and the stochastic 01 error. In this sense, uninformed priors
centered at zero, as suggested by [1], seem to perform better than random
uninformed priors. Note that $f_{\mathrm{bbb}}$ requires an additional sweep
over a KL trade-off coefficient (which attenuates the effect of the KL in the
training objective) which was done with values in
$[1\mathrm{e}^{-6},1\mathrm{e}^{-5},\ldots,1\mathrm{e}^{-2}]$, see [1]. The
results using the original $f_{\texttt{quad}}$ objective for learning the
prior led to vacuous bounds, so we experimented with a version of this
objective that also attenuates the KL (see [7]). The results in Table 4
include the results of $f_{\mathrm{quad}}$ with the same KL trade-off
coefficient sweep as in $f_{\mathrm{bbb}}$. In this case, the risk
certificates are not vacuous although not as tight as when learning the prior
with $f_{\mathrm{bbb}}$. We conjecture that $f_{\texttt{quad}}$ still puts
much more emphasis on the KL term than $f_{\mathrm{bbb}}$, which in this case
measures the divergence to an uninformed prior, which restricts the learning
significantly. Given that mixup regularisation was originally proposed for
vision datasets, convolutional architectures and data starvation scenarios
[30], we set to experiment with its use in MNIST with a CNN with $4$ layers
where we remove training data at random, reducing the training set to 25% of
the original. The results are in Table 6: note that dropout also helps in this
case to learn better PAC-Bayes priors. Additionally, $f_{\texttt{mixup}}$,
which can be seen as a type of data augmentation, leads to very promising
results in this case. Specially, it can be seen that $f_{\texttt{mixup}}$ is
data efficient, in the sense that we may need less data to build the prior to
achieve tight risk certificates and accurate classifiers (compare the results
of $f_{\texttt{mixup}}$ with 10% of the data to build the prior to
$f_{\texttt{erm}}$ with 50% of the data for the prior).

Perc. of data for prior | Prior objective | R. Cert. | S. 01 err.
---|---|---|---
10% | $f_{\texttt{erm}}$ | $0.110$ | $0.074$
$f_{\texttt{erm}}$ \+ dropout | $0.106$ | $0.067$
$f_{\texttt{mixup}}$ | $0.065$ | $0.041$
25% | $f_{\texttt{erm}}$ | $0.089$ | $0.056$
$f_{\texttt{erm}}$ \+ dropout | $0.044$ | $0.029$
$f_{\texttt{mixup}}$ | $0.041$ | $0.025$
50% | $f_{\texttt{erm}}$ | $0.061$ | $0.042$
$f_{\texttt{erm}}$ \+ dropout | $0.026$ | $0.020$
$f_{\texttt{mixup}}$ | $0.025$ | $0.019$
Table 6: Results with deterministic training objectives to learn the prior on
MNIST in a data starvation experiment (reducing the training set in MNIST to
25% of the original).

Figure 3: Consistency of the tightness of the risk certificates across
datasets.

##### Tightness across datasets.

Risk certificates are consistently tight across multiple datasets, as shown in
Figure 3, where we plot all results reported in this paper together.

##  6 Conclusion

This work empirically studies learning PAC-Bayes priors from data. Our results
show that data-dependent priors lead to consistently tight risk certificates
in $6$ datasets for probabilistic neural classifiers and over parameterised
networks, setting a stepping stone towards achieving self-certified learning.
We compare a wide range of training objectives for learning the prior
distribution, and show that regularisation of the prior is important and that
Bayesian inspired learning objectives hold potential for learning appropriate
priors, in this case learning the full prior distribution, as opposed to only
the mean. Our results also demonstrate that data augmentation may be desirable
during prior learning.

##  7 Acknowledgments

We gratefully acknowledge support and funding from the U.S. Army Research
Laboratory and the U. S. Army Research Office, and by the U.K. Ministry of
Defence and the U.K. Engineering and Physical Sciences Research Council
(EPSRC) under grant number EP/R013616/1.

We warmly acknowledge the Department of Computer Science and the AI Centre at
University College London, and the Department of Electric and Electronic
Engineering at the University of Surrey, for providing friendly and
stimulating work environments.

Omar Rivasplata gratefully acknowledges sponsorship from DeepMind for
conducting research studies in machine learning at University College London.

## References

  * [1] Charles Blundell, Julien Cornebise, Koray Kavukcuoglu, and Daan Wierstra.  Weight uncertainty in neural networks.  In International Conference on Machine Learning [ICML], pages 1613–1622. PMLR, 2015.
  * [2] Kazuki Osawa, Siddharth Swaroop, Mohammad Emtiyaz Khan, Anirudh Jain, Runa Eschenhagen, Richard E. Turner, and Rio Yokota.  Practical deep learning with Bayesian principles.  In Advances in Neural Information Processing Systems [NeurIPS], pages 4289–4301, 2019.
  * [3] Andrew Gordon Wilson and Pavel Izmailov.  Bayesian deep learning and a probabilistic perspective of generalization.  In Hugo Larochelle, Marc’Aurelio Ranzato, Raia Hadsell, Maria-Florina Balcan, and Hsuan-Tien Lin, editors, Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual, 2020.
  * [4] Gintare Karolina Dziugaite and Daniel M. Roy.  Computing nonvacuous generalization bounds for deep (stochastic) neural networks with many more parameters than training data.  In Conference on Uncertainty in Artificial Intelligence [UAI]. AUAI Press, 2017.
  * [5] Gaël Letarte, Pascal Germain, Benjamin Guedj, and François Laviolette.  Dichotomize and generalize: PAC-Bayesian binary activated deep neural networks.  In Advances in Neural Information Processing Systems [NeurIPS], pages 6872–6882, 2019.
  * [6] Felix Biggs and Benjamin Guedj.  Differentiable pac-bayes objectives with partially aggregated neural networks.  CoRR, abs/2006.12228, 2020.
  * [7] María Pérez-Ortiz, Omar Rivasplata, John Shawe-Taylor, and Csaba Szepesvári.  Tighter risk certificates for neural networks.  CoRR, abs/2007.12911, 2020.
  * [8] John Shawe-Taylor and Robert C. Williamson.  A PAC analysis of a Bayesian estimator.  In Computational Learning Theory [COLT], pages 2–9. ACM, 1997\.
  * [9] David A. McAllester.  Some PAC-Bayesian theorems.  In Computational Learning Theory [COLT], pages 230–234. ACM, 1998\.  Also one year later in Machine Learning 37(3), pages 355–363, 1999.
  * [10] David A. McAllester.  PAC-Bayesian model averaging.  In Computational Learning Theory [COLT], pages 164–170. ACM, 1999\.
  * [11] Benjamin Guedj.  A primer on PAC-Bayesian learning.  In Emmanuel Breuillard, editor, Congrès de la Société Mathématique de France, Collection SMF, volume 33, 2019.
  * [12] Olivier Catoni.  Statistical Learning Theory and Stochastic Optimization: Ecole d’Eté de Probabilités de Saint-Flour XXXI-2001, volume 1851 of Lecture Notes in Mathematics.  Springer, 2004.
  * [13] Olivier Catoni.  PAC-Bayesian Supervised Classification: The Thermodynamics of Statistical Learning, volume 56 of IMS Lecture Notes-Monograph Series.  Institute of Mathematical Statistics, 2007.
  * [14] Gintare Karolina Dziugaite, Kyle Hsu, Waseem Gharbieh, Gabriel Arpino, and Daniel M. Roy.  On the role of data in PAC-Bayes bounds.  In International Conference on Artificial Intelligence and Statistics [AISTATS], pages 604–612. PMLR, 2021.
  * [15] Pierre Alquier, James Ridgway, and Nicolas Chopin.  On the properties of variational approximations of Gibbs posteriors.  Journal of Machine Learning Research, 17(236):1–41, 2016.
  * [16] Alessandro Achille and Stefano Soatto.  Emergence of invariance and disentanglement in deep representations.  Journal of Machine Learning Research, 19(50):1–34, 2018.
  * [17] Omar Rivasplata, Ilja Kuzborskij, Csaba Szepesvári, and John Shawe-Taylor.  PAC-Bayes analysis beyond the usual bounds.  In Advances in Neural Information Processing Systems [NeurIPS], pages 16833–16845, 2020.
  * [18] Gintare Karolina Dziugaite and Daniel M. Roy.  Data-dependent PAC-Bayes priors via differential privacy.  In Advances in Neural Information Processing Systems [NeurIPS], pages 8440–8450, 2018.
  * [19] Amiran Ambroladze, Emilio Parrado-Hernández, and John Shawe-Taylor.  Tighter PAC-Bayes bounds.  Advances in Neural Information Processing Systems [NIPS], pages 9–16, 2007.
  * [20] Guy Lever, François Laviolette, and John Shawe-Taylor.  Distribution-dependent PAC-Bayes priors.  In International Conference on Algorithmic Learning Theory [ALT], pages 119–133. Springer, 2010.
  * [21] Emilio Parrado-Hernández, Amiran Ambroladze, John Shawe-Taylor, and Shiliang Sun.  PAC-Bayes bounds with data dependent priors.  Journal of Machine Learning Research, 13(112):3507–3531, 2012.
  * [22] Guy Lever, François Laviolette, and John Shawe-Taylor.  Tighter PAC-Bayes bounds through distribution-dependent priors.  Theoretical Computer Science, 473:4–28, 2013.
  * [23] Hongyi Zhang, Moustapha Cissé, Yann N. Dauphin, and David Lopez-Paz.  mixup: Beyond empirical risk minimization.  In 6th International Conference on Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track Proceedings. OpenReview.net, 2018.
  * [24] Linjun Zhang, Zhun Deng, Kenji Kawaguchi, Amirata Ghorbani, and James Zou.  How does mixup help with robustness and generalization?  In 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021. OpenReview.net, 2021.
  * [25] Robert Price.  A useful theorem for nonlinear devices having Gaussian inputs.  IRE Transactions on Information Theory, 4(2):69–72, 1958.
  * [26] Martin Jankowiak and Fritz Obermeyer.  Pathwise derivatives beyond the reparameterization trick.  In International Conference on Machine Learning [ICML], pages 2240–2249. PMLR, 2018.
  * [27] John Langford and Rich Caruana.  (Not) bounding the true error.  In Advances in Neural Information Processing Systems [NIPS], pages 809–816, 2001.
  * [28] Mikhail Belkin, Daniel Hsu, Siyuan Ma, and Soumik Mandal.  Reconciling modern machine-learning practice and the classical bias–variance trade-off.  Proceedings of the National Academy of Sciences, 116(32):15849–15854, 2019.
  * [29] Chiyuan Zhang, Samy Bengio, and Y. Singer.  Are all layers created equal?  ArXiv, abs/1902.01996, 2019.
  * [30] Masanari Kimura.  Why mixup improves the model performance, 2021.


# Third-harmonic generation of spatially structured light in a quasi-
periodically poled crystal

Yan-Chao Lou  Zi-Mo Cheng  Zhi-Hong Liu  Yu-Xiang Yang  Zhi-Cheng Ren
Jianping Ding  Xi-Lin Wang<EMAIL_ADDRESS>National Laboratory of Solid
State Microstructures and School of Physics, Nanjing University, Nanjing
210093, China  Collaborative Innovation Center of Advanced Microstructures,
Nanjing University, Nanjing 210093, China  Hui-Tian Wang<EMAIL_ADDRESS>National Laboratory of Solid State Microstructures and School of Physics,
Nanjing University, Nanjing 210093, China  Collaborative Innovation Center of
Advanced Microstructures, Nanjing University, Nanjing 210093, China

(September 3, 2024)

###### Abstract

Nonlinear optical processes of spatially structured light, including optical
vortex and vector optical fields, have stimulated a lot of interesting
physical effects and found a variety of important applications ranging from
optical imaging to quantum information processing. However, high harmonic
generation of vector optical fields with space-varying polarization states is
still a challenge. Here we demonstrate third harmonic generation of spatially
structured light including vector optical fields, in a nonlinear Sagnac
interferometer containing a carefully designed quasi-periodically poled
potassium titanyl phosphate for the first time. The experimental results are
in good agreement with the theoretical predictions. Our results will enable to
manipulate spatially structured light or photons at new wavelengths and
carrying higher orbital angular momentum. Our approach has the potential
applications for the research of optical skyrmions and may open up new
opportunities to produce spatially structured entangled photons for quantum
communication and computation.

Spatially structured light Rubinsztein-Dunlop2017 ; Wang2020 ; Forbes2021  has
attracted broad interest because of its novel effects and important
applications. For example, the subluminal effect has been observed
Giovannini2015 . Various novel optical fields could be created, such as
optical needle Wang2008 , optical polarization Möbius strips Bauer2015 , and
spatiotemporal optical vortices Chong2020 . Spatially structured light has
been found useful in many important applications from classical optics
Andrews2008  to quantum optics Fobres2019 .

Optical nonlinearity of spatially structured light has become an important and
hot topic in recent years. As a typical spatially structured light, an optical
vortex with helical phase front of $\exp(jm\varphi)$ (where $\varphi$ is the
azimuthal angle in polar coordinate system and $m$ is the topological charge)
could carry an orbital angular momentum (OAM) of $m\hbar$ per photon Allen1992
; Padgett2017 . A variety of nonlinear interactions involving OAM have been
demonstrated, including second harmonic generation (SHG) Dholakia1996 ;
Bloch2012 ; Shemer2013 ; Li2013 ; Zhou2014 ; Chaitanya2015 ; Ni2016 ; Chen2020
; Tang2020 ; Gui2021 ; Hancock2021 , third harmonic generation (THG) Fang2016
; Xu2018 , high-order harmonic generation Gauthier2017 , and spontaneous
parametric down-conversion (SPDC) Mair2001 ; Leach2010 ; Dada2011 . The
introduction of spatial phase structure in nonlinear interaction will result
in novel effects and promote new applications. For instance, besides the
energy and linear momentum conservations in nonlinear process, the OAM
conservation Dholakia1996 ; Mair2001 ; Leach2010 ; Dada2011 ; Li2013 ;
Fang2016 ; Xu2018 ; Gui2021 ; Hancock2021  has been explored in many nonlinear
interactions involving optical vortices. The applications of nonlinear
interaction with OAMs are prolific in many realms such as optical imaging
Qiu2018 , holography Fang2021 , and high-dimensional quantum information
processing Erhard2020 .

Vector optical fields, as another typical spatially structured light, have
space-varying polarization structures. Such kind of spatially structured light
can manipulate more photonic degrees of freedom simultaneously and then lead
to a variety of novel functionalities, such as spin-orbit photonics
Cardano2015 , alignment-free communication Aolita2007 ; DAmbrosio2012 , and
enhancing channel capacity Barreiro2008 ; Zhu2021 . It is more challenging to
explore the nonlinear process of vector optical fields, which involves various
polarization components, while the phase-matching condition generally
restricts nonlinear interaction only for a fixed combination of polarization
states. The SHG of vector optical fields has been reported by a nonlinear
interferometer YangOL2019  or two orthogonally configured nonlinear crystals
LiuOL2018  to achieve the nonlinear interaction between two components of
vector optical field. However, the experimental realization of THG of vector
optical field is still a challenge.

In this letter, we demonstrate the first THG of vector optical field in a
nonlinear Sagnac interferometer with a quasi-periodically poled potassium
titanyl phosphate (QPPKTP). The quasi-periodic structure of second-order
nonlinear coefficient can provide more abundant reciprocal vectors, which is a
flexible approach for achieving multiple second-order nonlinear interactions
in a single crystal. For example, the THG was demonstrated experimentally in
one-dimensional nonlinear quasicrystals Zhu1997 ; Zhang2001 ; Vernay2020  and
a theoretical design of two-dimensional photonic quasicrystals for multiple
nonlinear frequency conversion was proposed Lifshitz2005 . Here, the single
QPPKTP produces third harmonic (TH) wave by implementing two second-order
nonlinear optical processes simultaneously: the SHG of a fundamental wave
(FW), and the sum-frequency generation (SFG) of the FW and the second harmonic
(SH) wave.

Here we devote to the THG of vector optical field. To achieve our goal, a good
choice is to use a QPPKTP, which is designed by following the projection
method Zhang2001 ; Zia1985 . To use the largest second-order nonlinear
coefficient, the QPPKTP is designed as a type-0 phase-matching scheme, i.e.,
all the interacting waves are vertically polarized (denoted by $|V\rangle$).
As shown in Fig. 1(a), our QPPKTP consists of two building blocks A and B.
Each block (A or B) contains a pair of antiparallel domains. The width of A
(B) is $l_{A}\\!=\\!l_{A^{+}}+l_{A^{-}}$ ($l_{B}\\!=\\!l_{B^{+}}+l_{B^{-}}$).
The widths of domains, the projection angle $\theta$ (or
$\gamma\\!=\\!\tan\theta$), and the arrangement of the QPPKTP (with a length
of 5.532 mm) are listed in Fig. 1(a)

Figure 1: THG of vector optical field. (a) QPPKTP consists of two building
blocks A and B, and each block contains a pair of antiparallel domains. (b)
Configuration for THG of vector field in a double-pass nonlinear Sagnac
interferometer with a QPPKTP. There are four units: the generation and
detection units of the FW vector field, the THG unit, the TH detection unit.
Filter is used to block the SH at 780 nm.

Under the non-depletion approximation and the phase-matching conditions, when
the input FW field is a $|V\rangle$-polarized vortex with a topological charge
of $m$ (its state is written as $|V\rangle|m\rangle$), the generated
$|V\rangle$-polarized TH field in the type-0 QPPKTP should be written as
$|V\rangle|3m\rangle$ Zhu1997 . Clearly, the topological charge of the TH
vortex field is triple that of the FW field because three FW photons are
converted into one TH photon and the OAM is conserved in this THG process
Fang2016 ; Xu2018 . Our QPPKTP is only valid for the $|V\rangle$-polarized
case. For instance, for a FW vector field containing both $|V\rangle$\- and
$|H\rangle$-polarized (horizontally polarized) vortices, its state can be
written as

$|\phi\rangle=\alpha|H\rangle|\\!+\\!m\rangle+\beta|V\rangle|\\!-\\!m\rangle,$
(1)

where $\alpha$ and $\beta$ are the normalized factors with
$\alpha^{2}+\beta^{2}=1$. For convenience, we choose the two vector field
expression bases as $|H\rangle$ and $|V\rangle$, which are the two bases for
the vector field generation element utilized in our experiment: polarization
interferometric prism (PIP) Ren2021 . A single-pass configuration with only
one QPPKTP cannot realize the THG of the two components for a vector field,
hence a double-pass configuration should be a good solution. Here we utilize a
Sagnac interferometer with the QPPKTP shown in Fig. 1(b) to implement a
double-pass THG. For the FW vector field, its $|H\rangle$-polarized
($|V\rangle$-polarized) component is transmitted (reflected) by a polarizing
beam splitter (PBS), and then propagates along the clockwise
(counterclockwise) loop. It is of great importance to insert a half-wave plate
(HWP) at 45∘ into the Sagnac interferometer, which converts the
$|H\rangle$-polarized component into the $|V\rangle$-polarized one to generate
the $|V\rangle$-polarized TH field. As PBS and HWP in the Sagnac
interferometer should also work for the TH field, both are designed to be
dual-wavelength PBS (DPBS) and dual-wavelength HWP (DHWP). And reflection will
invert the sign of OAM or topological charge of the optical vortex Ren2021 ,
suggesting that the odd (even) reflections will invert (keep) the sign of
topological charge of the optical vortex.

In the counterclockwise loop, the $|V\rangle$-polarized FW vortex undergoes
only two reflections (by the DPBS and a mirror) and then generates the
$|V\rangle$-polarized TH vortex in the QPPKTP; after one reflection (by a
mirror) and passing through the DHWP at $45^{\circ}$, the TH vortex becomes
$|H\rangle$-polarized and its sign of OAM is also reversed; finally it outputs
from the Sagnac interferometer through the DPBS. In the clockwise loop, the
$|H\rangle$-polarized FW vortex is converted to be $|V\rangle$-polarized by
the DHWP at $45^{\circ}$; after one reflection by a mirror (the sign of OAM is
reversed), it is incident into the QPPKTP; the generated $|V\rangle$-polarized
TH vortex experiences two reflections (by a mirror and the DPBS) and then
outputs from the Sagnac interferometer. As a result, the two orthogonally
polarized components of the TH fields generated in the two loops are combined
at the DPBS to generate the TH vector field propagating along the direction
opposite to the FW one. As stated above, therefore, when the FW vector field
shown in Eq. (1) is incident into the Sagnac interferometer, correspondingly,
the output TH vector field from the THG unit should have the following form of
state

$\displaystyle|\phi^{\prime}\rangle=(\alpha^{6}+\beta^{6})^{-1/2}(\alpha^{3}|V\rangle|\\!-\\!3m\rangle+\beta^{3}|H\rangle|\\!+\\!3m\rangle).$
(2)

Clearly, the $|H\rangle$-polarized ($|V\rangle$-polarized) FW vortex at the
state $|H\rangle|\\!+\\!m\rangle$ ($|V\rangle|\\!-\\!m\rangle$) generates the
$|V\rangle$-polarized ($|H\rangle$-polarized) TH vortex at the state
$|V\rangle|\\!-\\!3m\rangle$ ($|H\rangle|\\!+\\!3m\rangle$) from the Sagnac
interferometer. The generated TH vector field is reflected by a dichroic
mirror (DM) to be separated from the FW one and then enters the TH detection
unit in Fig. 1(b).

In experiment, we can generate the FW vector fields by the FW generation unit
shown in Fig. 1(b). A femtosecond (fs) laser (FF ULTRA 1560, TOPTICA Photonics
Inc.) with a central wavelength of 1560 nm, a repetition frequency of 80 MHz,
and a pulse duration of 200 fs is incident on a spatial light modulator (SLM,
model P1920, Meadowlark Optics Inc.). The SLM is loaded the programmable
holographic grating to produce the FW vortex field with variable topological
charge, which is incident into PIP after passing through a HWP, finally the FW
vector fields can be generated (the details can be found in Ref. Ren2021 ). Of
course, our FW generation unit generates easily the scalar optical vortices
with homogeneous polarization. It should be pointed out that the profiles of
the FW fields we generated are very similar to the zero-radial-index Laguerre-
Gaussian (LG) beams.

Firstly, we would like to demonstrate the THG of scalar FW optical vortex in
the state of $|P\rangle|m\rangle$, where $|P\rangle$ denotes the polarization
state, such as $|H\rangle$-polarized ($|V\rangle$-polarized) state and right-
handed (left-handed) circularly polarized state
$|R\rangle\\!=\\!(|H\rangle-j|V\rangle)/\sqrt{2}$
($|L\rangle\\!=\\!(|H\rangle+j|V\rangle)/\sqrt{2}$). Here we prepare three FW
optical vortices in the states of $|V\rangle|\\!-\\!1\rangle$,
$|L\rangle|\\!-\\!2\rangle$ and $|H\rangle|\\!-\\!3\rangle$. The $|H\rangle$\-
or $|V\rangle$-polarized FW optical vortex will propagates along the clockwise
or counterclockwise loop of the nonlinear Sagnac interferometer independently,
while the circularly polarized FW optical vortex will pass through both the
clockwise and counterclockwise loops simultaneously.

To detect the FW optical vortices, a flip mirror is utilized to reflect it
from the main optical path and another mirror is utilized to preserve its
topological charge before entering the FW detection unit shown in Fig. 1(b).
The intensity patterns of the FW vortices at 1560 nm are measured by CCD1
(SP907-1550, Ophir-Spiricon Inc.). The experimentally measured patterns of the
FW vortex fields with the states of $|V\rangle|\\!-\\!1\rangle$,
$|L\rangle|\\!-\\!2\rangle$ and $|H\rangle|\\!-\\!3\rangle$ are shown in the
first, third and fifth columns of Fig. 2. CCD2 (LBP2-HR-VIS2, Newport
Corporation) is used to record the intensity patterns of the TH fields at 520
nm by using the TH detection unit shown in Fig. 1(b). The experimentally
measured patterns of the TH vortices with theoretically predicted states of
$|H\rangle|\\!+\\!3\rangle$, $|R\rangle|\\!+\\!6\rangle$ and
$|V\rangle|\\!+\\!9\rangle$ are shown in the second, fourth and sixth columns
of Fig. 2, respectively. The topological charges are measured by a tilt lens
(the last row of Fig. 2), in which the stripe patterns along the diagonal or
anti-diagonal direction indicate the positive or negative topological charge,
and the number of dark stripes give the absolute value of the topological
charge VaityPLA2013 . Clearly, the polarization states and topological charges
of the TH vortex fields are in good agreement as predicted by theory as shown
in Eq. (2).

Figure 2: Experimental results for THG of scalar vortices. All the red
patterns show the experimental results for the scalar FW vortices, while all
the green patterns show the corresponding TH results. First row shows the
total intensity patterns for the FW and TH scalar vortices. Second-fifth rows
are the measured intensity patterns of the components under the bases of
$|H\rangle$, $|V\rangle$, $|R\rangle$ and $|L\rangle$. The last row shows the
patterns of the measured vortices after passing through a tilt lens to reveal
topological charges.

Next, we would like to investigate the THG of vector fields. Using the FW
generation unit in Fig. 1(b), we prepare the FW vector fields in the three
states

$\displaystyle|\phi_{1}\rangle$
$\displaystyle=2^{-1/2}\left(|H\rangle|\\!-\\!1\rangle+|V\rangle|\\!+\\!1\rangle\right),$
(3) $\displaystyle|\phi_{2}\rangle$
$\displaystyle=2^{-1/2}\left(|H\rangle|\\!-\\!4\rangle-|V\rangle|\\!+\\!4\rangle\right),$
(4) $\displaystyle|\phi_{3}\rangle$ $\displaystyle=\sin
35^{\circ}|H\rangle|\\!-\\!1\rangle+\cos 35^{\circ}|V\rangle|\\!+\\!1\rangle.$
(5)

With Eq. (2), the three FW vector fields in Eqs. (3)-(5) will generate the
corresponding TH fields with the following states

$\displaystyle|\phi^{\prime}_{1}\rangle$
$\displaystyle=2^{-1/2}\left(|V\rangle|\\!+\\!3\rangle+|H\rangle|\\!-\\!3\rangle\right),$
(6) $\displaystyle|\phi^{\prime}_{2}\rangle$
$\displaystyle=2^{-1/2}\left(|V\rangle|\\!+\\!12\rangle-|H\rangle|\\!-\\!12\rangle\right),$
(7) $\displaystyle|\phi^{\prime}_{3}\rangle$
$\displaystyle=1.721(\sin^{3}35^{\circ}|V\rangle|\\!+\\!3\rangle+\cos^{3}35^{\circ}|H\rangle|\\!-\\!3\rangle).$
(8)

To evaluate the qualities of the prepared FW vector fields and the generated
TH vector fields, we will use the Stokes parameters normalized by the total
intensity, as $S_{1}\\!=\\!E_{x}E^{\ast}_{x}\\!-\\!E_{y}E^{\ast}_{y}$,
$S_{2}\\!=\\!E_{x}E^{\ast}_{y}\\!+\\!E^{\ast}_{x}E_{y}$ and
$S_{3}\\!=\\!j(E_{x}E^{\ast}_{y}\\!-\\!E^{\ast}_{x}E_{y})$ Born1999 ; Wang2010
, where $E_{x}$ and $E_{y}$ indicate the normalized complex amplitudes of the
$|H\rangle$\- and $|V\rangle$-polarized components, respectively. We measure
the total intensity patterns and the intensity patterns under the polarization
bases of $|H\rangle$, $|V\rangle$,
$|D\rangle\\!=\\!(|H\rangle-|V\rangle)/\sqrt{2}$,
$|A\rangle\\!=\\!(|H\rangle+|V\rangle)/\sqrt{2}$, $|R\rangle$ and $|L\rangle$.
The second-fourth rows of Fig. 3 show the normalized $S_{1}$, $S_{2}$ and
$S_{3}$ calculated from the corresponding patterns measured above. Clearly,
all the experimental results of three prepared FW vector fields and their TH
vector fields are in good agreement with our expectations in Eqs. (3)-(8).

Figure 3: Experimentally measured total intensity patterns and the calculated
normalized Stokes parameters with values ranging from $-1$ to $+1$ by the
difference of the measured intensity patterns for the prepared three FW vector
fields (first, third and fifth columns) and their TH vector fields (second,
fourth and sixth columns).

As for the conversion efficiency of our QPPKTP, when the FW field is a
fundamental Gaussian mode with an average power of 2 W and a radius of 100
$\mu$m, the SHG efficiency is 17.6$\%$ and the THG efficiency is 8.0$\%$,
suggesting that the QPPKTP we designed and prepared should be efficient.
However, the power of the generated FW vector fields is limited by the damage
threshold of the SLM in the near infrared band and the FW vector fields have
an average pump power of 91.2 mW in our experiment. In addition, to achieve
the high-quality TH vector fields, the FW vector fields are weakly focused
into the QPPKTP, resulting in the lower power density, hence the THG
efficiency are relatively low to be 1.2$\times$10-4, 3.0$\times$10-5,
1.0$\times$10-4 and 3.1$\times$10-5 for the FW fields of
$|V\rangle|\\!-\\!1\rangle$, $|L\rangle|\\!-\\!2\rangle$,
$|H\rangle|\\!-\\!3\rangle$ and $|\phi_{1}\rangle$ in Eq. (3), respectively.
The THG efficiencies of the two FW fields of $|L\rangle|\\!-\\!2\rangle$ and
$|\phi_{1}\rangle$ are $\sim$1/4 that of the other two FW fields, as the
reasons stated below. Differently from the THGs of $|V\rangle|\\!-\\!1\rangle$
and $|H\rangle|\\!-\\!3\rangle$ in only one loop, the two FW
$|L\rangle|\\!-\\!2\rangle$ and $|\phi_{1}\rangle$ have to be divided into the
two loops with a half of power to realize the THGs, resulting in that the THG
efficiency drops to $\sim$1/8 in one loop and the sum of THG efficiency in the
two loops becomes $\sim$1/4. In future, the THG efficiency could be increased
by increasing the pump power density through preparing high-power FW vector
field and more tightly focusing. If the dual-wavelength elements are replaced
by triple-wavelength elements, the experimental setup can output both the SH
wave and TH wave of spatially structured light simultaneously.

In summary, we successfully demonstrated the THG of spatially structured light
with only one QPPKTP, and in future, generalizing nonlinear crystal design to
3 dimensions may enable to generate SH and TH spatially structured light by
using 3D nonlinear crystals Zhang2021  instead of shaping the FW with an SLM.
By using a nonlinear Sagnac interferometer, we realized the THG of vector
fields with space-varying polarization states for the first time. Our approach
may find potential applications to investigate optical skyrmions Gao2020 ;
Karnieli2021  by exploring pseudo-spin textures of vector fields in nonlinear
optics. The inverse process, SPDC, is expected to produce spatially structured
entangled photons and may open up new opportunity for high-dimension Mair2001
; Dada2011 ; Erhard2020  quantum information processing.

Funding. National Natural Science Foundation of China (No. 11922406 and No.
91750202); National Key R&D Program of China (No. 2019YFA0308700, No.
2020YFA0309500, and No. 2018YFA0306200); Key R&D Program of Guangdong Province
(No. 2020B0303010001).

Acknowledgment. The authors would like to thank the support by Collaborative
Innovation Center of Extreme Optics.

## References

  * (1) H. Rubinsztein-Dunlop, A. Forbes, M. V. Berry, M. R. Dennis, D. L. Andrews, M. Mansuripur, C. Denz, C. Alpmann, P. Banzer, T. Bauer, E. Karimi, L. Marrucci, M. J. Padgett, M. Ritsch-Marte, N. M. Litchinitser, N. P. Bigelow, C. Rosales-Guzmán, A. Belmonte, J. P. Torres, T. W. Neely, M. Baker, R. Gordon, A. B. Stilgoe, J. Romero, A. G. White, R. Fickler, A. E. Willner, G. Xie, B. J. McMorran, and A. M. Weiner, J. Opt. 19, 013001 (2017).
  * (2) J. Wang, F. Castellucci, and S. Franke-Arnold, AVS Quantum Sci. 2, 031702 (2020).
  * (3) A. Forbes, M. Oliveira, and M. R. Dennis, Nature Photon. 15, 253 (2021).
  * (4) D. Giovannini, J. Romero, V. Potocek, G. Ferenczi, F. Speirits, S. M. Barnett, D. Faccio, and M. J. Padgett, Science 347, 857 (2015).
  * (5) H. Wang, L. Shi, B. Lukyanchuk, C. Sheppard, and C. T. Chong, Nature Photon. 2, 501 (2008).
  * (6) T. Bauer, P. Banzer, E. Karimi, S. Orlov, A. Rubano, L. Marrucci, E. Santamato, R. W. Boyd, and G. Leuchs, Science 347, 964 (2015).
  * (7) A. Chong, C. Wan, J. Chen, and Q. Zhan, Nature Photon. 14, 350 (2020).
  * (8) D. L. Andrews, Structured light and its applications: an introduction to phase-structured beams and nanoscale optical forces (Academic Press-Elsevier, Burlington, 2008).
  * (9) A. Forbes and I. Nape, AVS Quantum Sci. 1, 011701 (2019).
  * (10) L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, Phys. Rev. A. 45, 8185 (1992).
  * (11) M. J. Padgett, Opt. Express 25, 11265 (2017).
  * (12) K. Dholakia, N. B. Simpson, M. J. Padgett, and L. Allen, Phys. Rev. A 54, R3742 (1996).
  * (13) N. Voloch-Bloch, K. Shemer, A. Shapira, R. Shiloh, I. Juwiler, and A. Arie, Phys. Rev. Lett. 108, 233902 (2012).
  * (14) K. Shemer, N. Voloch-Bloch, A. Shapira, A. Libster, I. Juwiler, and A. Arie, Opt. Lett. 38, 5470 (2013).
  * (15) S. M Li, L. J Kong, Z. C. Ren, Y. N. Li, C. H. Tu, and H. T. Wang, Phys. Rev. A. 88, 035801 (2013).
  * (16) Z. Y. Zhou, Y. Li, D.-S. Ding, W. Zhang, S. Shi, B. S. Shi, and G. C. Guo, Opt. Express 22, 23673 (2014).
  * (17) N. A. Chaitanya, A. Aadhi, M. V. Jabir, and G. K. Samanta, Opt. Lett. 40, 2614 (2015).
  * (18) R. Ni, Y. F. Niu, L. Du, X. P. Hu, Y. Zhang, and S. N. Zhu, Appl. Phys. Lett. 109, 151103 (2016).
  * (19) Y. Chen, R. Ni, Y. Wu, L. Du, X. Hu, D. Wei, Y. Zhang, and S. N. Zhu, Phys. Rev. Lett. 125, 143901 (2020).
  * (20) Y. Tang, K. Li, X. Zhang, J. Deng, G. Li, and E. Brasselet, Nature Photon. 14 658 (2020).
  * (21) G. Gui, N. J. Brooks, H. C. Kapteyn, M. M. Murnane, and C.-T. Liao, Nature Photon. 14, 658 (2021).
  * (22) S. W. Hancock, S. Zahedpour, and H. M. Milchberg, Optica 8, 594 (2021).
  * (23) X. Y. Fang, G. Yang, D. Z. Wei, D. Wei, R. Ni, W. Ji, Y. Zhang, X. P. Hu, W. Hu, Y. Q. Lu, S. N. Zhu, and M. Xiao, Opt. Lett. 41, 1169 (2016).
  * (24) Z. Xu, Z. Y. Lin, Z. L. Ye, Y. Chen, X. P. Hu, Y. D. Wu, Y. Zhang, P. Chen, W. Hu, Y. Q. Lu, M. Xiao, S. N. Zhu, Opt. Express 26, 17563 (2018).
  * (25) D. Gauthier, P. R. Ribic, G. Adhikary,A. Camper,C. Chappuis, R. Cucini, L. F. DiMauro,G. Dovillaire,F. Frassetto, R. Geneaux, P. Miotti ,L. Poletto, B. Ressel, C. Spezzani, M. Stupar, T. Ruchon, and G. De Ninno, Nature Commun. 8, 14971 (2017).
  * (26) A. Mair, A. Vaziri, G. Weihs, and A. Zeilinger, Nature 412, 313 (2001).
  * (27) J. Leach, B. Jack, J. Romero, A. M. Jha, A. M. Yao, S. Franke-Arnold, D. G. Ireland, R. W. Boyd, S. M. Barnett, and M. J. Padgett, Science 329, 662 (2010).
  * (28) A. C. Dada, J. Leach,G. S. Buller, M. J. Padgett, and E. Andersson, Nat. Phys. 7, 677 (2011).
  * (29) X. Qiu, F. Li, W. Zhang, Z. Zhu, and L. Chen, Optica 5, 208 (2018).
  * (30) X. Y. Fang, H. C. Yang, W. Z. Yao, T. X. Wang, Y. Zhang, M. Gu, and M. Xiao, Adv. Photon. 3, 015001 (2021).
  * (31) M. Erhard, M. Krenn, and A. Zeilinger, Nature Rev. Phys. 2, 365 (2020).
  * (32) F. Cardano and L. Marrucci, Nature Photon. 9, 776 (2015).
  * (33) L. Aolita and S. P. Walborn, Phys. Rev. Lett. 98, 100501 (2007).
  * (34) V. D’Ambrosio, E. Nagali, S. P. Walborn, L. Aolita, S. Slussarenko, L. Marrucci, and F. Sciarrino, Nature Commun. 3, 961 (2012).
  * (35) J. T. Barreiro, T.-C Wei, and P. G. Kwiat, Nature Phys. 4, 282 (2008).
  * (36) Z. Y. Zhu, M. Janasik, A. Fyffe, D. Hay, Y. Y. Zhou, B. Kantor, T. Winder, R. W. Boyd, G. Leuchs, and Z. M. Shi, Nature Commun. 12, 1666 (2021).
  * (37) C. Yang, Z. Y. Zhou, Y. Li, Y. H. Li, S. L. Liu, S. K. Liu, Z. H. Xu, G. C. Guo, and B. S. Shi, Opt. Lett. 44, 219 (2019).
  * (38) H. G. Liu, H. Li, Y. L. Zheng, and X. F. Chen, Opt. Lett. 43, 5981 (2018).
  * (39) S. N. Zhu, Y. Y. Zhu, and N. B. Ming, Science 278, 843 (1997).
  * (40) C. Zhang, H. Wei, Y. Y. Zhu, H. T. Wang, S. N. Zhu and N. B. Ming, Opt. Lett. 26, 899 (2001).
  * (41) A. Vernay, L. Bonnet-Gamard, V. Boutou, S. Trajtenberg-Mills, A. Arie, and B. Boulanger, OSA Continuum. 3, 1536 (2020).
  * (42) R. Lifshitz, A. Arie, and A. Bahabad, Phys. Rev. Lett. 95, 133901 (2005).
  * (43) R. K. P. Zia and W. J. Dallas, J. Phys. A 18, L341 (1985).
  * (44) Z. C. Ren, Z. M. Cheng, X. L. Wang, J. P. Ding and H. T. Wang, Appl. Phys. Lett. 118, 011105 (2021).
  * (45) P. Vaity, J. Banerji, and R. P. Singh, Phys. Lett. A 377, 1154 (2013).
  * (46) M. Born and E. Wolf, Principles of Optics, 7th ed. (Cambridge U. Press, 1999).
  * (47) X. L. Wang, Y. Li, J. Chen, C. S. Guo, J. Ding, and H. T. Wang, Opt. Express 18, 10786 (2010).
  * (48) Y. Zhang, Y. Sheng, S. N. Zhu, M. Xiao, and W. Krolikowski, Optica 8, 372 (2021).
  * (49) S. J Gao, F. C. Speirits, F. Castellucci, S. Franke-Arnold, S. M. Barnett, and J. B. Götte, Phys. Rev. A 102, 053513(2020).
  * (50) A. Karnieli, S. Tsesses, G. Bartal, and A. Arie, Nature Commun. 12, 1092 (2021).