* [57] Bernhard Schölkopf, Ralf Herbrich, and Alex J Smola.  A generalized representer theorem.  In International conference on computational learning theory, pages 416–426. Springer, 2001.
  * [58] Charles A Micchelli, Yuesheng Xu, and Haizhang Zhang.  Universal kernels.  Journal of Machine Learning Research, 7(12), 2006.
  * [59] Ryan Shaffer, Lucas Kocia, and Mohan Sarovar.  Surrogate-based optimization for variational quantum algorithms, 2022\.
  * [60] Maria Schuld and Francesco Petruccione.  Quantum Models as Kernel Methods, pages 217–245.  Springer International Publishing, Cham, 2021.
  * [61] Adriano Barenco, Charles H. Bennett, Richard Cleve, David P. DiVincenzo, Norman Margolus, Peter Shor, Tycho Sleator, John A. Smolin, and Harald Weinfurter.  Elementary gates for quantum computation.  Phys. Rev. A, 52:3457–3467, Nov 1995.
  * [62] Alexei Gilchrist, Daniel R. Terno, and Christopher J. Wood.  Vectorization of quantum operations and its use, 2011.
  * [63] Christopher K Williams and Carl Edward Rasmussen.  Gaussian processes for machine learning, volume 2.  MIT press Cambridge, MA, 2006.
  * [64] Karol Życzkowski and Hans-Jürgen Sommers.  Average fidelity between random quantum states.  Phys. Rev. A, 71:032313, Mar 2005.
  * [65] Yunchao Liu, Srinivasan Arunachalam, and Kristan Temme.  A rigorous and robust quantum speed-up in supervised machine learning.  Nature Physics, Jul 2021.
  * [66] Lukasz Cincio, Yiğit Subaşı, Andrew T Sornborger, and Patrick J Coles.  Learning the quantum algorithm for state overlap.  New Journal of Physics, 20(11):113022, 2018.
  * [67] Donald R Jones, Matthias Schonlau, and William J Welch.  Efficient global optimization of expensive black-box functions.  Journal of Global optimization, 13(4):455–492, 1998.
  * [68] Daniel James Lizotte.  Practical bayesian optimization.  PhD thesis, University of Alberta, 2008.
  * [69] Eric Brochu, Vlad M. Cora, and Nando de Freitas.  A tutorial on bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning, 2010\.
  * [70] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.
  * [71] IBM Quantum.  Qiskit: An open-source framework for quantum computing, 2021.
  * [72] J.C. Spall.  Multivariate stochastic approximation using a simultaneous perturbation gradient approximation.  IEEE Transactions on Automatic Control, 37(3):332–341, 1992.
  * [73] IBM Quantum team.  ibmq_quito, 2021.
  * [74] Johnnie Gray.  quimb: a python library for quantum information and many-body calculations.  Journal of Open Source Software, 3(29):819, 2018.
  * [75] Lucas Slattery, Benjamin Villalonga, and Bryan K. Clark.  Unitary block optimization for variational quantum algorithms, 2021.
  * [76] Andrea Skolik, Jarrod R. McClean, Masoud Mohseni, Patrick van der Smagt, and Martin Leib.  Layerwise learning for quantum neural networks.  Quantum Machine Intelligence, 3(1):5, Jan 2021.
  * [77] D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac.  Matrix product state representations.  Quantum Info. Comput., 7(5):401–430, jul 2007.
  * [78] F. Verstraete, V. Murg, and J.I. Cirac.  Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems.  Advances in Physics, 57(2):143–224, 2008.
  * [79] Ulrich Schollwöck.  The density-matrix renormalization group in the age of matrix product states.  Annals of Physics, 326(1):96–192, 2011.  January 2011 Special Issue.
  * [80] Román Orús.  A practical introduction to tensor networks: Matrix product states and projected entangled pair states.  Annals of Physics, 349:117–158, 2014.
  * [81] Sebastian Paeckel, Thomas Köhler, Andreas Swoboda, Salvatore R. Manmana, Ulrich Schollwöck, and Claudius Hubig.  Time-evolution methods for matrix-product states.  Annals of Physics, 411:167998, 2019.
  * [82] Jacob C Bridgeman and Christopher T Chubb.  Hand-waving and interpretive dance: an introductory course on tensor networks.  Journal of Physics A: Mathematical and Theoretical, 50(22):223001, may 2017.
  * [83] Robert E. Kass and Adrian E. Raftery.  Bayes factors.  Journal of the American Statistical Association, 90(430):773–795, 1995.
  * [84] Hector J. Garcia, Igor L. Markov, and Andrew W. Cross.  Efficient inner-product algorithm for stabilizer states, 2012.
  * [85] Scott Aaronson and Daniel Gottesman.  Improved simulation of stabilizer circuits.  Physical Review A, 70(5), nov 2004.
  * [86] Sergey Bravyi, Dan Browne, Padraic Calpin, Earl Campbell, David Gosset, and Mark Howard.  Simulation of quantum circuits by low-rank stabilizer decompositions.  Quantum, 3:181, September 2019.
  * [87] Richard H Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou Zhu.  A limited memory algorithm for bound constrained optimization.  SIAM Journal on scientific computing, 16(5):1190–1208, 1995.
  * [88] Jakob Bossek, Carola Doerr, and Pascal Kerschke.  Initial design strategies and their effects on sequential model-based optimization: An exploratory case study based on bbob.  In Proceedings of the 2020 Genetic and Evolutionary Computation Conference, GECCO ’20, page 778–786, New York, NY, USA, 2020. Association for Computing Machinery.
  * [89] Guifré Vidal.  Efficient classical simulation of slightly entangled quantum computations.  Phys. Rev. Lett., 91:147902, Oct 2003.

##  Appendix A Derivation of feature spaces

The state kernel for the ansatz in (4) is given by:

$\begin{split}k_{s}(\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})&=\Tr{U(\boldsymbol{\theta}^{\prime})\rho_{0}U^{\dagger}(\boldsymbol{\theta}^{\prime})U(\boldsymbol{\theta})\rho_{0}U^{\dagger}(\boldsymbol{\theta})}\\\
&=\llangle
U(\boldsymbol{\theta}^{\prime})\rho_{0}U^{\dagger}(\boldsymbol{\theta}^{\prime})|U(\boldsymbol{\theta})\rho_{0}U^{\dagger}(\boldsymbol{\theta})\rrangle\\\
&=\llangle\rho_{0}\rvert(U^{\dagger}(\boldsymbol{\theta}^{\prime})\otimes
U^{T}(\boldsymbol{\theta}^{\prime}))(U(\boldsymbol{\theta})\otimes
U^{*}(\boldsymbol{\theta}))\lvert\rho_{0}\rrangle,\end{split}$ (30)

where $\rho_{0}=\ket{\boldsymbol{0}}\bra{\boldsymbol{0}}$, $\lvert A\rrangle$
denotes the vectorization of operator $A$ and we have made use of the
identities $\llangle A|B\rrangle=\Tr{A^{\dagger}B}$ and $\lvert
ABC\rrangle=A\otimes C^{T}\lvert B\rrangle$. At this point, we could stop and
identify the state kernel’s feature vectors as quantum states on
$\mathcal{H}\otimes\mathcal{H}^{*}$ (where $\mathcal{H}$ is the Hilbert space
of the qubits) of the form $U(\boldsymbol{\theta})\otimes
U^{*}(\boldsymbol{\theta})\lvert\rho_{0}\rrangle$. However in doing so would
have to construct these states explicitly, making analysis and comparison to
classical kernels difficult. By expanding each Pauli rotation as
$e^{-iP_{i}\theta_{i}/2}=\cos(\frac{\theta_{i}}{2})I-i\sin(\frac{\theta_{i}}{2})P_{i}$
we can express $U(\boldsymbol{\theta})$ in a convenient form as the product of
$k$ inner-products of 2-dimensional vectors
$\boldsymbol{u}_{q}=(R_{q},-iP_{q}R_{q})^{T}$ and
$\boldsymbol{a}(\theta)=(\cos{\theta},\sin{\theta})^{T}$

$U(\boldsymbol{\theta})=R_{p+1}(\boldsymbol{u}^{T}_{p}\boldsymbol{a}(\theta_{p}/2))\dots(\boldsymbol{u}^{T}_{1}\boldsymbol{a}(\theta_{1}/2)).$
(31)

This allows us to separate out the parameter-dependent part of the kernel from
the ansatz-dependent part. We can then write the vectorized unitary
$U(\boldsymbol{\theta})\otimes U^{*}(\boldsymbol{\theta})$ as

$U(\boldsymbol{\theta})\otimes U^{*}(\boldsymbol{\theta})=R_{p+1}\otimes
R_{p+1}^{*}(\boldsymbol{u}^{T}_{p}\boldsymbol{a}(\theta_{p}/2)\otimes(\boldsymbol{u}_{p}^{*})^{T}\boldsymbol{a}(\theta_{p}/2))\dots(\boldsymbol{u}^{T}_{1}\boldsymbol{a}(\theta_{1}/2)\otimes(\boldsymbol{u}_{1}^{*})^{T}\boldsymbol{a}(\theta_{1}/2)).$
(32)

The tensor products of the form
$\boldsymbol{u}^{T}_{q}\boldsymbol{a}(\theta_{q}/2)\otimes(\boldsymbol{u}_{q}^{*})^{T}\boldsymbol{a}(\theta_{q}/2)$
can be expanded and simplified into the following inner product

$\boldsymbol{u}^{T}_{q}\boldsymbol{a}(\theta_{q}/2)\otimes(\boldsymbol{u}_{q}^{*})^{T}\boldsymbol{a}(\theta_{q}/2)=\frac{1}{2}\begin{pmatrix}(I\otimes
I+P_{q}\otimes P_{q}^{*})R_{q}\otimes R_{q}^{*}\\\ (I\otimes
iP_{q}^{*}-iP_{q}\otimes I)R_{q}\otimes R_{q}^{*}\\\ (I\otimes I-P_{q}\otimes
P_{q}^{*})R_{q}\otimes R_{q}^{*}\end{pmatrix}^{T}\begin{pmatrix}1\\\
\sin{\theta_{q}}\\\
\cos{\theta_{q}}\end{pmatrix}=\boldsymbol{s}_{q}^{T}\boldsymbol{v}_{q}(\theta_{q}),$
(33)

where $\boldsymbol{v}_{q}=(1,\sin{\theta_{q}},\cos{\theta_{q}})^{T}$ and
$\boldsymbol{s}_{q}$ is the vector of sums of unitaries above. The overall
vectorized unitary is

$U(\boldsymbol{\theta})\otimes U^{*}(\boldsymbol{\theta})=R_{p+1}\otimes
R_{p+1}^{*}(\boldsymbol{s}_{p}^{T}\boldsymbol{v}_{p}(\theta_{p}))\dots(\boldsymbol{s}_{1}^{T}\boldsymbol{v}_{1}(\theta_{1}))=R_{p+1}\otimes
R_{p+1}^{*}\boldsymbol{s}^{T}\boldsymbol{v}(\boldsymbol{\theta}),$ (34)

where $\boldsymbol{s}$ is given by

$\boldsymbol{s}=\frac{1}{2^{p}}\left[\begin{pmatrix}(I\otimes I+P_{p}\otimes
P_{p}^{*})R_{p}\otimes R_{p}^{*}\\\ (I\otimes iP_{p}^{*}-iP_{p}\otimes
I)R_{p}\otimes R_{p}^{*}\\\ (I\otimes I-P_{p}\otimes P_{p}^{*})R_{p}\otimes
R_{p}^{*}\end{pmatrix}\otimes_{K}\dots\otimes_{K}\begin{pmatrix}(I\otimes
I+P_{1}\otimes P_{1}^{*})R_{1}\otimes R_{1}^{*}\\\ (I\otimes
iP_{1}^{*}-iP_{1}\otimes I)R_{1}\otimes R_{1}^{*}\\\ (I\otimes I-P_{1}\otimes
P_{1}^{*})R_{1}\otimes R_{1}^{*}\end{pmatrix}\right]$ (35)

and $\boldsymbol{v}(\boldsymbol{\theta})$ by

$\boldsymbol{v}(\boldsymbol{\theta})=\begin{pmatrix}1\\\ \sin(\theta_{p})\\\
\cos(\theta_{p})\end{pmatrix}\otimes_{K}\dots\otimes_{K}\begin{pmatrix}1\\\
\sin(\theta_{1})\\\ \cos(\theta_{1})\end{pmatrix}.$ (36)

Here we have made a distinction between the tensor products $\otimes$ that
denote operators on the physical and conjugated multi-qubit Hilbert spaces and
the Kronecker products (for vectors with matrix entries) $\otimes_{K}$ between
the $k$-“sub-vectors” which determine the dimension of the resulting kernels’
feature spaces. The state kernel is then

$\begin{split}k_{s}(\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})&=\llangle\rho_{0}\rvert\boldsymbol{v}^{T}(\boldsymbol{\theta}^{\prime})\boldsymbol{s}^{*}\boldsymbol{s}^{T}\boldsymbol{v}(\boldsymbol{\theta})\lvert\rho_{0}\rrangle\\\
&=\boldsymbol{v}^{T}(\boldsymbol{\theta}^{\prime})\llangle\rho_{0}\rvert\boldsymbol{s}^{*}\boldsymbol{s}^{T}\lvert\rho_{0}\rrangle\boldsymbol{v}(\boldsymbol{\theta})\\\
&=\boldsymbol{v}^{T}(\boldsymbol{\theta}^{\prime})\boldsymbol{S}\boldsymbol{v}(\boldsymbol{\theta}),\end{split}$
(37)

where each element of $\boldsymbol{s}^{*}$ is the Hermitian conjugate of the
corresponding element in $\boldsymbol{s}$ so that
$(\boldsymbol{S})_{ij}=\llangle\rho_{0}\rvert
s_{i}^{\dagger}s_{j}\lvert\rho_{0}\rrangle$. (Note that the unitary
$R_{p+1}\otimes R_{p+1}^{*}$ has been cancelled by its Hermitian conjugate.)

The unitary kernel evaluated for this ansatz takes the following form:

$\begin{split}k_{u}(\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})&=\frac{1}{d^{2}}\absolutevalue{\Tr{U^{\dagger}(\boldsymbol{\theta}^{\prime})U(\boldsymbol{\theta})}}^{2}\\\
&=\frac{1}{d^{2}}\Tr{U^{\dagger}(\boldsymbol{\theta}^{\prime})U(\boldsymbol{\theta})}\Tr{U^{T}(\boldsymbol{\theta}^{\prime})U^{*}(\boldsymbol{\theta})}\\\
&=\frac{1}{d^{2}}\Tr{[U^{\dagger}(\boldsymbol{\theta}^{\prime})\otimes
U^{T}(\boldsymbol{\theta}^{\prime})][U(\boldsymbol{\theta})\otimes
U^{*}(\boldsymbol{\theta})]}.\end{split}$ (38)

Using the decomposition of this unitary given in (34), we can write the
unitary kernel as

$\begin{split}k_{u}(\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})&=\frac{1}{d^{2}}\Tr{\boldsymbol{v}^{T}\boldsymbol{s}^{*}\boldsymbol{s}^{T}\boldsymbol{v}(\boldsymbol{\theta})}\\\
&=\boldsymbol{v}^{T}(\boldsymbol{\theta}^{\prime})\frac{\Tr{\boldsymbol{s}^{*}\boldsymbol{s}^{T}}}{d^{2}}\boldsymbol{v}(\boldsymbol{\theta})\\\
&=\boldsymbol{v}^{T}(\boldsymbol{\theta}^{\prime})\boldsymbol{T}\boldsymbol{v}(\boldsymbol{\theta}),\end{split}$
(39)

where again each element of $(\boldsymbol{s}^{*})_{i}=s_{i}^{\dagger}$ so that
$\boldsymbol{T}={\Tr{\boldsymbol{s}^{*}\boldsymbol{s}^{T}}}/{d^{2}}$. (Again
note that the $R_{p+1}$ terms have cancelled.)

Using the form for $U(\boldsymbol{\theta})\otimes U^{*}(\boldsymbol{\theta})$
given in (34) we also can now write the noiseless energy function
$E(\boldsymbol{\theta})$ in terms of the $\boldsymbol{v}(\boldsymbol{\theta})$
vector, and thereby the feature space vectors, as

$\begin{split}E(\boldsymbol{\theta})&=\Tr{HU(\boldsymbol{\theta})\rho_{0}U^{\dagger}(\boldsymbol{\theta})}=\llangle
H\rvert U(\boldsymbol{\theta})\otimes
U^{*}(\boldsymbol{\theta})\lvert\rho_{0}\rrangle\\\ &=\llangle
R_{p+1}^{\dagger}HR_{p+1}\rvert\boldsymbol{s}^{T}\boldsymbol{v}(\boldsymbol{\theta})\lvert\text{vec}(\rho_{0})\rrangle\\\
&=\llangle
R_{p+1}^{\dagger}HR_{p+1}\rvert\boldsymbol{s}^{T}\lvert\text{vec}(\rho_{0})\rrangle\boldsymbol{v}(\boldsymbol{\theta})\\\
&=\boldsymbol{h}^{T}\boldsymbol{v}(\boldsymbol{\theta}),\text{ where
}\boldsymbol{h}=\llangle
R_{p+1}^{\dagger}HR_{p+1}\rvert\boldsymbol{s}\lvert\text{vec}(\rho_{0})\rrangle.\end{split}$
(40)

We can also use (40) to calculate the mean value of the energy function
$\mathbb{E}_{\boldsymbol{\theta}}[E(\boldsymbol{\theta})]$ with respect to the
gate angles, which is required to properly set the Gaussian process’ mean.
Because $E(\boldsymbol{\theta})$ is a linear function of
$\boldsymbol{v}(\theta)$ we have

$\mathbb{E}_{\boldsymbol{\theta}}[E(\boldsymbol{\theta})]=\boldsymbol{h}^{T}\mathbb{E}_{\boldsymbol{\theta}}[\boldsymbol{v}(\boldsymbol{\theta})].$
(41)

$\boldsymbol{v}(\boldsymbol{\theta})=\boldsymbol{v}^{(k)}\otimes_{K}\dots\otimes_{K}\boldsymbol{v}^{(1)}$
is a Kronecker product of $k$ subvectors with
$\boldsymbol{v}^{(j)}=(1,\sin(\theta_{j}),\cos(\theta_{j}))^{T}$ each of which
only depends on a single gate angle (we assume all gate angles are
independent). This means that

$\begin{split}\mathbb{E}_{\boldsymbol{\theta}}[\boldsymbol{v}(\boldsymbol{\theta})]=\frac{1}{{(2\pi)}^{p}}\int
d^{p}\\!\boldsymbol{\theta}\
\boldsymbol{v}(\boldsymbol{\theta})&=\frac{1}{{(2\pi)}^{p}}\left(\int_{0}^{2\pi}\\!\boldsymbol{v}^{(k)}d\theta_{p}\right)\otimes_{K}\
\dots\otimes_{K}\left(\int_{0}^{2\pi}\\!\boldsymbol{v}^{(1)}d\theta_{1}\right)\\\
&=\begin{pmatrix}1\\\ 0\\\
0\end{pmatrix}\otimes_{K}\dots\otimes_{K}\begin{pmatrix}1\\\ 0\\\
0\end{pmatrix}.\end{split}$ (42)

So the only the first element in
$\mathbb{E}_{\boldsymbol{\theta}}[\boldsymbol{v}(\boldsymbol{\theta})]$ is
nonzero, greatly simplifying the expected energy to

$\mathbb{E}_{\boldsymbol{\theta}}[E(\boldsymbol{\theta})]=h_{0}=\llangle
R_{p+1}^{\dagger}HR_{p+1}\rvert\frac{(R_{p}\otimes R_{p}^{*}+P_{p}R_{p}\otimes
P_{p}^{*}R_{p}^{*})}{2}\dots\frac{(R_{1}\otimes R_{1}^{*}+P_{1}R_{1}\otimes
P_{1}^{*}R_{1}^{*})}{2}\lvert\text{vec}(\rho_{0})\rrangle.$ (43)

One way to interpret $\frac{1}{2}{(R_{q}\otimes R_{q}^{*}+P_{q}R_{q}\otimes
P_{q}^{*}R_{q}^{*})}$ is as a projector onto the positive eigenspace of the
Pauli operator $P_{q}\otimes P_{q}^{*}$. Alternatively it can be seen as
taking the part of its unvectorized input (following the application of
$R_{q}$) that commutes with $P_{q}$. We can decompose an input
$\rho^{\prime}=R_{q}\rho R_{q}^{\dagger}$ into parts (sums of Pauli operators)
that commute and anticommute with $P$, i.e. $\rho^{\prime}=c_{q}+a_{q}$ with
$[P_{q},c_{q}]=0$ and $\\{P_{q},a_{q}\\}=0$ then $\frac{1}{2}{(R_{q}\otimes
R_{q}^{*}+P_{q}R_{q}\otimes
P_{q}^{*}R_{q}^{*})}\lvert\rho\rrangle=\frac{1}{2}\lvert R_{q}\rho
R_{q}^{\dagger}+P_{q}R_{q}\rho R_{q}^{\dagger}P_{q}\rrangle=\lvert
c_{q}\rrangle$. Applying each of these projectors in this way we see that

$\mathbb{E}_{\boldsymbol{\theta}}[E(\boldsymbol{\theta})]=h_{0}=\Tr{R^{\dagger}_{p+1}HR_{p+1}c_{1,\dots,k}}$
(44)

where $c_{1,\dots,k}$ is what remains of $\rho_{0}$ once the process of
applying projector (applying each $R_{q}$ and removing the part of the result
that anticommutes with $P_{q}$). In general calculating this is complicated as
it depends on the fixed unitaries $\\{R_{q}\\}$ present in the ansatz however
we suggest that for most ansatzes of interest the remainder will be
$c_{1,\dots,k}=I/2^{n}$ as this factor is present in all density matrices
(ensuring $\Tr{\rho}=1$) and commutes with any $P_{q}$ regardless of the fixed
unitaries applied. This leads to our suggested approximation to the mean
energy of
$\mathbb{E}_{\boldsymbol{\theta}}[E(\boldsymbol{\theta})]=h_{0}\approx\Tr{R_{p+1}^{\dagger}HR_{p+1}I/2^{n}}=\Tr{H}$.
It takes only a few anti-commuting parameterized Pauli rotations on each qubit
(e.g. two sequential non-commuting single-qubit PPRs on each qubit) for this
approximation to be exact.

##  Appendix B Limits on state kernel feature space dimension

In section VI we identified the unnormalized states
$\\{s_{i}\lvert\rho_{0}\rrangle\\}$ as vectorizations of Hermitian operators
on $\mathcal{H}$. We can show this explicitly by considering the action of the
$q^{\mathrm{th}}$ subvector of $\boldsymbol{s}$ (5), $\boldsymbol{s}^{(q)}$,
on a vectorized input $\lvert\rho\rrangle$ (with Hermitian $\rho$):

$\boldsymbol{s}^{(q)}\lvert\rho\rrangle=\frac{1}{2}\begin{pmatrix}(I\otimes
I+P_{q}\otimes P_{q}^{*})R_{q}\otimes R_{q}^{*}\\\ (I\otimes
iP_{q}^{*}-iP_{q}\otimes I)R_{q}\otimes R_{q}^{*}\\\ (I\otimes I-P_{q}\otimes
P_{q}^{*})R_{q}\otimes
R_{q}^{*}\end{pmatrix}\lvert\rho\rrangle=\frac{1}{2}\begin{pmatrix}\lvert
R_{q}\rho R_{q}^{\dagger}+P_{q}R_{q}\rho R_{q}^{\dagger}P_{q}\rrangle\\\
\lvert iR_{q}\rho R_{q}^{\dagger}P-iP_{q}R_{q}\rho R_{q}^{\dagger}\rrangle\\\
\lvert R_{q}\rho R_{q}^{\dagger}-P_{q}R_{q}\rho
R_{q}^{\dagger}P_{q}\rrangle\end{pmatrix}$ (45)

As in the discussion given in given in Appendix A on the mean energy, if we
write $\rho^{\prime}=R_{q}\rho R_{q}^{\dagger}$ as a sum of a commuting and an
anti-commuting part $P$, i.e. $\rho^{\prime}=c_{q}+a_{q}$ with
$[P_{q},c_{q}]=0$ and $\\{P_{q},a_{q}\\}=0$, then we can rewrite the resulting
vector as

$\boldsymbol{s}^{(q)}\lvert\rho\rrangle=\begin{pmatrix}\lvert c_{q}\rrangle\\\
\lvert ia_{q}P_{q}\rrangle\\\ \lvert a_{q}\rrangle\end{pmatrix}.$ (46)

So the $q^{\mathrm{th}}$ subvector of $\boldsymbol{s}$ acts to separate the
input $\rho^{\prime}$ (transformed by $R_{q}$) into the parts which commute
(first entry) and anti-commute (last entry) with $P_{q}$ while the middle
entry contains the anti-commuting part post-multiplied by $iP_{q}$. Like the
input, all of these operators are Hermitian. The first and last entries are
formed by separating Pauli operators in the decomposition of the input into
commuting and anti-commuting subsets; this means that if $\rho^{\prime}$ has
support on (can be decomposed into) a finite set of $m$ Pauli operators,
$\rho^{\prime}=\sum_{i=1}^{m}P_{i}$, then Pauli operators that do not feature
in this decomposition can only be introduced in the middle entry. We also note
that $c_{q}$, $ia_{q}P_{q}$, and $a_{q}$ are necessarily linearly independent
as it can be shown using the cylic property of the trace, $PP=I$, and
$P_{q}a_{q}=-a_{q}P_{q}$, that they are mututally orthogonal under the
Hilbert-Schmidt norm $\langle A,B\rangle_{HS}=\Tr{A^{\dagger}B}$; $\langle
a_{q},c_{q}\rangle_{HS}=-\langle a_{q},c_{q}\rangle_{HS}=0$, $\langle
ia_{q}P_{q},c_{q}\rangle_{HS}=-\langle ia_{q}P_{q},c_{q}\rangle_{HS}=0$,
$\langle ia_{q}P_{q},a_{q}\rangle_{HS}=-\langle
ia_{q}P_{q},a_{q}\rangle_{HS}=0$.

The dimension of the feature space is determined by the number of
parameterized rotations (increasing the dimension of
$\boldsymbol{v}(\boldsymbol{\theta})$) and the rank of the Gram matrix
$\boldsymbol{S}$. Each unnormalized state $s_{i}\lvert\rho_{0}\rrangle=\lvert
O_{i}\rrangle$ is a vectorization of a Hermitian operator $O_{i}$ meaning the
inner products forming $S_{ij}=\llangle\rho_{0}\rvert
s_{i}^{\dagger}s_{j}\lvert\rho_{0}\rrangle$ can be written
$S_{ij}=\Tr{O_{i}O_{j}}$, the Hilbert-Schmidt inner-product between $O_{i}$
and $O_{j}$. To maximise the dimension of $\mathcal{F}_{s}$ we must ensure
that the $\\{O_{i}\\}$ for different entries of
$\boldsymbol{s}\lvert\rho_{0}\rrangle$ are linearly independent (spanning a
subspace of Hermitian operators of dimension $\geq 3^{p}$, the number of
entries).

Suppose after the $j^{\mathrm{th}}$ subvector has been applied, the vector
$s^{(j)}\otimes_{K}\dots\otimes_{K}s^{(1)}\lvert\rho_{0}\rrangle$ has $m$
linearly independent elements/operators. We can write each element as a sum of
an operator that commutes with the next parameterized Pauli rotation $P_{j+1}$
and one that anti-commutes:
$s^{(j)}\otimes_{K}\dots\otimes_{K}s^{(1)}\lvert\rho_{0}\rrangle=(\lvert
a_{1}+c_{1}\rrangle,\dots,\lvert a_{m}+c_{m}\rrangle)^{T}$. Applying the
subvector for $P_{j+1}$ splits each element into $3$ new elements, e.g.
$\lvert a_{1}+c_{1}\rrangle\to(\lvert c_{1}\rrangle,\lvert
ia_{1}P_{j+1}\rrangle,\lvert a_{1}\rrangle)$, increasing the feature space
dimension by a factor of $3$. However, these elements will not always be
linearly independent. Clearly, the contribution to the number of linearly
independent elements from the commuting terms $c_{1}\dots c_{m}$ is at most of
size $m$. From the anti-commuting terms we separate out introduce the elements
$a_{1},\dots,a_{m}$ and $ia_{1}P_{j+1},\dots,ia_{m}P_{j+1}$, which also anti-
commute with $P_{j+1}$. An $n$-qubit Pauli operator anti-commutes with
$4^{n}/2$ other elements of the $n$-qubit Pauli group (ignoring phases) and
commutes with the other remaining $4^{n}/2$ elements. As there are only
$4^{n}/2$ Pauli operators that anti-commute with $P_{j+1}$, the span of the
combined set of $a_{1},\dots,a_{m},ia_{1}P_{j+1},\dots,ia_{m}P_{j+1}$ is at
most $4^{n}/2$ dimensional. This means that the maximal increase in the
feature space dimension possible when adding another parameterized Pauli
rotation to the circuit is $4^{n}/2$ (as we can go from $m$ independent
elements to at most $m+4^{n}/2$). Denoting the maximal feature space dimension
of $\mathcal{F}_{s}$ for $n$ qubits and a $p$ parameterized Pauli rotation
ansatz as $d_{s}^{(n,p)}$ we have the recursion relation

$d_{s}^{(n,p)}=\min(4^{n},3d_{s}^{(n,p-1)},4^{n}/2+d_{s}^{(n,p-1)}).$ (47)

We define the base-case $d_{s}^{n,0}=1$ for all $n$ – this is purely for
completeness as a kernel on a zero-dimensional input (no parameterized
rotations) has no meaning.

##  Appendix C Limits on unitary kernel feature space dimension

Like the state kernel, the unitary kernel’s feature space can be expressed as
a quadratic form involving a parameter-dependent vector
$\boldsymbol{v}(\boldsymbol{\theta})$ and a Gram matrix $\boldsymbol{T}$ with
elements $T_{ij}=\Tr{s_{i}^{\dagger}s_{j}}$. Like with the state kernel, the
dimension of the feature space $\mathcal{F}_{u}$ is determined by the number
of parameters (determines the size of $\boldsymbol{v}$) and the rank of the
$\boldsymbol{T}$ matrix.

The $T$ matrix is comprised of Hilbert-Schmidt inner products between the
elements of the vector $\boldsymbol{s}$. If the elements of $\boldsymbol{s}$
were completely generic $2^{2n}\times 2^{2n}$ matrices then the maximal rank
of $\boldsymbol{T}$ would be $4^{2n}$ as there are this many orthogonal (under
the Hilbert-Schmidt inner product) basis matrices, e.g. the $2n$-qubit Pauli
matrices $\\{I,X,Y,Z\\}^{\otimes 2n}$ or the standard basis matrices
$e_{ij}=\outerproduct{i}{j}$. However there is additional structure present in
the elements of $\boldsymbol{s}$ that reduces the maximal rank of
$\boldsymbol{T}$. For example, all elements in $\boldsymbol{s}$ are symmetric
under interchange of the physical ($\mathcal{H}$) and conjugated
($\mathcal{H}^{*}$) subsystems followed by complex conjugation and involve
sums and products of various unitary and Hermitian operators. To derive the
maximal scaling of the unitary kernel’s feature space dimension we will first
examine weighted sums of unitary matrices of the form $A\otimes A^{*}$,
matrices $\boldsymbol{M}(\boldsymbol{w},\boldsymbol{A})$ given by

$\boldsymbol{M}(\boldsymbol{w},\boldsymbol{A})=\sum_{k}w_{k}A^{(k)}\otimes{A^{(k)}}^{*}$
(48)

with ${A^{(k)}}^{\dagger}A^{(k)}=I_{d}$ (and each $A$ has dimension $d\times
d$). In the following we will demonstrate that the maximal dimension of
$\mathcal{F}_{u}$ is ultimately set by the number of linearly independent
entries of such matrices.

For a single $A\otimes A^{*}$, the unitarity constraint implies that
$\sum_{k}A_{ki}A^{*}_{kj}=\sum_{k}A_{ik}A^{*}_{jk}=\delta_{ij}$ – these
constraints can be easily related to $A\otimes A^{*}$ as they correspond to
contractions over pairs of indices. A general element of this matrix is
$(A\otimes A^{*})_{ijkl}=A_{ij}A^{*}_{kl}$ so contraction over the first and
third indices gives $\sum_{p}(A\otimes
A^{*})_{kikj}=\sum_{p}A_{ki}A^{*}_{kj}=\delta_{ij}$ (the first set of
unitarity constraints) while contraction over the second and fourth indices
gives $\sum_{p}(A\otimes A^{*})_{ikjk}=\sum_{p}A_{ik}A^{*}_{jk}=\delta_{ij}$
(the second set of constraints). This appears to give us $2d^{2}$ equations,
each of which will reduce the number of linearly independent elements by $1$
giving a total of $d^{4}-2d^{2}$ independent elements. However, because the
constraints for $i=j$ all sum to the same value,
$\sum_{p}A_{ik}A^{*}_{ik}=\sum_{p}A_{ki}A^{*}_{ki}=1$ for all $i$, this fixed
value acts as another free parameter in each set of constraints, increasing
the number of linearly independent elements by $2$ to a total of
$d^{4}-2(d^{2}-1)$. When considering weighted sums of these matrices as given
in (48) these constraints become

$\sum_{p}\sum_{p}(w_{p}A^{(p)}\otimes{A^{(p)}}^{*})_{ikjk}=\sum_{p}\sum_{p}(w_{p}A^{(p)}\otimes{A^{(p)}}^{*})_{kikj}=(\sum_{p}w_{p})\delta_{ij}$
(49)

The linear dependence of the terms in the constraints with $i=j$ is only
defined up to specification of the sum of the weights meaning that the
quantity $\sum_{p}w_{p}$ is an additional free parameter, giving the same
result.

The initial fixed unitary $R_{1}\otimes R_{1}^{*}$ is in the form we have just
discussed and for an $n$-qubit ansatz will in general have
$d^{4}-2(d^{2}-1)=4^{2n}-2(4^{n}-1)$ linearly independent elements. For
unitary operators $B$ and $C$, the product $(B\otimes B^{*})(C\otimes
C^{*})=A\otimes A^{*}$ (with $A=BC$) is also a matrix of the form discussed
above. As a result, the Pauli-containing part of the top/bottom subvector
elements $I\otimes I\pm P_{1}\otimes P_{1}^{*}$ acts on $R_{1}\otimes
R_{1}^{*}$ to give matrices that admit a decomposition as in (48) and so have
the same number of linearly independent elements as before (as will the result
from applying the next fixed unitary $R_{2}\otimes R_{2}^{*}$ and the
subsequent top/bottom elements of the next subvector). All that is left is to
consider how the application of $i(I\otimes P^{*}-P\otimes I)$, from the
middle elements of each subvector, to a general unitary $A\otimes A^{*}$
(which could be part of a sum as in (48)) affects the unitarity constraints
(49) that determine the number of linearly independent elements. We have:

$\begin{split}\sum_{k}(i(I\otimes P^{*}-P\otimes I)(A\otimes
A^{*}))_{ikjk}&=i\sum_{p}\sum_{m}({A}_{ik}P^{*}_{jm}A^{*}_{mk}-P_{im}A_{mk}A^{*}_{jk})\\\
&=i\sum_{m}(\delta_{im}P^{*}_{jm}-P_{im}\delta_{mj})=i(P^{*}_{ji}-P_{ij})\\\
&=i(P^{\dagger}-P)_{ij}=0,\end{split}$ (50)

while for the other constraint we get the same result

$\begin{split}\sum_{p}(i(I\otimes P^{*}-P\otimes I)(A\otimes
A^{*}))_{kikj}&=i\sum_{p}\sum_{m}({A}_{ki}P^{*}_{km}A^{*}_{mj}-P_{km}A_{mi}A^{*}_{kj})\\\
&=i\sum_{p}\sum_{m}(P_{mk}{A}_{ki}A^{*}_{mj})-i\sum_{m}\sum_{p}(P_{km}A_{mi}A^{*}_{kj})\\\
&=0\end{split}$ (51)

(in both cases we have used that $P^{\dagger}=P$). By linearity, the same
results will be obtained for an arbitrary weighted sum as in (48). These are
identical to the constraints in (49) for $i\neq j$ but give $0$ rather than a
fixed value ($1$ for a single $A\otimes A$ or $\sum_{p}w_{p}$ for a weighted
sum) for $i=j$. Having the $i=j$ case evaluate to $0$ means that the number of
linearly independent elements in matrices formed of weighted sums of terms
like this is $d^{4}-2(d^{2})=4^{2n}-2(4^{n})$. However, because these new
constraints are a subset of the original constraints (49) with
$\sum_{p}w_{p}=0$, the over-all number of linearly independent elements in the
matrices that form the entries of $\boldsymbol{s}$ and therefore the maximal
rank of $\boldsymbol{T}$ is given by $4^{2n}-2(4^{n}-1)$. We verified this
scaling with simulations for $n=1,2,\text{ and }3$ qubits using Haar random
unitaries for $\\{R_{q}\\}$ to ensure minimal structure in the ansatz that
might limit the feature space dimension.