# Impact of Variable Photospheric Radius on Exoplanet Atmospheric Retrievals

Jake Taylor1
1Institute for Research on Exoplanets, Department of Physics, University of
Montréal, Montréal, H2V 0B3, Canada E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Inverse techniques are used to extract information about an exoplanet’s
atmosphere. These techniques are prone to biased results if the appropriate
forward model is not used. One assumption used in a forward model is to assume
that the radius of the planet is constant with wavelength, however a more
realistic assumption is that the photospheric radius varies with each
wavelength. We explore the bias induced when attempting to extract the
molecular abundance from an emission spectrum which was generated with a
variable radius. We find that for low gravity planets, the retrieval model is
not able to fit the data if a constant radius model is used. We find that
biased results are obtained when studying a typical hot Jupiter in the MIRI
LRS wavelength range. Finally, we show that high gravity planets do not suffer
a bias. We recommend that future spectral retrievals that interpret exoplanet
emission spectra should take into account a variable radius.

###### keywords:

methods: data analysis –techniques: spectroscopic

††pubyear: 2022††pagerange: Impact of Variable Photospheric Radius on
Exoplanet Atmospheric Retrievals–Impact of Variable Photospheric Radius on
Exoplanet Atmospheric Retrievals

##  1 Introduction

Studying the emission spectra of exoplanets using spectral retrievals has been
the pioneering technique over the past decade to infer planetary properties
(Lee et al., 2012; Line et al., 2014; Stevenson et al., 2014; Line et al.,
2016; Evans et al., 2017; Edwards et al., 2020; Changeat, 2021). Current work
has utilised both the Hubble Space Telescope (HST) and the Spitzer Space
Telescope. With such limited spectral coverage, assumptions about the forward
model used in the retrieval framework need to be made. These assumptions
include, the parameterization of the temperature-pressure structure, the
molecules in the atmosphere, and the geometry and the cloud composition.

With the successful launch of the James Webb Space Telescope (JWST) on the
25th of December 2021 the exoplanet community will have access to spectral
observations of exoplanets from $0.6-12$ microns at unprecedented spectral
resolution. The higher signal to noise observations will require more complex
forward models such as geometries that extend past the current 1D assumption
(Feng et al., 2016; Taylor et al., 2020; Feng et al., 2020; Irwin et al.,
2020; Changeat & Al-Refaie, 2020) and the parameterisation of clouds (Taylor
et al., 2021).

A common model assumption is that the radius of the planet is constant with
wavelength. Fortney et al. (2019) compared the effect of computing the radius
at each wavelength has on the emission spectra. They found that while the
effect is not important for current HST observations, it can cause a $10-20$ %
difference in the wavelength dependent eclipse depth for JWST quality data.
Therefore in this study we explore how the variable radius assumption can bias
our retrieved results when studying JWST quality observations.

##  2 Methods

In this section we outline the model set up and retrieval framework used.

The outline of the methodology of this study are as follows:

  * •

Generate a model of a hot Jupiter which has a radius that varies with
wavelength

  * •

Consider 3 different observing scenarios: NIRSpec PRISM, MIRI LRS and NIRSpec
PRISM + MIRI LRS

  * •

For each of the observing scenarios, considering 3 different error envelopes:
100ppm, 60ppm and 30ppm.

  * •

Retrieve on these simulated datasets with a retrieval set up which assumes
that the radius is fixed at each wavelength

###  2.1 Model Setup

Figure 1: Left: Forward models used in this study with variable radius shown
by the dashed line and constant radius with the solid line. In purple we
present the model based on WASP-43b. In red we present the low gravity
scenario. In blue we present the high gravity scenario. Right: The residuals
when taking the difference between the variable and constant radius models.

It was found by Fortney et al. (2019) that the impact of the variable radius
is dependent on the surface gravity of the planet. Hence, we explore 3
atmospheric scenarios:

  1. 1.

A cloud free model based on the physical properties of WASP-43b.

  2. 2.

Same as (1) however with 10% of the mass. We call this the low gravity
scenario. This is analogous to a low gravity planet such as WASP-107b.

  3. 3.

Same as (1) however with 2$\times$ the mass. We call this the high gravity
scenario. This is analogous to a high gravity planet such as WASP-18b.

These scenarios let us explore a range of planetary scenarios observable with
JWST. To generate these model spectra we use the CHIMERA code (see Line et al.
(2013) for more detailed explanation of the capabilities of the framework, the
open source code can be found here: https://github.com/mrline/CHIMERA), we
show the main system parameters used in this study in Table 1. We use a
parameterised temperature-pressure profile developed by Parmentier & Guillot
(2014) which balances incoming shortwave radiation with the outgoing longwave
radiation. For simplicity we have used a blackbody spectrum calculated at T =
4400K for the stellar model. We present the models in Figure 1, it can see
seen that when comparing the difference between a model with a variable radius
and one with a constant radius, that there is a difference which increases
with wavelength. We see that the difference is larger for lower gravity cases,
which is consistent with Fortney et al. (2019). We note that our study only
considers hot Jupiters. Fortney et al. (2019) points out that the impact may
be different for ultra hot Jupiters.

Parameter | Value | Parameter | Value
---|---|---|---
log(H2O)1 | -3 | log(CO2)4 | -8
log(CO)2 | -3 | log(CH4)3 | -8
Teq [K] | 1400 | Rs [R${}_{\text{sun}}$]5 | 0.667
log($\kappa_{IR}$) | -1.5 | Tstar [K] | 4400
log($\gamma$) | -0.7 | SMA [AU]5 | 0.015
Mp [Mj]5 | 2.034 | Rp [Rj]5 | 1.036
Table 1: Atmospheric and planetary properties used to create the model
atmospheres. 1Polyansky et al. (2018), 2Li et al. (2015), 3Yurchenko &
Tennyson (2014), 4Yurchenko et al. (2020), 5Gillon et al. (2012).

###  2.2 Variable Radius Implementation

We upgraded the CHIMERA code to be able to calculate the flux when considering
a change in the photospheric radius (defined as the radius at which optical
depth $\tau(\lambda)=2/3$) at each wavelength.

The planet-to-star flux is given by the equation:

$F=\Big{(}\frac{R_{p}(\lambda)}{R_{s}(\lambda)}\Big{)}^{2}\frac{F_{p}(\lambda)d\lambda}{F_{s}(\lambda)d\lambda}$
(1)

where $R_{p}(\lambda)$ is the planetary radius, which is shown to have a
wavelength dependence. The planetary radius is often assumed to be constant
with wavelength.

To account for the wavelength dependence we rewrite Equation 1 to be:

$F=\Big{(}\frac{R_{p}+dR_{p}(\lambda)}{R_{s}(\lambda)}\Big{)}^{2}\frac{F_{p}(\lambda)d\lambda}{F_{s}(\lambda)d\lambda}$
(2)

where $dR_{p}(\lambda)$ is the change in radius at each wavelength. We can
compute this $dR_{p}(\lambda)$ by first extracting the altitude levels in the
atmosphere ($Z(\lambda)$) and dividing this by the planetary radius ($R_{0}$),
$Z(\lambda)/R_{0}$. By calculating the cumulative optical depth
($\tau(\lambda)$) as a function of $Z(\lambda)$, it is possible to determine
the altitude at $\tau(\lambda)=2/3$, hence this would provide the
$dR_{p}(\lambda)$ at each wavelength. We note that in this study we are
considering this effect in 1D and the $\tau(\lambda)=2/3$ should vary for
different emission angles.

###  2.3 Sampling Method

To perform the retrievals we used the dynesty package (Speagle, 2020) which
implements a nested sampling algorithm (Skilling, 2004, 2006). This technique
has been effective in extracting atmospheric information from spectra (Rathcke
et al., 2021; Lustig-Yaeger et al., 2021; Ahrer et al., 2022).

##  3 Results

In this section we will outline the results from the three different planetary
scenarios. We retrieve on the variable radius models shown in Figure 1 using a
fixed radius assumption in our retrieval model. We perform this for three
different wavelength ranges obtainable with JWST. These are: NIRSpec PRISM,
MIRI LRS and then combining the two instrument modes. For each of these
wavelength ranges a 100ppm, 60ppm and 30ppm error envelope is assumed to
assess the limits in which a bias may arise.

###  3.1 WASP-43b

Figure 1 shows that the difference between the variable and constant radius
models do not deviate until mid-IR, with around 100ppm difference at 4.5
microns. Figure 2 shows the best fitting spectra for all the observing modes
for the 100ppm case. Table 2 shows the retrieved abundances for each modelling
set up, accompanied by the reduced $\chi^{2}$.

Figure 2: The best fitting model and simulated data for WASP-43b. The top,
middle, and bottom panel show the NIRSpec observing mode, the NIRSpec + MIRI
LRS observing mode, the MIRI LRS observing mode, respectively. The reduced
$\chi^{2}$ of the fit is shown inset.

Observations with NIRSpec PRISM show that the variable radius does not bias
the retrieved chemical abundances. For both the 100ppm and 60ppm cases, the
abundance of CO2 is unconstrained, however at the 30ppm level it is possible
to constrain the abundance of CO2. For observations with MIRI LRS, the
abundance of H2O is retrieved to the correct value within 2$\sigma$ for the
100ppm and 60ppm case. A biased result is retrieved for the 30ppm case. When
combining instrument modes H2O and CO are retrieved to be greater than
2$\sigma$ away from the input value. It is interesting that the combining
NIRSpec PRISM with MIRI LRS results in a biased abundance of CO, as there is
no CO information in this wavelength region. When combining instrument modes
we begin to detect CO2 at the 60ppm level and then constrain the abundance to
the correct value at the 30ppm level.

###  3.2 Low Gravity Scenario

Figure 1 shows that the difference between the variable and constant radius
models is large throughout the whole infrared. Although the initial 100ppm
simulation produced the incorrect chemistry, this is not the interesting
takeaway point. The constant radius model is not able to fit the variable
radius simulated data, we show this in Figure 3. Therefore low gravity planets
will require a variable radius model to be able to fit the observations. After
17 days the MIRI LRS still had not converged, compared to 0.477 days and 0.282
days from the WASP-43b and High Gravity cases. Each simulation was run on 12
central processing units (CPU) cores. The forward model, which assumes a
constant radius, is too simple. The sampling algorithm cannot converge to
produce an adequate reduced $\chi^{2}$. As a result of the biases being
observed at the 100ppm level, we do not explore the 60ppm or 30ppm level.

Figure 3: The best fitting model and simulated data for the low gravity case.
The top and bottom panel show the NIRSpec observing mode and the NIRSpec +
MIRI LRS observing mode, respectively. The reduced $\chi^{2}$ of the fit is
shown inset.

###  3.3 High Gravity Scenario

Figure 1 shows that the difference in the models is only around 50ppm at 4.5
microns. Figure 4 shows the best fitting models from the 100ppm case and
present the retrieved abundances in Table 4. For each of the NIRSpec PRISM and
NIRSpec PRISM + MIRI LRS observing modes, the abundance of H2O and CO were
constrained within 2$\sigma$ of the input value. For the MIRI LRS observing
mode the H2O abundance was constrained within 2$\sigma$ and the rest of the
chemistry is unconstrained, which was to be expected. We can summarise that a
constant radius model is sufficient to interpret the atmosphere of a high
gravity exoplanet.

Figure 4: The best fitting model and simulated data for the high gravity case.
The top, middle, and bottom panel show the NIRSpec observing mode, the NIRSpec
+ MIRI LRS observing mode, the MIRI LRS observing mode, respectively. The
reduced $\chi^{2}$ of the fit is shown inset.

##  4 Conclusion

In this study we explored the impact of a variable photospheric radius on a
retrieval analysis of emission spectra. We explored three model scenarios
which span a range of gravities. Our conclusions can be summarised as follows:

  1. 1.

For a typical hot Jupiter, the variable radius assumption is not important
when interpreting observations in the NIRSpec PRISM wavelength range. However,
the impact of a variable radius becomes important when observing over the MIRI
LRS wavelength region. If the variable radius is not assumed in the retrieval
model then biased chemistry will be found.

  2. 2.

For a low gravity hot Jupiter it is not possible to fit the emission spectrum
with a constant radius model. Hence for all observing scenarios of a low
gravity hot Jupiter, a variable radius model needs to be used for accurate
interpretation.

  3. 3.

For a high gravity hot Jupiter no biases were found when retrieving on the
variable radius model with a constant radius model.

  4. 4.

It is possible to constrain the abundance of CO2 unambiguously in the WASP-43b
case with an error envelope of 30ppm. A similar constraint is placed on the
abundance of CO2 in the high gravity case. The upper limit is well
constrained, however a large tail is seen in the lower limit.

From this study we recommend that the community include a variable radius
model in their retrieval frameworks as it will be important in the JWST era.
This can be implemented in a similar way as discussed in Section 2.2. Given
that the optical depths are already computed in the radiative transfer
calculation, determining and then adding $dR_{p}(\lambda)$ does not increase
computation time.

##  5 Acknowledgements

Jake Taylor thanks the Canadian Space Agency for financially supporting this
work. He thanks Dr Peter McGill for proof reading the manuscript. Jake thanks
Dr Michael Line for discussions about the photospheric variable radius and the
use of CHIMERA. Jake thanks Dr Joanna Barstow for the review of this letter
which greatly improved its clarity.

## Data Availability

Any data generated in this study is available on request.

## References

  * Ahrer et al. (2022) Ahrer E., Wheatley P. J., Kirk J., Gandhi S., King G. W., Louden T., 2022, MNRAS, 510, 4857
  * Changeat (2021) Changeat Q., 2021, arXiv e-prints, p. arXiv:2112.09973
  * Changeat & Al-Refaie (2020) Changeat Q., Al-Refaie A., 2020, ApJ, 898, 155
  * Edwards et al. (2020) Edwards B., et al., 2020, AJ, 160, 8
  * Evans et al. (2017) Evans T. M., et al., 2017, Nature, 548, 58
  * Feng et al. (2016) Feng Y. K., Line M. R., Fortney J. J., Stevenson K. B., Bean J., Kreidberg L., Parmentier V., 2016, ApJ, 829, 52
  * Feng et al. (2020) Feng Y. K., Line M. R., Fortney J. J., 2020, AJ, 160, 137
  * Fortney et al. (2019) Fortney J. J., Lupu R. E., Morley C. V., Freedman R. S., Hood C., 2019, ApJ, 880, L16
  * Gillon et al. (2012) Gillon M., et al., 2012, A&A, 542, A4
  * Irwin et al. (2020) Irwin P. G. J., Parmentier V., Taylor J., Barstow J., Aigrain S., Lee E. K. H., Garland R., 2020, MNRAS, 493, 106
  * Lee et al. (2012) Lee J. M., Fletcher L. N., Irwin P. G. J., 2012, MNRAS, 420, 170
  * Li et al. (2015) Li G., Gordon I. E., Rothman L. S., Tan Y., Hu S.-M., Kassi S., Campargue A., Medvedev E. S., 2015, The Astrophysical Journal Supplement Series, 216, 15
  * Line et al. (2013) Line M. R., et al., 2013, ApJ, 775, 137
  * Line et al. (2014) Line M. R., Knutson H., Wolf A. S., Yung Y. L., 2014, ApJ, 783, 70
  * Line et al. (2016) Line M. R., et al., 2016, AJ, 152, 203
  * Lustig-Yaeger et al. (2021) Lustig-Yaeger J., Stevenson K. B., Mayorga L. C., Showalter Sotzen K., May E. M., Izenberg N. R., Mandt K., 2021, ApJ, 921, L4
  * Parmentier & Guillot (2014) Parmentier V., Guillot T., 2014, Astronomy & Astrophysics, 562, A133
  * Polyansky et al. (2018) Polyansky O. L., Kyuberis A. A., Zobov N. F., Tennyson J., Yurchenko S. N., Lodi L., 2018, Monthly Notices of the Royal Astronomical Society, 480, 2597
  * Rathcke et al. (2021) Rathcke A. D., et al., 2021, AJ, 162, 138
  * Skilling (2004) Skilling J., 2004, in Fischer R., Preuss R., Toussaint U. V., eds, American Institute of Physics Conference Series Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering. pp 395–405, doi:10.1063/1.1835238
  * Skilling (2006) Skilling J., 2006, Bayesian Analysis, 1, 833
  * Speagle (2020) Speagle J. S., 2020, MNRAS, 493, 3132
  * Stevenson et al. (2014) Stevenson K. B., Bean J. L., Madhusudhan N., Harrington J., 2014, ApJ, 791, 36
  * Taylor et al. (2020) Taylor J., Parmentier V., Irwin P. G. J., Aigrain S., Lee E. K. H., Krissansen-Totton J., 2020, MNRAS, 493, 4342
  * Taylor et al. (2021) Taylor J., Parmentier V., Line M. R., Lee E. K. H., Irwin P. G. J., Aigrain S., 2021, MNRAS, 506, 1309
  * Yurchenko & Tennyson (2014) Yurchenko S. N., Tennyson J., 2014, Monthly Notices of the Royal Astronomical Society, 440, 1649
  * Yurchenko et al. (2020) Yurchenko S. N., Mellor T. M., Freedman R. S., Tennyson J., 2020, MNRAS, 496, 5282

NIRSpec PRISM |  |  |  |  |
---|---|---|---|---|---
Model | logH2O (-3) | logCO (-3) | logCO2 (-8) | logCH4 (-8) |  $\chi^{2}$/n
100ppm | -2.99${}^{+0.40}_{-0.24}$ | -3.05${}^{+0.63}_{-0.43}$ | -9.67${}^{+2.10}_{-2.20}$ | -9.35${}^{+2.08}_{-1.99}$ | 0.18
60ppm | -3.02${}^{+0.19}_{-0.14}$ | -3.10${}^{+0.32}_{-0.26}$ | -9.68${}^{+1.91}_{-2.21}$ | -9.46${}^{+2.56}_{-2.42}$ | 0.51
30ppm | -3.03${}^{+0.09}_{-0.07}$ | -3.12${}^{+0.14}_{-0.13}$ | -8.74${}^{+0.85}_{-3.10}$ | -9.59${}^{+2.35}_{-2.29}$ | 1.89
MIRI LRS |  |  |  |  |
Model | logH2O (-3) | logCO (-3) | logCO2 (-8) | logCH4 (-8) |  $\chi^{2}$/n
100ppm | -2.24${}^{+1.01}_{-1.56}$ | -6.65${}^{+4.80}_{-5.05}$ | -8.25${}^{+3.49}_{-3.53}$ | -8.68${}^{+3.16}_{-3.10}$ | 0.15
60ppm | -2.24${}^{+0.80}_{-0.99}$ | -6.75${}^{+4.90}_{-4.88}$ | -8.60${}^{+3.67}_{-3.23}$ | -8.81${}^{+2.97}_{-3.01}$ | 0.39
30ppm | -2.21${}^{+0.37}_{-0.38}$ | -6.75${}^{+4.76}_{-4.92}$ | -8.61${}^{+3.16}_{-3.08}$ | -9.09${}^{+2.87}_{-2.69}$ | 1.58
NIRSpec PRISM + MIRI LRS |  |  |  |  |
Model | logH2O (-3) | logCO (-3) | logCO2 (-8) | logCH4 (-8) |  $\chi^{2}$/n
100ppm | -2.68${}^{+0.38}_{-0.33}$ | -2.52${}^{+0.44}_{-0.42}$ | -9.20${}^{+2.05}_{-2.65}$ | -9.03${}^{+2.67}_{-2.81}$ | 1.25
60ppm | -2.68${}^{+0.22}_{-0.21}$ | -2.53${}^{+0.26}_{-0.27}$ | -8.29${}^{+1.03}_{-3.51}$ | -9.19${}^{+2.61}_{-2.67}$ | 3.43
30ppm | -2.68${}^{+0.11}_{-0.11}$ | -2.54${}^{+0.14}_{-0.14}$ | -7.60${}^{+0.22}_{-0.41}$ | -9.13${}^{+2.41}_{-2.71}$ | 12.24
Table 2: The retrieved chemical abundances for the WASP-43b model. We show each observing mode and error envelope. We also present the reduced $\chi^{2}$ where n = 7. The error ranges are the 2$\sigma$ ranges as outputted from dynesty. We present the true model values in brackets next to the parameters name. We have colour coordinated the results so they are easier for the reader to interpret. We present results within 2$\sigma$ in purple, outside of 2$\sigma$ in red, results with an upper limit inside 2$\sigma$ in blue and unconstrained results in black. NIRSpec PRISM |  |  |  |  |
---|---|---|---|---|---
Model | logH2O (-3) | logCO (-3) | logCO2 (-8) | logCH4 (-8) |  $\chi^{2}$/n
100ppm | -2.00${}^{+0.18}_{-0.25}$ | -8.16${}^{+3.81}_{-3.66}$ | -9.93${}^{+2.28}_{-1.96}$ | -8.77${}^{+3.18}_{-3.05}$ | 81.96
NIRSpec PRISM + MIRI LRS |  |  |  |  |
Model | logH2O (-3) | logCO (-3) | logCO2 (-8) | logCH4 (-8) |  $\chi^{2}$/n
100ppm | -9.26${}^{+2.65}_{-2.57}$ | -1.76${}^{+0.55}_{-2.14}$ | -7.03${}^{+1.54}_{-4.61}$ | -1.83${}^{+0.28}_{-0.27}$ | 865.62
Table 3: The retrieved chemical abundances for the Low Gravity model. We show each observing mode and error envelope. We also present the reduced $\chi^{2}$ where n = 7. The error ranges are the 2$\sigma$ ranges as outputted from dynesty.We present the true model values in brackets next to the parameters name. We have colour coordinated the results so they are easier for the reader to interpret. We present results within 2$\sigma$ in purple, outside of 2$\sigma$ in red, results with an upper limit inside 2$\sigma$ in blue and unconstrained results in black. NIRSpec PRISM |  |  |  |  |
---|---|---|---|---|---
Model | logH2O (-3) | logCO (-3) | logCO2 (-8) | logCH4 (-8) |  $\chi^{2}$/n
100ppm | -2.99${}^{+0.30}_{-0.20}$ | -3.01${}^{+0.50}_{-0.42}$ | -9.61${}^{+2.05}_{-2.28}$ | -9.29${}^{+2.75}_{-2.55}$ | 0.06
60ppm | -3.01${}^{+0.15}_{-0.13}$ | -3.04${}^{+0.28}_{-0.25}$ | -9.62${}^{+1.94}_{-2.25}$ | -9.37${}^{+2.53}_{-2.50}$ | 0.15
30ppm | -3.01${}^{+0.07}_{-0.07}$ | -3.05${}^{+0.13}_{-0.13}$ | -9.01${}^{+1.18}_{-2.82}$ | -9.54${}^{+2.41}_{-2.35}$ | 0.50
MIRI LRS |  |  |  |  |
Model | logH2O (-3) | logCO (-3) | logCO2 (-8) | logCH4 (-8) |  $\chi^{2}$/n
100ppm | -2.47${}^{+1.26}_{-1.40}$ | -6.79${}^{+5.09}_{-4.90}$ | -8.40${}^{+3.70}_{-3.39}$ | -8.74${}^{+3.28}_{-3.08}$ | 0.03
60ppm | -2.50${}^{+0.95}_{-1.21}$ | -6.84${}^{+4.89}_{-4.73}$ | -8.19${}^{+3.31}_{-3.50}$ | -8.84${}^{+3.06}_{-2.92}$ | 0.10
30ppm | -2.49${}^{+0.50}_{-0.69}$ | -6.84${}^{+4.81}_{-4.80}$ | -8.64${}^{+3.28}_{-3.07}$ | -9.10${}^{+3.02}_{-2.72}$ | 0.39
NIRSpec PRISM + MIRI LRS |  |  |  |  |
Model | logH2O (-3) | logCO (-3) | logCO2 (-8) | logCH4 (-8) |  $\chi^{2}$/n
100ppm | -2.93${}^{+0.30}_{-0.22}$ | -2.87${}^{+0.39}_{-0.35}$ | -9.59${}^{+2.20}_{-2.29}$ | -9.06${}^{+2.61}_{-2.81}$ | 0.31
60ppm | -2.95${}^{+0.15}_{-0.14}$ | -2.89${}^{+0.22}_{-0.21}$ | -9.39${}^{+1.81}_{-2.48}$ | -9.14${}^{+2.53}_{-2.71}$ | 0.91
30ppm | -2.96${}^{+0.07}_{-0.07}$ | -2.90${}^{+0.10}_{-0.10}$ | -8.11${}^{+0.44}_{-3.57}$ | -9.07${}^{+2.30}_{-2.79}$ | 3.29
Table 4: The retrieved chemical abundances for the High Gravity model. We show
each observing mode and error envelope. We also present the reduced $\chi^{2}$
where n = 7. The error ranges are the 2$\sigma$ ranges as outputted from
dynesty. We present the true model values in brackets next to the parameters
name. We have colour coordinated the results so they are easier for the reader
to interpret. We present results within 2$\sigma$ in purple, outside of
2$\sigma$ in red, results with an upper limit inside 2$\sigma$ in blue and
unconstrained results in black.


# Policy Gradient Approach to Compilation of Variational Quantum Circuits

D. A. Herrera-Martí  Université Grenoble Alpes, CEA List, 38000 Grenoble,
France

###### Abstract

We propose a method for finding approximate compilations of quantum circuits,
based on techniques from policy gradient reinforcement learning. The choice of
a stochastic policy allows us to rephrase the optimization problem in terms of
probability distributions, rather than variational parameters. This implies
that searching for the optimal configuration is done by optimizing over the
distribution parameters, rather than over the circuit free angles. The upshot
of this is that we can always compute a gradient, provided that the policy is
differentiable. We show numerically that this approach is more competitive
than those using gradient-free methods, even in the presence of depolarizing
noise, and argue analytically why this is the case. Another interesting
feature of this approach to variational compilation is that it does not need a
separate register and long-range interactions to estimate the end-point
fidelity. We expect these techniques to be relevant for training variational
circuit in different contexts.

##  I Introduction

The general problem of quantum compilation is to approximate any unitary
transform as a sequence of elements selected from a fixed set of universal
quantum gates. The existence of an approximate sequence of quantum gates for a
single qubit is guaranteed by the Solovay-Kitaev theorem [1]. It states that
it is possible to approximate any circuit with an overhead logarithmic in the
original number of gates, i.e. polylogarithmic as $O(\log^{c}(1/\epsilon))$,
where $\epsilon$ is the approximation accuracy and c is a constant lower-
bounded by 1 [2].

Although the Solovay-Kitaev theorem proves that any computation can be
efficiently approximated to within an arbitrary tolerance, it does not tell us
how to find the optimal sequence of gates. The standard Solovay-Kitaev
algorithm uses an exhaustive search technique to build a library of gate
sequences in its lowest level of recursion, and then builds on it recursively.
In general, the longer the sequence of library gates (and their inverses), the
better the approximation to the target unitary [3].

Finding the optimal compilation of a quantum unitary is equivalent to finding
the shortest path between two nodes (the geodesic) in a (hyperbolic) Cayley
graph [4]. Hyperbolic graphs resemble tree graphs in the sense that for an
overwhelming majority of node pairs, there is only one path linking both nodes
[5]. Therefore, the geometric intuition for the hardness of exact compilation
is that, in a hyperbolic graph, looking for a shortest path involves
evaluating an exponential number of nodes at each step. Indeed, performing an
optimal compilation of a given quantum circuit is believed to be a hard
problem under reasonable assumptions [6, 4].

With the advent of sub-threshold quantum architectures [7], research on
variational quantum algorithms has become central to the field of quantum
computing, giving rise, among others, to alternative routes to universal
quantum computation, based on variational circuits [8, 9, 10]. The existence
of variational tasks which are believed to be intractable for classical
computers [11, 12, 13, 14, 15, 16] is an encouraging motivation for research
in variational quantum algorithms. By considering a discretization of the
rotation angles, which are the free variational parameters, it can be seen
that the set of circuits that can be built using these gates has a Cayley
graph which retains its hyperbolic character [4], and therefore this method of
compilation is likely to be suboptimal for large enough unitaries. Thus we
deem our task “Approximate Compilation”.

Furthermore, several bottlenecks must be addressed before variational quantum
algorithms on a NISQ processor can be properly trained at scale. These include
the ability to evaluate efficiently cost functions which are made of non-
commuting observables[17, 18], and the design of new algorithms with no known
efficient classical simulation [19, 20, 21], in which the optimizations of the
variational parameters may suffer from the curse of dimensionality as the
number of qubits increases.

A further impediment to VQA training is due to the fact that there are regions
in parameter space, commonly known as _barren plateaus_ , where the cost
gradient vanishes exponentially in the number of qubits. This phenomenon,
which exponentially increases the resources required to train large scale
quantum neural networks, has been demonstrated in a number of proposed
architectures and classes of cost functions [22, 23]. Even gradient-free
optimizers do not solve the barren plateau problem, as cost function
differences, which are the basis for making decisions in a gradient-free
optimization, are exponentially suppressed in a barren plateau [24].

In this work, we use reinforcement learning, and in particular policy gradient
methods, to mitigate the effects of barren plateaus in the training of
variational quantum algorithms of low depth (logarithmic in the number of
qubits), and we apply it to the particular case of approximate compilation.
The intuition behind this choice comes given by the fact that in policy
gradient algorithms, the cost function can be written as an expected value of
a parameterized analytic function. It is therefore possible to classically
estimate the gradient at all times, albeit at the expense of incurring some
bias, as will be argued later.

##  II Compilation of Variational Quantum Algorithms

Our working hypothesis is that the interactions giving rise to the unknown
unitary, i.e. the qubit connectivity graph, are known. In this work, we have
assumed that they are nearest-neighbors, but more general interactions are
straighforward, provided that the interaction graph is known to the compiler.
In the absence of all information, an all-to-all circuit should be trained,
which would entail a quadratic growth of number of free parameters in the
number of qubits.

There are several approaches [8, 10, 26, 25] to assess the performance of an
approximate compilation method. We will adopt a variation of one metric
introduced in Ref.[26, 25] , motivated by its experimental feasibility:

$\hat{F}(\theta)=\frac{1}{m}\sum^{m}_{k}|\langle
k|V(\theta)^{\dagger}U|k\rangle|^{2}$ (1)

which corresponds to the fidelity between the initial and the final state. To
fully characterize $U$, a tomographically complete characterization demands
$\mathcal{O}(4^{n})$ different such initial states. Variational quantum
compiling of a full unitary matrix U by considering the action of U on a
complete basis can be computationally very demanding. A simpler task is to
learn to prepare only a particular column of the matrix U by considering the
action of U on a fixed input state. The Hilbert-Schmidt test in [26, 25] has
the advantage that it estimates the gate quantity $Tr[V(\theta)^{\dagger}U]$,
at the cost of doubling the number of qubits and introducing highly non-local
interactions. In this work, we will estimate the fidelity using only $m\sim
poly(n)$ initial states and $n$ qubits (see appendix), a fact which is denoted
by the hat punctuation in previous equation.

Approximate compilation, phrased as a variational algorithm, can get stuck in
barren plateaus [27], which have been shown to severely hinder trainability.
We will try to mitigate this problem using a stochastic version of the cost
function stated on Eq. (1), by means of optimization of a stochastic policy.

Figure 1: Circuit Setup. The goal is to retrieve a state with a maximal
overlap with the initial state $|k\rangle$, for all initial states, randomly
sampled from a fixed basis. The unknown unitary $U$ is followed by
$V(\theta)^{\dagger}$, resulting in a state $|\psi\rangle$. At the end of the
circuit, the measurement projects the state back onto the initial state with a
probability that depends on the overlap, which constitutes the reward, i.e.
$r^{(k)}_{\theta}=|\langle k|V(\theta)^{\dagger}U|k\rangle|^{2}$. The
parameterized unitary is made of single qubit rotations of the form
$R_{Y}(\theta_{i})=\exp-iY\theta_{i}$ and two-qubit rotations are of the form
$R_{ZZ}(\theta_{i})=\exp-iZZ\theta_{i}$

##  III Policy Gradient for Quantum Compilation

Policy Gradient (PG) Reinforcement Learning (RL) operates on the premise that
it is possible to optimize a parametric policy by probing the environment,
without the need of continuous update of value functions, and it constitutes
an alternative to Q-learning algorithms [28]. PGRL is naturally well suited to
handle continuous actions in stochastic environments and, provided that the
chosen policy is differentiable, the gradient of a cost function can always be
estimated. The bias incurred by this method will also be related to the
expressive power of the chosen policy (see Fig. 2). RL has found multiple
applications for quantum tasks [29, 30, 31, 32, 36]. In particular, several
works have addressed the possibility of automatically learning how to optimize
variational quantum algorithms [33, 34, 35, 37]. PG algorithms for variational
quantum algorithms have been shown to be robust against noise, and they can
systematically outperform other methods [34].

As a motivation to understand why policy gradient techniques may outperform
quasi-Newton and gradient-free optimizers as the number of qubits grows, we
trace some parallels between these different approaches and the existing
tradeoff in approximate inference methods [38, 39]. Two complementary
philosophies to the determination of a posterior probability given some
evidence are given by Markov Chain Monte-Carlo and Variational Inference
Methods. Variational inference is inherently biased, whereas MCMC’s bias
asymptotically vanishes. On the other hand, while variational inference is
fast, since it relies on the evaluation of a statistical distance (a proxy
thereof) at each iteration, convergence of MCMC is difficult to diagnose and
it yields spurious correlations in multimodal distributions unless it is
significantly downsampled, and as a result it can be painstainkingly slow.

Both in our variational compilation context and in the case of approximate
inference, a tradeoff between convergence rate and exactness must be chosen.
However, we argue that exactness is not as important in our case as we are
intrinsically dealing with an approximation problem, and because we can reduce
the bias by introducing more complicated parameterizations, which will
increase the expressivity of our model. Moreover, as it turns out, policy
gradient methods naturally lend themselves to paralelization, which allows for
a significant speed up compared to the time needed for typical numerical
optimizers (COBYLA, Powell and Nelder-Mead).

Figure 2: Expressivity of the model as a function of the chosen policy. A
given set of probability distributions parameterized by a vector $\xi$ is
contained in the larger set of all probability distributions. The statistical
distance which separates the learned distribution from the ground truth has
three components. The training error, which underlies the optimization method,
is the statistical distance between the current distribution
$\hat{p}_{\xi}(\theta)$ and the empirically optimal hypothesis
$\hat{h}_{\xi}(\theta)$. The distance between the empirically optimal
hypothesis and the optimal hypothesis $h^{*}_{\xi}(\theta)$ comes given by the
finiteness of sampling. Finally, the distance between $h^{*}(\theta)$ and the
ground truth $p^{*}(\theta)=\delta(\theta-\theta^{*})$ depends on the model
expressivity of the chosen parameterization.

In our case, the RL state will correspond to the quantum wavefunction (see
Appendix A). Fortunately, a simple PGRL algorithm known as REINFORCE [40] can
be slighly modified so that the reward is only related to the measured
fidelity at the end of each MonteCarlo rollout, meaning that we need not worry
about performing intermediate measurements which would collapse the
wavefunction and interrupt the computation. The reward will be the overlap
between the initial and the final state, given by Eq.(1). During this first
exploration of this approach, the set of possible actions under the current
policy, i.e. the angles of 1- and 2-qubit gates, will be randomly sampled from
a Gaussian distribution of the form:

$\theta\sim\pi(x;\mu,\Sigma)=\frac{1}{\sqrt{2\pi|\Sigma|}}e^{-(x-\mu)\Sigma^{-1}(x-\mu)^{T}}$
(2)

where the covariance matrix $\Sigma$ can either be fixed or obey some
exploration-exploitation schedule, or it can even be learnt (see Appendix A).
The objective function that we will use is therefore :

$J=\mathbb{E}_{\pi_{\mu}}[F]=\sum^{m}_{k}p_{k}\sum_{\theta}\pi(\theta|\mu,\Sigma)|\langle
k|V(\theta)^{\dagger}U|k\rangle|^{2}$ (3)

which corresponds to an average of the end reward, i.e the asymptotic
fidelity, over initial states (each sampled with probability $p_{k}$) and all
possible actions (given by the current policy $\pi(\theta|\mu,\Sigma)$). In
our case, maximizing $J$ corresponds to minimizing the associated cost
function
$J=\mathbb{E}_{\pi_{\mu}}[F]\stackrel{{\scriptstyle[}}{{|}}\Sigma|\rightarrow
0]{}{\longrightarrow}1-C$.

The gradient of this objective function in Eq.(3) is therefore:

$\nabla_{\mu}J=\sum^{m}_{k}p_{k}\sum_{\theta}\pi(\theta|\mu,\Sigma)\nabla_{\mu}\log\pi(\theta|\mu,\Sigma)|\langle
k|V(\theta)^{\dagger}U|k\rangle|^{2}$ (4)

where we have used the so-called “policy gradient theorem”, which amounts to
applying the chain rule to the policy function, and allows to write the
gradient of an expectation value as the expectation of a loglikelihood times a
cost function [28]. Estimating the gradient thus reduces to sampling a few
episodes, i.e. performing several Monte-Carlo tree searches, and evaluating
the cost function at the end of each one (see Appendix A).

In a barren plateau, $\nabla_{\theta}F\rightarrow 0$ [22, 24] and the value of
the cost function is exponentially flattened in the number of qubits. In our
case, we reformulate the problem in such a way that the space of solutions is
no longer a set of possible configurations, but rather a space of (Gaussian)
distributions. In the limiting case $|\Sigma|\rightarrow 0$, both
representations become equivalent. Therefore one expects that, in our RL
formulation, the gradient will consist of a (constant) cost times a gradient
of a parameterized function:

$\displaystyle\nabla_{\mu}J$ $\displaystyle=$
$\displaystyle\sum_{\theta}\sum^{m}_{k}p_{k}\pi(\theta|\mu,\Sigma)\nabla_{\mu}\log\pi(\theta|\mu,\Sigma)|\langle
k|V(\theta)^{\dagger}U|k\rangle|^{2}$ (5) $\displaystyle\approx$
$\displaystyle\mathbb{E}_{\pi_{\mu}}[\nabla_{\mu}\log\pi(\theta|\mu,\Sigma)]\times\textrm{constant}$

which shows that, even in the presence of a barren plateau, it should be
possible to utilize parameterized policies, like it is the case in policy
gradient, to further optimize the cost. It is equally clear that the minimum
of the objective function in Eq.(3) is not necessarily the same as in Eq.(1),
which means that this technique is inherently biased.

Our aim will therefore be to empirically verify that we can learn competitive
configurations of the variational parameters by optimizing a stochastic
policy, i.e. that we can reduce the sampling error as well as the training
error (see Fig.(2)) and to show that they are robust agains noise.

##  IV Numerical Experiments and Results

Figure 3: Reward (and cost) as a function of the iterations for different
methods and number of qubits. For a fixed number of iterations, policy
gradient achieves better performances than gradient-free optimizers in our
approximate compilation problem. (a) Typical training trajectories for PG-
based training of the variational circuit for compilation, for different
numbers of qubits.The mean and standard deviations over all Monte-Carlo
rollouts are shown. (b) Training instances for gradient-free optimizers for 5
qubits. In the noiseless case, Powell optimizer achieves higher rewards than
COBYLA on average, albeit it takes more than four times as many iterations.
Figure 4: Circuit trainability, evaluated through the asymptotic values of the
reward $J_{\infty}$ ($1-C_{\infty}$), as a function of the number of qubits
and the gathered statistics (5k, 10k, 50k and 100k numerical repetitions at
each iteration) and a fixed number of initial states,
$|\\{|k\rangle\\}|=n^{2}$. The solid lines correspond to the performance for
100k repetitions. Compilations of unitaries of logarithmic depth (2,3,4 and 5
for 5,10,15 and 20 qubits, respectively) While the onset of a barren plateau
is evident for gradient-free optimizers, REINFORCE remains competitive even
for up to 20 qubits. The reward (averaged fidelity) is well above $95\%$ even
for 15 qubits and we do not observe an exponential degradation in the
trainability of the algorithm. Increasing the experiment repetition rate from
5k to 100k only yields a moderate increase in fidelity, which indicates that
the sampling error is low (see Fig.(2)).

To assess the performance of PG methods applied to variational compilation, we
have run three different numerical experiments. We generated several random
shallow quantum circuits with depth logarithmic in the number of qubits and
known connectivity graph, which acted as the target unitary $U$, followed by a
circuit with the same connectivity graph and depth, and randomized parameters
implementing another unitary $V$. This setup is physically motivated because
in the absence of error correction, circuit depth of NISQ algorithms is
bounded by the inverse effective noise rate, which means that only shallow
circuits, i.e. of constant depth, can be realistically considered [14, 15].
Practically, our choice for this setup stems from the need to evaluate how
close the performance gets to its theoretical maximum. Given that most
unitaries have exponentially long circuits [1], sampling operators in
$SU(2^{n})$, instead of explicitely defining a quantum circuit, would almost
surely result in the optimization getting stuck at indeterminate values of the
cost function, which would in turn lead to a poor characterization of the
performance. Finally, since barren plateaus will arise frequently in
logarithmic and linear depth quantum circuits [23], we expect this choice of
circuit depths to still constitute a good testbed for our method.

Our methodology consisted of testing different gradient-free methods, such as
Powell and COBYLA, and a simple variation of REINFORCE PG, and to compare
their relative performance in training $V(\theta)$ to act as much as possible
as the inverse of the target unitary, i.e to determine $\theta^{*}$ such that
$V(\theta^{*})=U^{\dagger}$. We used the Cirq simulator for this purpose [41].

The first experiment was aimed at evaluating the trainability of variational
circuits, measured in terms of the obtained reward, as a function of the
number of qubits and for different repetition rates. We explored a logarithmic
depth regime, which is expected to suffer from the barren plateau effect [23].
The ability of PG to go deeper in the cost landscape (or, equivalently, to
achieve higher rewards) compared to gradient-free optimizers can be
intuitively understood as a consequence of an efficient reducition of the
training and the sampling errors for a fixed model expressivity (see Fig.(2)).
The finite sampling error is reduced by increasing the repetition rate at each
iteration, and it is straighforward to diagnose whether or not adding more
repetitions is beneficial. Reducing the training error involves tuning all the
hyperparameters of the policy gradient algorithm, for which there are several
prescriptions to follow [42] (see Appendix A). It is important to realize that
in REINFORCE, gradient ascent happens in the space of Gaussian means, rather
than in the space of angles. Indeed, these two formulations become equivalent
if the Gaussian distributions feature vanishing covariances. This crucial
difference makes it possible to calculate a gradient, since the policy can be
differentiated, and to take steps in the right direction (on average) in the
cost landscape. This should be compared with the situation of gradient-free
optimizers, which quickly do get stuck in barren plateaus [24]. We have
verified empirically that PG with Gaussian policies has better convergence
properties than gradient-free optimizers, as can be gleaned from Fig.(3): for
instance, it PG requires around 600 iterations to reach a fidelity of $90\%$,
whereas COBYLA needs slightly over 1000 iteration and Powell needs almost 2000
iterations.

Fig.(4) shows that the performance of policy gradient is degraded as the
number of qubits increases, albeit significantly less than for gradient-free
optimizers. Two sources of error will contribute mostly to this degradation:
hyperparameter selection and model expressivity. The meta-optimization of
hyperparameters, such as the exploration-exploitation schedule or the discount
factor for RMSprop, can become difficult as the number of parameters grows.
Importantly, the expressivity of the model, which depends on policies
capturing all relevant correlations between angles, will necessarily be
reduced, as Gaussian policies capture an ever-diminishing subset of all
possible distributions as the dimension of the problem increases.

Figure 5: Circuit trainability in the presence of noise. (a) Logarithmic depth
circuit simulations, as a function of the number of qubits and a depolarizing
noise strength of $n=0.01$ after each gate (5k repetitions per iteration). We
observe that the PG-based unitary compilation shares with gradient-free
methods a robustness to incoherent noise [43], and as a consequence larger
unitaries can be learnt. (b) Variance ratios of the asymptotic fluctuations in
$J_{\infty}$,
$\sigma^{\textrm{n}}_{J_{\infty},\textrm{R}}/\sigma^{\textrm{0}}_{J_{\infty},\textrm{R}}$
compared to the noiseless case, as a function of the number of qubits and the
repetition rate.

The second experiment seeked to estimate the robustness of the method in the
presence of noise, and to compare it to its gradient-free counterparts. Our
noise simulations involved appending a single-qubit depolarizing channel after
each gate. PG approximate compilation, being a stochastic optimization method,
handles fluctuations in the training process naturally. Recall that in the
case of REINFORCE applied to optimization of angles, the gradient is estimated
by performing a Monte-Carlo “tree search” over continuous actions and
computing the average reward. The existence of noise will amount to an
increased variance in the gradient estimator. To see this, consider that noise
can be modelled as a deviation from nominal values
$\theta_{\textrm{noisy}}=\theta_{\textrm{noiseless}}+\delta\theta$, with
$\delta\theta$ being a vector whose components are i.i.d random variables such
that $\mathbb{E}_{\theta}[\delta\theta]=0$. Then, a Taylor expansion gives:

$\displaystyle\nabla_{\xi}J|_{\textrm{noisy}}$ $\displaystyle=$
$\displaystyle\nabla_{\xi}J|_{\textrm{noiseless}}+$ (6) $\displaystyle
2\delta\theta$
$\displaystyle\sum_{\theta}\pi\nabla_{\xi}\log\pi\sum^{m}_{k}p_{k}\nabla_{\theta}|\langle
V^{\dagger}U\rangle|+\mathcal{O}(\delta^{2}\theta)$
$\displaystyle\stackrel{{\scriptstyle 1\ll m}}{{\longrightarrow}}$
$\displaystyle\nabla_{\xi}J|_{\textrm{noiseless}}$

in which $\nabla_{\theta}|\langle V^{\dagger}U\rangle|$ can account for both
coherent and incoherent noise and $\xi\in\\{\mu,\Sigma\\}$. Increasing the
number of measures can be shown to improve the trainability of the circuits.
In particular, we compared the standard deviations
$\sigma^{\textrm{n}}_{J_{\infty},\textrm{R}}$ of the asymptotic rewards
$J_{\infty}$ for different depolarizing noise strengths and repetition rates
(see Fig.(5)(b)).

Finally, in order to empirically validate that a polynomial number of
independent initial states suffices to correctly train the circuit, we used
several test sets to assess the performance of the learnt policy. We used a
reduced number of initial states (ranging from 100 to 100k) to train the
circuit and a set of test states, of fixed size 500, with little or no overlap
with the training set. The training set was obtained by randomly applying
Pauli operators, rotations around X and Y axis and Hadamard gates to every
qubit on the register, which amounts to local deformations of the $|0\rangle$
state. A tomographically complete set of measurements would require on the
order of $4^{10}\approx 10^{6}$ initial states. We tested the optimal policy
on three different test cases (see Fig. 6): (i) the $|0\rangle$ state, which
was not necesarily contained in the training set, and two sets of random
quantum states: (ii) one obtained from localy sampling Bloch vectors and (iii)
another by random sampling of the whole Hilbert space. This shows that
relatively few random local deformations of $|0\rangle$ state suffice to
attain a subtantial generalization of the algorithm for initial states that
has never seen before. Further optimizations of the training process are
needed to reduce the generalization error.

As a last remark, we would like to comment that fixed-state compilation can be
shown to be compatible with the requirements for fault tolerance. Indeed,
experiments with fixed-state compilation (i.e. $|\\{|k\rangle\\}|=1$) of
blocks of up to 5 qubits could be done within the fidelity requirements for
fault tolerance, for circuit depths of even linear depth in the number of
qubits. A similar calculation for 10 qubits gives a (fixed-state) fidelity of
$>99\%$ in contrast to the $<90\%$ that was achieved with gradient-free
methods.

Figure 6: Performance on various test sets, defined as
$F_{\textrm{test}}(\theta^{*})=1/T\sum^{T}_{k}|\langle
k|V(\theta^{*})^{\dagger}U|k\rangle|^{2}$, with
$\theta^{*}\sim\pi(x;\mu^{*},\Sigma\rightarrow 0)$ for fixed $T=500$, as a
function of size of training set $|\\{|k\rangle\\}|$. The numerical experiment
was carried out for 10 qubits. The thick dashed line represents the
performance on the training set, while the other three correspond to (i) state
$|k\rangle=|0\rangle$, which was not necessarily in the training set, (ii) a
set of locally randomized states, i.e.
$|k\rangle=|\eta_{1}\rangle\otimes|\eta_{2}\rangle\dots\otimes|\eta_{10}\rangle$,
with $|\eta_{j}\rangle$ a random Bloch vector in the XZ plane, and (iii) a set
of randomly sampled vectors of dimension $2^{10}$ (all of them with real
amplitudes). Since the initial states in the training set are local
modifications of $|0\rangle$, it is expected that the test fidelity of case
(ii) exceed that of random state vectors (iii), which are almost surely non-
local.

##  V Conclusion

We have introduced a method to compile variational quantum circuits using
techniques from reinforcement learning. This approach constitutes an
alternative to gradient-based and gradient-free methods, and has therefore the
potential to mitigate adverse effects at the onset of a barren plateau. This
is because the optimization is not performed over $\theta$-space, but rather
in the space of statistical distributions over $\theta$, for which it is
possible to implement classical gradient ascent methods for the distribution
parameters ($\mu$ and $\Sigma$ in our Gaussian case). These two descriptions
converge as the covariance of the Gaussian vanishes. Optimizing over
distributions rather than over variational parameters allowed us to find
competitive configurations even when gradient-free methods cease to be useful.

Another salient feature of PG for variational compilation is that its
performance is not degraded significantly in the presence of noise. This
crucial property allows us to achieve better configurations in simulations
with depolarizing noise strengths which are commesurate with state-of-the-art
gate fidelities [44, 45, 46, 47]. This method has the potential to be used as
an efficient way to train shallow quantum circuits (for up to about 10 qubits)
in the presence of noise, which could in turn be used as building blocks for
larger circuits. An intriguing question is whether this method could be used
as the lowest level of recursion in compilers reliant on the Solovay-Kitaev
construction, where it could be used to replace the need for a library of gate
sequences. We leave this approach for further work.

One obvious advantage of this method relates to training efficacy. An
iteration of the COBYLA or Powell optimizers typically takes far longer time
(between 10 and 100 times longer) than a REINFORCE iteration (which uses
parallel threads to evaluate the gradient). While it is difficult to compare
the runtimes of different approaches, convergence of PG-based approximate
compilation is typically faster both in terms of absolute time and number of
iterations needed to attain a given performance threshold.

We conjecture that the reward, i.e. the fidelity, will be ultimately limited
by the expressivity of the policy. Further optimizating of the hyperparameters
and using more complicated policies would be a way to verify this claim. Some
obvious improvements for policy gradient-based variational compilation would
include the adoption of a more expressive policy than the Gaussian policy we
have used, either by letting the algorithm learn the covariance matrix, or by
considering a mixture of Gaussians (which itself can take longer to train),
for example. Optimizing the exploration-exploitation schedule is another
possibility, such as it is done in simulated annealing optimization, as well
as introducing a moving-average calibration of the baseline which might
improve the convergence rate ([29], see appendix A), and perhaps using more
advanced gradient evaluation methods, such as proximal policy optimization.
These modification all have the potential to further increase the gap between
gradient-free optimizers and policy gradient methods for variational
compilation. Another interesting question is how to assess the bias incurred
by different parametric policies.

Finally, we expect this RL-based approach to circuit training to be beneficial
in quantum variational tasks other than unitary compilation that suffer from
the barren plateau effect.

## References

  * [1] Nielsen, M. A., & Chuang, I. Quantum Computation and Quantum Information (2002).
  * [2] Harrow, A. W., Recht, B., & Chuang, I. L. Efficient discrete approximations of quantum gates. Journal of Mathematical Physics, 43(9), 4445-4451 (2002).
  * [3] Dawson, C. M., & Nielsen, M. A. The Solovay-Kitaev algorithm. arXiv preprint quant-ph/0505030 (2005)
  * [4] Lin, H. W. Cayley graphs and complexity geometry. Journal of High Energy Physics, 2019(2), 1-15 (2019).
  * [5] Krioukov, D., Papadopoulos, F., Kitsak, M., Vahdat, A., & Boguná, M. Hyperbolic geometry of complex networks. Physical Review E, 82(3), 036106 (2010).
  * [6] Nielsen, M. A., Dowling, M. R., Gu, M., & Doherty, A. C. Quantum computation as geometry. Science, 311(5764), 1133-1135 (2006).
  * [7] Preskill, J. Quantum computing in the NISQ era and beyond. Quantum, 2, 79 (2018).
  * [8] Lloyd, S. Quantum approximate optimization is computationally universal. arXiv preprint arXiv:1812.11075 (2018).
  * [9] Morales, M. E., Biamonte, J. D., & Zimborás, Z. On the universality of the quantum approximate optimization algorithm. Quantum Information Processing, 19(9), 1-26 (2020).
  * [10] Kiani, B., Maity, R., & Lloyd, S. Learning unitaries via gradient descent optimization. Bulletin of the American Physical Society, 65 (2020).
  * [11] Farhi, E., & Harrow, A. W. Quantum supremacy through the quantum approximate optimization algorithm. arXiv preprint arXiv:1602.07674 (2016).
  * [12] Arute, F., Arya, K., Babbush, R., Bacon, D., Bardin, J. C., Barends, R., … & Martinis, J. M. Quantum supremacy using a programmable superconducting processor. Nature, 574(7779), 505-510 (2019).
  * [13] Zhu, Q., Cao, S., Chen, F., Chen, M. C., Chen, X., Chung, T. H., … & Pan, J. W. Quantum Computational Advantage via 60-Qubit 24-Cycle Random Circuit Sampling. arXiv preprint arXiv:2109.03494 (2021).
  * [14] Bravyi, S., Gosset, D., & König, R. Quantum advantage with shallow circuits. Science, 362(6412), 308-311 (2018).
  * [15] Bravyi, S., Gosset, D., Koenig, R., & Tomamichel, M. Quantum advantage with noisy shallow circuits. Nature Physics, 16(10), 1040-1045 (2020).
  * [16] Bauer, B., Bravyi, S., Motta, M., & Chan, G. K. L. Quantum algorithms for quantum chemistry and quantum materials science. Chemical Reviews, 120(22), 12685-12717 (2020).
  * [17] O’Malley, P. J., Babbush, R., Kivlichan, I. D., Romero, J., McClean, J. R., Barends, R., … & Martinis, J. M. Scalable quantum simulation of molecular energies. Physical Review X, 6(3), 031007 (2016).
  * [18] Ralli, A., Love, P. J., Tranter, A., & Coveney, P. V. Implementation of measurement reduction for the variational quantum eigensolver. Physical Review Research, 3(3), 033195 (2021).
  * [19] Hastings, M. B. Classical and quantum bounded depth approximation algorithms. arXiv preprint arXiv:1905.07047 (2019).
  * [20] Bravyi, S., Kliesch, A., Koenig, R., & Tang, E. Obstacles to variational quantum optimization from symmetry protection. Physical Review Letters, 125(26), 260505 (2020).
  * [21] Bravyi, S., Kliesch, A., Koenig, R., & Tang, E. Hybrid quantum-classical algorithms for approximate graph coloring. arXiv preprint arXiv:2011.13420 (2020).
  * [22] McClean, J. R., Boixo, S., Smelyanskiy, V. N., Babbush, R., & Neven, H. Barren plateaus in quantum neural network training landscapes. Nature communications, 9(1), 1-6 (2018).
  * [23] Cerezo, M., Sone, A., Volkoff, T., Cincio, L., & Coles, P. J. Cost-function-dependent barren plateaus in shallow quantum neural networks. arXiv e-prints, arXiv-2001 (2020).
  * [24] Arrasmith, A., Cerezo, M., Czarnik, P., Cincio, L., & Coles, P. J. Effect of barren plateaus on gradient-free optimization. Quantum, 5, 558 (2021).
  * [25] Khatri, S., LaRose, R., Poremba, A., Cincio, L., Sornborger, A. T., & Coles, P. J. Quantum-assisted quantum compiling. Quantum, 3, 140 (2019).
  * [26] Sharma, K., Khatri, S., Cerezo, M., & Coles, P. J. Noise resilience of variational quantum compiling. New Journal of Physics, 22(4), 043006 (2020).
  * [27] Holmes, Z., Arrasmith, A., Yan, B., Coles, P. J., Albrecht, A., & Sornborger, A. T. Barren plateaus preclude learning scramblers. Physical Review Letters, 126(19), 190501 (2021).
  * [28] Sutton, R. S., & Barto, A. G. Reinforcement learning: An introduction. MIT press (2018).
  * [29] Fösel, T., Tighineanu, P., Weiss, T., & Marquardt, F. Reinforcement learning with neural networks for quantum feedback. Physical Review X, 8(3), 031084 (2018).
  * [30] Porotti, R., Essig, A., Huard, B., & Marquardt, F. Deep Reinforcement Learning for Quantum State Preparation with Weak Nonlinear Measurements. arXiv preprint arXiv:2107.08816 (2021).
  * [31] Nautrup, H. P., Delfosse, N., Dunjko, V., Briegel, H. J., & Friis, N. Optimizing quantum error correction codes with reinforcement learning. Quantum, 3, 215 (2019).
  * [32] August, M., & Hernández-Lobato, J. M. . Taking gradients through experiments: LSTMs and memory proximal policy optimization for black-box quantum control. In International Conference on High Performance Computing (pp. 591-613). Springer (2018)
  * [33] Garcia-Saez, A., & Riu, J. Quantum observables for continuous control of the quantum approximate optimization algorithm via reinforcement learning. arXiv preprint arXiv:1911.09682 (2019).
  * [34] Yao, J., Bukov, M., & Lin, L.. Policy gradient based quantum approximate optimization algorithm. In Mathematical and Scientific Machine Learning (pp. 605-634). PMLR (2020)
  * [35] Yao, J., Lin, L., & Bukov, M. Reinforcement Learning for Many-Body Ground State Preparation based on Counter-Diabatic Driving. arXiv preprint arXiv:2010.03655 (2020).
  * [36] Moro, L., Paris, M. G., Restelli, M., & Prati, E. Quantum Compiling by Deep Reinforcement Learning. arXiv preprint arXiv:2105.15048 (2021).
  * [37] He, Z., Li, L., Zheng, S., Li, Y., & Situ, H. Variational quantum compiling with double Q-learning. New Journal of Physics, 23(3), 033002 (2021).
  * [38] Blei, D. M., Kucukelbir, A., & McAuliffe, J. D. Variational inference: A review for statisticians. Journal of the American statistical Association, 112(518), 859-877 (2017).
  * [39] Koller, D., & Friedman, N. Probabilistic graphical models: principles and techniques. MIT press (2009).
  * [40] Williams, R. J. Simple statistical gradient-following algorithms for connectionist reinforcement learning. Machine learning, 8(3), 229-256 (1992).
  * [41] Cirq, a python framework for creating, editing, and invoking noisy intermediate scale quantum NISQ circuits
github.com/quantumlib/Cirq

  * [42] Shahriari, B., Swersky, K., Wang, Z., Adams, R. P., & De Freitas, N. Taking the human out of the loop: A review of Bayesian optimization. Proceedings of the IEEE, 104(1), 148-175 (2015).
  * [43] Colless, J. I., Ramasesh, V. V., Dahlen, D., Blok, M. S., Kimchi-Schwartz, M. E., McClean, J. R., … & Siddiqi, I. Computation of molecular spectra on a quantum processor with an error-resilient algorithm. Physical Review X, 8(1), 011021 (2018).
  * [44] Barends, R., Kelly, J., Megrant, A., Veitia, A., Sank, D., Jeffrey, E., … & Martinis, J. M. Superconducting quantum circuits at the surface code threshold for fault tolerance. Nature, 508(7497), 500-503 (2014).
  * [45] Yang, C. H., Chan, K. W., Harper, R., Huang, W., Evans, T., Hwang, J. C. C., … & Dzurak, A. S. Silicon qubit fidelities approaching incoherent noise limits via pulse engineering. Nature Electronics, 2(4), 151-158 (2019).
  * [46] Huang, W., Yang, C. H., Chan, K. W., Tanttu, T., Hensen, B., Leon, R. C. C., … & Dzurak, A. S. Fidelity benchmarks for two-qubit gates in silicon. Nature, 569(7757), 532-536 (2019).
  * [47] Schäfer, V. M., Ballance, C. J., Thirumalai, K., Stephenson, L. J., Ballance, T. G., Steane, A. M., & Lucas, D. M. Fast quantum logic gates with trapped-ion qubits. Nature, 555(7694), 75-78 (2018).
  * [48] Tieleman, T., Hinton, G. Lecture 6.5-rmsprop: Divide the gradient by a running average of its recent magnitude. COURSERA: Neural networks for machine learning (2012).

Appendices

## Appendix A: Policy Gradient with Gaussian Policies and RMSprop

The first step towards phrasing our approximate compilation problem as a
reinforcement learning problem is to map the process of training a variational
quantum circuit to a Markov Decision Process (MDP). Different mappings will
lead to different learning scenarios in which distinct facets or challenges of
the original quantum task will become apparent. An MDP consists of a tuple
$(\mathcal{S},\mathcal{A},R_{s,a},P_{s^{\prime},(s,a)},)$, with $\mathcal{S}$
the set of states, $\mathcal{A}$ the set of actions, $R$ is the reward
obtained by taking action $a$ in state $s$, and finally $P$ is a stochastic
matrix giving the probability of transitioning to state $s^{\prime}$ given
that the current state is $s$ and the current action is $a$. Generally, P is
so large that it can only be sampled by an agent exploring an environment. An
agent seeking to maximize the long-term reward of an MDP can do so by
optimizing a policy $\pi(a,s)$, which associates a probability to each
available action-state pair.

There exist several approaches to optimizing a policy. Temporal Difference
methods aim at measuring the reward after each transition and update the value
of each state under the current policy (the estimated long-term reward
associated to that state). Optimality of the corresponding policies is ensured
by the Bellmann Optimality Condition [28]. Another approach is given by the
direct optimization of the policy, thus relying little, or not at all, on
value iteration. One of the simplest policy gradient algorithms, conceptually
as well as in terms of implementation, is the REINFORCE algorithm [40]. In
REINFORCE, the policy is parameterized and belongs to a variational family of
distributions, such that it is possible to differentiate it with respect to
the variational parameters.

In this work we have used a variation of the REINFORCE algorithm in which the
reward will only be obtained at the end of a given sequence of actions (an
episode). The rewards gathered at the end of each episode are used to estimate
the gradient as follows: at each iteration the agent performs a Monte-Carlo
tree search, i.e. it explores the space of actions for a fixed amount of time.

$J=\mathbb{E}_{\pi_{\mu,\Sigma}}[R]=\sum_{\theta}\sum^{m}_{s}p_{s}\pi(\theta|\mu,\Sigma)r^{(s)}_{\theta},$
(7)

where we start at state s with probability $p_{s}$ and the reward
$r^{(s)}_{\theta}=|\langle s|V(\theta)^{\dagger}U|s\rangle|^{2}$ is related to
the end-state fidelity. Now the problem is to compute the gradient of this
expectation value. Each particular trajectory in the Monte-Carlo tree search,
a _rollout_ in the RL jargon, contributes with a weight which is proportional
to its associated final reward. This can be seen by applying the chain rule to
the policy function:

$\nabla_{\xi}J=\sum_{\theta}\sum^{m}_{s}p_{s}\pi(\theta|\mu,\Sigma)\nabla_{\xi}\log\pi(\theta|\mu,\Sigma)r^{(s)}_{\theta}$
(8)

where $\xi\in\\{\mu,\Sigma\\}$ and we have used the Gaussian policy:

$\pi(x;\mu,\Sigma)=\frac{1}{\sqrt{2\pi|\Sigma|}}e^{-(x-\mu)\Sigma^{-1}(x-\mu)^{T}}$
(9)

The “Log-likelihood trick” allows us to express a gradient of an expectation
value as the expectation value of a different gradient, which can be estimated
numerically. The logarithm of the Gaussian policy has the following gradients:

$\displaystyle\nabla_{\mu}\log\pi(x;\mu,\Sigma)$ $\displaystyle=$
$\displaystyle\Sigma^{-1}(x-\mu)$ (10)
$\displaystyle\nabla_{\Sigma}\log\pi(x;\mu,\Sigma)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\Sigma^{-1}(\mathbf{1}-(x-\mu)^{T}(x-\mu)\Sigma^{-1})$
(11)

which allow for learning both the mean $\mu$ and the covariance $\Sigma$.
However, in this work we do not learn $\Sigma$ but rather fix a simple
exploration-exploitation schedule $\Sigma(t)=(1-t/T)\Sigma_{i}+t/T\Sigma_{f}$,
such that $\Sigma_{i}\gg\Sigma_{f}\longrightarrow 0$.

A further consideration is the choice of the baseline in the computation of
the gradient, which allows to reduce the variance of the estimator. Provided
the baseline is positively correlated with the end reward, i.e. if
$Cov[R,b]>0$, the advantage $R-b$ will have less variance as the baseline will
compensate the fluctuations of R. We have used a simple mean, i.e.
$b=\frac{1}{m}\Sigma_{\theta}r_{\theta}$, which intuitively gives an advantage
for the rollouts which perform better than average.

Once the gradient has been estimated, the updates are performed following a
_RMSprop_ gradient update rule [48]. RMSprop is an adaptive learning rate
method with better convergence properties than simple gradient ascent methods.
It works by (i) computing a moving discounted average of the gradient
variances and (ii) dividing the update step by the discounted variance. The
result is that the learning rate will increase in relatively flat landscape
directions and it will decrease in steep landscapes. Let $\sigma^{2}_{g}$ be
the variance of the computed gradients over different rollouts.

$\displaystyle\sigma^{(t)}_{g}$ $\displaystyle=$
$\displaystyle\gamma\sigma^{(t-1)}_{g}+(1-\gamma)(\nabla_{\xi}J|_{g})^{2}$
(12) $\displaystyle\xi$ $\displaystyle\leftarrow$
$\displaystyle\xi+\eta\frac{\nabla_{\xi}J|_{g}}{\sqrt{\sigma^{(t)}_{g}+\varepsilon}}$
(13)

where $\gamma$ is the discount factor. It has been empirically shown to allow
for a more efficient exploration of complex cost lanscapes. A list with the
parameters used for this work is provided below:

Parameter | Value
---|---
$\Sigma_{i}$ | $diag(10^{-2})$
$\Sigma_{f}$ | $diag(10^{-5})$
$N_{rollouts}$ | 20-25
$p_{s}$ | 1/m
m | $\max(15n,n^{2})$
$\gamma$ | 0.9
$\eta$ | $(2.5-3)\times 10^{-3}$
$\varepsilon$ | $10^{-8}$
$maxiter_{\textrm{REINFORCE}}$ | $2000-5000$
$maxiter_{\textrm{COBYLA}}$ | $10000$
$maxiter_{\textrm{Powell}}$ | $50000$

## Appendix B: Efficient Trainability without Tomographically Complete
Measurements

It is in principle possible to efficiently calculate the fidelity at the end
of each Monte-Carlo rollout by means of the so-called Hilbert-Schmidt test
[25, 26], which makes use of an array of ancilla qubits and a series of
potentially long range interactions.

In the absence of a secondary qubit register, one needs to resort to an
estimation of the fidelity via Eq.(1). In the noiseless case, it is possible
to recover the exact value of the fidelity
$F(\theta)=\frac{1}{4^{n}}\sum^{\sim 4^{n}}_{k}|\langle
k|V(\theta)^{\dagger}U|k\rangle|^{2}$ using a tomographically complete
characterization.

An important consideration, therefore, is to ensure that our estimate of the
fidelity $\hat{F}(\theta)$, which makes use of a massively downsampled
subensemble of orthogonal states, is sufficient to train our variational
circuit.

We would like to make sure that our error is dominated by the expressivity of
the policy, rather than by an intrinsically flawed estimation of the fidelity.
The set of initial states, $\\{|k\rangle\\}^{m}$ with $\langle
k_{i}|k_{j}\rangle=\delta_{ij}$, allows to estimate the fidelity across
independent directions of the operator space, i.e it allows us to build an
estimator $\hat{F}=\frac{1}{m}\sum^{m}_{k}f_{i}$, where $0\leq f_{i}\leq 1$
can be interpreted as independent samples of the true fidelity. In the limit
where the fluctuations due to finite sampling vanish (infinitely many
repetitions for each Monte-Carlo rollout), the Hoeffding inequality implies:

$P(|\hat{F}-F|\geq\epsilon)\leq 2\exp-2\epsilon^{2}m$ (14)

In the limit of a large number of qubits, small differences of fidelity become
increasinly significative and difficult to obtain through numerical
optimization, which amounts to imposing a scaling on the estimator accuracy,
$\epsilon\sim n^{-t}$ for some positive real $t$. Being able to attain the
limit $|\hat{F}-F^{*}|\rightarrow 0$, i.e. to satisfy the regime
$\epsilon^{2}m\gg 1$ for any number of qubits gives a scaling of
$m\sim\Omega(n^{2t})$.