\section{\label{sec:level1}Introduction}

The field of quantum machine learning (QML) is currently in its infancy and there exist, not only questions of advantages over its classical counter-part, but also questions of its realisability in real-world applications. Many early QML approaches \cite{Wiebe_2012, lloyd2013quantum} utilised the HHL algorithm \cite{Harrow_2009} to obtain speed-ups on the linear algebraic computations behind many learning methods.
However, such methods presume efficient preparation of quantum states or \textit{quantum access} to data, limiting practicality \cite{10.1038/nphys3272}.

Currently, the two leading contenders for near-term supervised QML methods are quantum neural networks (QNNs) \cite{10.1103/physreva.98.032309, farhi2018classification} and quantum kernel methods \cite{10.1038/s41586-019-0980-2, 10.1103/physrevlett.122.040504}.
QNNs constructed from parameterised quantum circuits embed data into a \textit{quantum feature space} and train parameters to minimise a loss function of some observable of the circuit. However, such variational algorithms face problems of barren plateaus in optimisation \cite{10.1038/s41467-018-07090-4, 10.1038/s41467-021-21728-w, PRXQuantum.3.010313} which hinder the feasibility of training the model. In contrast, quantum kernel methods suggest a non-parametric approach whereby only the inner product of quantum embedded data points are estimated on a quantum device -- a process known as \textit{quantum kernel estimation} (QKE) \cite{10.1038/s41586-019-0980-2, 10.1103/physrevlett.122.040504}. We refer to the model utilising QKE for constructing a \textit{support vector machine} (SVM), as a \textit{quantum}-SVM (QSVM) -- the most commonly explored quantum kernel method.
Despite rigorous proofs of quantum advantage using a QSVM \cite{10.1038/s41567-021-01287-z}, they are not without their limitations for arbitrary \textit{real-world} datasets. Vanishing kernel elements \cite{10.1038/s41534-021-00498-9} subsequently require a large number of quantum circuit samples (suggesting a deep connection with barren plateaus \cite{inductive_bias_of_qk}) and indications that quantum kernels fail to generalise with the addition of qubits without the careful encoding of the correct problem-dependent \textit{inductive bias} \cite{inductive_bias_of_qk}.
Furthermore, the kernel approach quadratically grows in complexity with the number of training samples, as opposed to the linear scaling of QNNs. This becomes untenable for large datasets.
These limitations indicate that further work is required to practically deploy QML techniques.

The close similarity of QNNs and quantum kernel methods \cite{qnns_are_kernel} highlight the importance of developing alternate QML models that present distinct model hypotheses.
In this work, we propose a \textit{quantum random forest} (QRF) that fundamentally cannot be expressed through a kernel machine due to its discontinuous decision tree structure. Such a model is inspired by the classical \textit{random forest} (CRF) \cite{10.1023/A:1010933404324}.

The CRF has shown to be an extremely successful model in machine learning applications. However, there exist many barriers permitting the entry of an analogous QRF. The main challenge is the determination of an algorithm that is not limited by either the intrinsically quantum complication of the \textit{no-cloning theorem} \cite{10.1017/cbo9780511976667.010}, or the hard problem of efficiently loading data in superposition \cite{10.1038/nphys3272}.
As opposed to previous works \cite{10.1007/s11128-013-0687-5, https://doi.org/10.48550/arxiv.2112.13346} that employ Grover's algorithm with input states of specific form, the QRF proposed in this work is ultimately an ensemble of weak classifiers referred to as \textit{quantum decision trees} (QDTs) with an approximate QSVM forming the \textit{split function}  -- shown in Figure \ref{fig:qrf_diagram}. Crucially, each QDT has random components that distinguishes it from the ensemble. The aim is to form a collection of uncorrelated weak QDT classifiers so that each classifier provides distinct additional information. The possible benefits of an ensemble of quantum classifiers have been expressed \cite{10.1038/s41598-018-20403-3}, with recent suggestions of boosting techniques for ensembles of QSVMs \cite{boosting_qsvm}. However, the distinctive tree structure has not yet been explored. Not only does this allow for discontinuous model hypotheses, we will see that the structure naturally accommodates multi-class supervised problems without further computation. In comparison, QSVMs would require \textit{one-against-all} (OAA) and \textit{one-against-one} (OAO) \cite{10.7551/mitpress/4175.001.0001} strategies which would respectively require constructing $|\mathcal{C}|$ and $|\mathcal{C}|(|\mathcal{C}|-1)/2$ separate QSVMs, where $\mathcal{C}$ is the set of classes in the problem.

The expressive advantage of the QRF is predicated on the distinguishability of data when mapped to a higher dimensional feature space. This proposition is supported by Cover's Theorem \cite{10.1109/pgec.1965.264137}, which states that data is more likely to be separable when projected to a sparse higher dimensional space. Cover's result is also a principal motivator for developing both quantum and classical kernel methods. However, it is often the case that generalisation suffers at the expense of a more expressive model -- a phenomena commonly known through the bias-variance trade-off in statistical learning. To counter the expressivity of the (already quite expressive) quantum model, we make a randomised low-rank Nystr\"{o}m approximation \cite{using_nystr_speed_kern} of the kernel matrix supplied to the SVM, limiting the ability of the model to overfit. We refer to this strategy as \textit{Nystr\"{o}m}-QKE (NQKE). We show that the inclusion of this approximation further allows for a reduced circuit sampling complexity.

The paper is structured as follows: in Section \ref{section:qrf_main} we address the exact form of the QRF model, before providing theoretical results in \ref{sec:theory}. This is followed by numerical results and discussion in Section \ref{section:results}.
Background to many of the ML techniques are presented in Appendix \ref{section:background}.

\section{The Quantum Random Forest model}\label{section:qrf_main}

\begin{figure*}
    \centering
    \includegraphics[width=500pt]{Figures/qrf_diagram.pdf}

    \caption{The quantum random forest (QRF) constitutes an ensemble of classifiers, referred to as quantum decision trees (QDTs). This tree structure is in turn a directed graph structure of nodes defined by the split function  $\mathcal{N}_{\Phi, L}:\mathbb{R}^D \xrightarrow{}\{-1,+1\}$ and $\Phi, L$ labels the type of function. The split function is an SVM that admits a Nystr\"{o}m approximated quantum kernel, $\widetilde{K}$. The kernel is defined through a chosen embedding $\Phi$ where $K_{ij}=|\langle \Phi(x_j) |\Phi(x_i)\rangle|^2$ and $\ket{\Phi(x_i)}=\mathcal{U}(x_i)\ket{0}$ for parameterised unitary $\mathcal{U}$. The randomly selected $L$ points for Nystr\"{o}m approximation determine the hyperplane generated in the feature space, $\mathcal{H}_\Phi$.}
    \label{fig:qrf_diagram}
\end{figure*}

The Quantum Random Forest (QRF) is composed of a set of $T$ distinct quantum decision trees (QDTs), $\{ \mathbf{Q}_t \}_{t=1}^{T}$, which form the weak independent classifiers of the QRF ensemble. This is illustrated in Figure \ref{fig:qrf_diagram}. As with all supervised ML methods, the QDTs are given the set of $N$ training instances, $x^{\mathrm{train}}_i \in \mathbb{R}^D$, and their associated $n_c$ class labels, $y^{\mathrm{train}}_i \in \mathcal{C} = \{0, 1, ..., n_c - 1 \}$, to learn from the annotated data set, $\mathcal{S} = \{ (x^{\mathrm{train}}_i, y^{\mathrm{train}}_i)\}_{i=1}^{N}$ sampled from some underlying distribution $\mathcal{D}$.
Each tree is trained using $N_p \leq N$ -- which we will refer to as the partition size -- instances sampled from the original dataset $\mathcal{S}$. This acts as both a regularisation strategy, as well as a way to reduce the time complexity of the model. Such a technique is known as bagging \cite{10.5555/3086952} or bootstrap aggregation and it has the ability to greatly minimise the probability of overfitting.

Once each classifier is trained, the model is tested with a separate set, $\mathcal{S}^{\mathrm{test}}$ of size $N'$, from which we compare the predicted class distribution with the real. To obtain the overall prediction of the QRF model, predictions of each QDT is averaged across the ensemble. Specifically, each QDT returns a probability distribution, $\mathbf{Q}_t(\vec{x}; c) = \mathrm{Pr}(c | \vec{x})$ for $c \in \mathcal{C}$, with the overall prediction as the class with the greatest mean across the ensemble. Hence, given an input, $\vec{x}$, we have prediction of the form,
\begin{equation}\label{eq:voting_qrf}
    \widetilde{y} (\vec{x}) = \argmax_{c \in \mathcal{C}} \Bigg\{  \frac{1}{T}\sum_{t=1}^T  \mathbf{Q}_t(\vec{x}; c)  \Bigg\} .
\end{equation}
\noindent The accuracy of the model can be obtained with the comparison of predictions with the real labels, $1/N'\sum_{(\vec{x}, y) \in \mathcal{S}^{\mathrm{test}}}  \delta_{y, \widetilde{y}(\vec{x})}$ where $\delta_{i, j}$ is the Kronecker delta. The distinction between training and testing data is generally critical to ensure that the model is not overfitting to the trained dataset and hence generalises to the prediction of previously unknown data points. We now elaborate on the structure of the QDT model.

\subsection{A Quantum Decision Tree}

The Quantum Decision Tree (QDT) proposed in this work is an independent  classification algorithm that has a directed graph structure identical to that of a binary decision tree -- illustrated in Figure \ref{fig:qrf_diagram} and elaborated further in Appendix \ref{section:rf_model}. Each vertex, also referred to as a \textit{node}, can either be a split or leaf node, distinguished by colour in Figure \ref{fig:qrf_diagram}. A leaf node determines the classification output of the QDT, while a split node separates inserted data points into two partitions that subsequently continue down the tree. The split effectiveness is measured by the reduction in entropy, referred to as the \textit{information gain} (IG). Given a labelled data set, $\mathcal{S}$, that is split into partitions $\mathcal{S}_L$ and $\mathcal{S}_R$ by the split function, we have,
\begin{equation}\label{eq:ig}
    \mathrm{IG}(\mathcal{S}; \mathcal{S}_L , \mathcal{S}_R ) = H(\mathcal{S}) - \sum_{i \in \{L, R\}} \frac{|\mathcal{S}_i|}{|\mathcal{S}|} H(\mathcal{S}_i),
\end{equation}
\noindent where $H$ is the entropy over the class distribution defined as $H (\mathcal{S}) = -\sum_{c \in \mathcal{C}} \mathrm{Pr}(c) \log_2 \mathrm{Pr}(c)$ where class $c \in \mathcal{C}$ occurs with probability $\mathrm{Pr}(c)$ in the set $\mathcal{S}$. Clearly the IG will increase as the splitting more distinctly separates instances of different classes. Using information gain is advantageous especially when more than two classes are involved and an accuracy-based effectiveness measure fails.

Given a partitioned training set at the root node, the QDT carries out this splitting process to the \textit{leaves} of the tree where a prediction of class distribution is made from the proportions of remaining data points in each class. Mathematically, with a subset of data points $\mathcal{S}^{(l)}$ supplied at training to leaf $\ell^{(l)}$ indexed by $l$, we set its prediction as the probability distribution,
\begin{equation}\label{eq:leaf_prob}
  \ell^{(l)}(\mathcal{S}^{(l)}; c) = \frac{1}{\lvert \mathcal{S}^{(l)} \rvert}\sum_{(\vec{x}, y) \in \mathcal{S}^{(l)}} \big[ y = c  \big] ,
\end{equation}
\noindent where $[p ]$ is the Iversion bracket that returns 1 when the proposition $p$ is true and 0 otherwise. When training the model, a node is defined a \textit{leaf} if any of the following conditions are met, (i) the training instances supplied to the node are of a single class -- in which case further splitting is unnecessary, (ii) when the number of data points left for splitting is smaller than a user defined value, $m_s$, ie. $|\mathcal{S}^{(l)}| \leq m_s$, or (iii) the node is at the maximum depth of the tree, $d$, defined by the user. Once trained, prediction occurs by following an instance down the tree until it reaches a leaf node. The unique path from root to leaf is referred to as the \textit{evaluation path}. At this end point, the prediction of the QDT is the probability distribution held within the leaf, defined at training in equation \eqref{eq:leaf_prob}. The ensemble of predictions made by distinct QDTs are subsequently combined with Eq. \eqref{eq:voting_qrf}.

The crux of the QDT model lies at the performance of split nodes, with the $l^{\text{th}}$ node defined through a split function, $\mathcal{N}^{(l)}_{\theta}: \mathcal{X} \xrightarrow{} \{-1, +1\}$, where $\theta$ are (hyper-) parameters of the split function, that are either manually selected prior to, or algorithmically tuned during, the training phase of the model -- the specifics are elaborated in the next section. Instances are divided into partitions $\mathcal{S}_- = \{ (\vec{x}, y) \in \mathcal{S} | \mathcal{N}^{(l)}_{\theta}(\vec{x}) = -1 \}$ and $\mathcal{S}_+ = \{ (\vec{x}, y) \in \mathcal{S} | \mathcal{N}^{(l)}_{\theta}(\vec{x}) = 1 \}$. At this point it is natural to refer to the information gain of a split function as, $\mathrm{IG}(\mathcal{S}|\mathcal{N}^{(l)}_{\theta}):=\mathrm{IG}(\mathcal{S}; \mathcal{S}_- , \mathcal{S}_+ )$ with the aim being to maximise $\mathrm{IG}(\mathcal{S}|\mathcal{N}_{\theta})$ at each node.
In essence, we desire a tree that is able to separate points so that instances of different classes terminate at different leaf nodes. The intention being that a test instance that follows the unique evaluation path to a particular node is most likely to be similar to the instances that followed the same path during training. Hence, effective split functions amplify the separation of classes down the tree. However, it remains important that the split function generalises, as repeated splitting over a deep tree can easily lead to overfitting.
It is for this reason decision trees often employ simple split functions to avoid having a highly expressive model. Namely, the commonly used CART algorithm \cite{li1984classification} has a threshold function for a single attribute at each node.

In the next section we specify the split function to have the form of a support vector machine (SVM) that employs a quantum kernel. The core motivation being to generate a separating hyperplane in a higher dimensional quantum feature space that can more readily distinguish instances of different classes. However, this is a complex method that will inevitably lead to overfitting.
Hence, we take an approximate approach to reduce its effectiveness without stifling possible quantum advantage.

\subsection{Nystr\"{o}m Quantum Kernel Estimation}

The proposed split function generates a separating hyperplane through a kernelised SVM. The kernel matrix (or Gram matrix) is partially computed using Quantum Kernel Estimation (QKE) \cite{10.1103/physrevlett.113.130503} and completed using the Nystr\"{o}m approximation method. We therefore refer to this combined process as \textit{Nystr\"{o}m Quantum Kernel Estimation} (NQKE) with its parameters determining the core performance of a splitting node. Background on SVMs and a deeper discussion on the Nystr\"{o}m method are supplied in Appendix \ref{section:svm}.

Given a dataset $\mathcal{S}^{(i)} = \{(x_j , y_j)\}_{j=1}^{N^{(i)}}$ to the $i^{\text{th}}$ split node, the Nystr\"{o}m method randomly selects a set of $L$ landmarks points $\mathcal{S}^{(i)}_L \subseteq \mathcal{S}^{(i)}$ . Without loss of generality we may assume $\mathcal{S}^{(i)}_L = \{(x_j , y_j)\}_{j=1}^{L}$. The inner product between elements in $\mathcal{S}^{(i)}$ and $\mathcal{S}^{(i)}_L$ are computed using a quantum kernel defined through the inner product of parameterised density matrices,
\begin{equation}
    k_{\Phi}(x', x'') = \Tr [\rho_{\Phi}(x') \rho_{\Phi}(x'')],
\end{equation}
where we have the spectral decomposition $\rho_{\Phi}(x) = \sum_j \lambda_j \ket{\Phi_j(x)} \bra{ \Phi_j(x)}$ with parameterised pure states $\{\ket{\Phi_j(x)}\}_j$ that determine the quantum feature map, $x\xrightarrow{} \rho(x)$. In practice, $\rho_{\Phi}$ are pure states reducing the trace to $k(x', x'')= |\langle\Phi(x'')| \Phi(x')\rangle|^2 = |\bra{0} \mathcal{U}(x'')^{\dagger}\mathcal{U}(x')\ket{0}|^2$, allowing one to obtain an estimate of the kernel by sampling the probability of measuring $\ket{0}$ on the $\mathcal{U}(x'')^{\dagger}\mathcal{U}(x')\ket{0}$ state.

The Nystr\"{o}m method can be used for the matrix completion of positive semi-definite matrices. Using the subset of columns measured through the inner product of elements between $\mathcal{S}$ and $\mathcal{S}_L$ for an arbitrary node (dropping the the node label for brevity), which we define  as the $N\times L$ matrix $G := [W, B]^\top$ where $G_{ij} = k(x_i, x_j)$, $i\leq N$, $j\leq L$ and $W\in\mathbb{R}^{L\times L}$, $B\in\mathbb{R}^{L \times (N-L)}$, we complete the $N\times N$ kernel matrix by making the approximation $K \approx G W^{-1} G^\top$. Expanding, we have,
\begin{equation}
     K \approx \widehat{K} :=
    \begin{bmatrix}
    W      &  B \\
    B^\top      & B^\top W^{-1} B
    \end{bmatrix},
\end{equation}
\noindent where, in general, $W^{-1}$ is the Moore-Penrose gneralised inverse of matrix $W$. This becomes important in cases where $W$ is singular. Intuitively, the Nystr\"{o}m method utilises correlations between sampled columns of the kernel matrix to form a low-rank approximation of the full matrix. This means that the approximation suffers when the underlying matrix $K$ is close to full-rank.
However, the \textit{manifold hypothesis} \cite{10.5555/3086952} suggests that choosing $L<<N$ is not entirely futile as data generally does not  explore all degrees of freedom and often lie on a sub-manifold, where a low-rank approximation is reasonable.

The Nystr\"{o}m approximation does not change the positive semi-definiteness of the kernel matrix, $\widehat{K} \succeq 0$, and the SVM optimisation problem remains convex. Furthermore, the reproducing kernel Hilbert space (RKHS), induced by the adjusted kernel, $\hat{k}$, is a nonparametric class of functions of the form,
$
    \mathbb{H}_{\Phi} = \left\{ f_{\Phi} \big| f_{\Phi}(\cdot) = \sum_i^N \alpha_i'  \hat{k}_\Phi(\cdot, x_i) , \alpha_i' \in \mathbb{R}\right\}
$,
\noindent which holds as a result of the celebrated \textit{representer theorem} \cite{kimeldorf1971some, 10.1007/3-540-44581-1_27}. The specific set of $\{\alpha_i'\}_{i=1}^N$ for a given dataset is obtained through solving the SVM quadratic program. This constructs a split function that has the form,
\begin{equation}\label{eq:split_function}
    \mathcal{N}_{\Phi; \alpha} (x)=  \mathrm{sign}\left[ \sum_{i=1}^{N} \alpha_i \widetilde{y}_i \hat{k}_{\Phi}(x, x_i) \right],
\end{equation}
\noindent where $\alpha_i' = \alpha_i \widetilde{y}_i$, $\alpha_i \geq 0$ and $\widetilde{y}_i = F(y_i)$ with the function $F:\mathcal{C} \xrightarrow{} \{-1, 1\}$ mapping the original class labels, $y_i\in\mathcal{C}$ -- which in general can be from a set of many classes, i.e. $|\mathcal{C}|>2$, to a binary class problem. In Section \ref{section:results}, we provide numerical comparisons between two possible approaches of defining the function $F$, namely, (i) one-against-all (OAA), and (ii) even-split (ES) strategies -- see Appendix \ref{sec:multiclass_classification} for more details.

The construction of a split function of this form adopts a \textit{randomised node optimisation} (RNO) \cite{10.1007/978-1-4471-4929-3} approach to ensuring that the correlation between QDTs are kept to a minimum. The randomness is injected through the selection of landmark data points, as the optimised hyperplane will differ depending on the subset chosen. Furthermore, there exist possibilities to vary hyperparameters, $\Phi$ and $L$ both, down and across trees. This is depicted in Figure \ref{fig:qrf_diagram} with split functions notated with a separate tuple $(\Phi_i, L_i)$ for depths $i=1,...,D-1$ of the tree.
The specific kernel defined by an embedding $\Phi$ implies a unique Hilbert space of functions for which $k_{\Phi}$ is the reproducing kernel \cite{10.1090/s0002-9947-1950-0051437-7}. The two types of embeddings numerically explored are defined in Appendix \ref{section:embedding}. The QRF approach of employing distinct kernels at each depth of the tree, gives a more expressive enhancement to the QSVM method. A simple ensemble of independent QSVMs would require greater computational effort and would also lose the tree structure present in the QRF.

\section{Theoretical Results}\label{sec:theory}

QML literature has seen a push towards understanding the performance of models through the well established perspective of statistical learning theory. There has recently been a great deal of work determining the generalisation error for various quantum models \cite{Caro_2020, 10.1038/s41467-021-22539-9,https://doi.org/10.48550/arxiv.2103.03139, PRXQuantum.2.040321} that illuminate optimal strategies for constructing quantum learning algorithms. This is often achieved by bounding the generalisation error (or risk) of a model producing some hypothesis $h \in \mathbb{H}$, $\mathcal{R}(h) = \Pr_{(x,y)\sim \mathcal{D}} [h(x) \neq y]$, where $\mathcal{D}$ indicates the underlying distribution of the dataset. However, $\mathcal{D}$ is unknown and one simply has a set of samples $\mathcal{S}\sim \mathcal{D}^N$. We therefore define the empirical error of a particular hypothesis $h$ as the error over known points, $\widehat{\mathcal{R}}(h) = \frac{1}{N}\sum_{(x_i , y_i) \in \mathcal{S}} \mathbf{1}_{h(x_i)\neq y_i}$. This is often referred to as the training error associated with a model and is equivalent to $(1-\text{accuracy})$. Achieving a small $\widehat{\mathcal{R}}(h)$ however does not guarantee a model capable of predicting new instances. Hence providing an upper bound to the generalisation error $\mathcal{R}(h)$ is imperative in understanding the performance of a learning model.

The QDT proposed in this work can be shown to generalise well in the case where the margins at split nodes are large. This is guaranteed by the following Lemma which builds on the generalisation error of \textit{perceptron decision trees} \cite{10.1023/a:1007600130808}, with proof supplied in Appendix \ref{sec:error_bounds}.
\begin{lemma}\label{lemma:gen_err_qdt} (Generalisation error of Quantum Decision Trees)
Let $\mathbb{H}_J$ be the hypothesis set of all QDTs composed of $J$ split nodes. Suppose we have a QDT, $h\in \mathbb{H}_J$, with kernel matrix $K^{(i)}$ and labels $\{y^{(i)}_j\}_{j=1}^{N^{(i)}}$ given at the $i^{\text{th}}$ node. Given that $m$ instances are correctly classified, with high probability we can bound the generalisation error,
\begin{align}
 \mathcal{R}(h) \leq  \widetilde{\mathcal{O}}&\Bigg(\frac{1}{m} \Bigg[ J \log\left(4mJ^2\right) \nonumber \\
 &+\log(4m)^2 \sum_{i=1}^J \Big|\sum_{j,k}^{N^{(i)}} y^{(i)}_{j} y^{(i)}_{k}(K^{(i)+})_{jk} \Big| \Bigg]\Bigg),
\end{align}
where $K^{(i)}_{jk}=|\langle \Phi(x^{(i)}_j) | \Phi(x^{(i)}_k) \rangle|^2$ and $\{y_j^{(i)}\}_{j=1}^{N^{(i)}}$ are respectively the Gram matrix and binary labels given to the $i^{\text{th}}$ node in the tree, and $K^+$ denotes the Moore-Penrose generalised inverse of matrix $K$.
\end{lemma}
\noindent The term $s_K = \sum_{j,k}^{N^{(i)}} y^{(i)}_{j} y^{(i)}_{k}(K^{(i)+})_{jk}$ is referred to as the model complexity and appears in the generalisation error bounds for kernel methods \cite{10.1038/s41467-021-22539-9}. The term is both inversely related to the kernel target alignment \cite{Cristianini2006} measure that indicates how well a specific kernel function \textit{aligns} with a given dataset, as well as being an upper bound to the reciprocal of the geometric margin of the SVM. Therefore, Lemma \ref{lemma:gen_err_qdt} highlights the fact that, for a given training error, a QDT with split nodes exhibiting larger margins (and smaller model complexities) are more likely to generalise. Crucially, the generalisation bound is not dependent on the dimension of the feature space but rather the margins produced. This is in fact a common strategy for bounding kernel methods \cite{10.1007/s00362-019-01124-9}, as there are cases in which kernels represent inner products of vectors in infinite dimensional spaces \cite{10.7551/mitpress/4175.001.0001}. In the case of quantum kernels, bounds based on the Vapnik-Chervonenkis (VC) dimension would grow exponentially with the number of qubits. The result of Lemma \ref{lemma:gen_err_qdt} therefore suggests that large margins are analogous to working in a lower VC class. However, the disadvantage is that the bound is a function of the dataset and there are no \textit{a priori} guarantees that the problem exhibits large margins.

The Nystr\"{o}m approximated quantum kernel matrix used in the optimisation of the SVM is also a low rank approximation, with estimated elements that are subject to finite sampling error. This prompts an analysis on the error introduced -- stated in the following Lemma.
\begin{lemma}\label{lemma:k_k_tilde_main}
Given a set of data points $\mathcal{S} = \{x_i\}_{i=1}^N$ and a quantum kernel function $k$ with an associated Gram matrix $K_{ij} = k(x_i, x_j)$ for $i,j=1,...,N$, we choose a random subset $\mathcal{S}_L \subseteq \mathcal{S}$ of size $L\leq N$ which we define as $\mathcal{S}_L = \{x_i\}_{i=1}^L$ without loss of generality. Estimating each matrix element with $M$ Bernoulli trials such that $\widetilde{K}_{ij} = (1/M)\sum_{p}^M \widetilde{K}_{ij}^{(p)} $ with $\widetilde{K}_{ij}^{(p)} \sim \mathrm{Bernoulli}(k(x_i, x_j))$ for $i\leq N$, $j\leq L$, the error in the Nystr\"{o}m completed $N\times N$ matrix, $\widetilde{K}$, can be bounded with high probability,
\begin{equation}\label{eq:k_k_tilde_bound}
    || K - \widetilde{K} ||_2 \leq \widetilde{\mathcal{O}} \Bigg(\frac{NL}{M} + \frac{N}{\sqrt{L}}\Bigg),
\end{equation}
\noindent where $\widetilde{\mathcal{O}}$ hides the $\log$ terms.
\end{lemma}
\noindent The proof is given in Appendix \ref{sec:error_bounds}. This Lemma indicates that there are two competing interests when it comes to reducing the error with respect to the expected kernel. Increasing the number of landmark points $L$ reduces the number of elements estimated and hence the error introduced by finite sampling noise. On the other hand, the Nystr\"{o}m approximation becomes less accurate as expressed through the second term in Eq. \eqref{eq:k_k_tilde_bound}.
However, it is important to realise that the approximation error of the Nystr\"{o}m method was by design, in an attempt to weaken the effectiveness of the split function. To understand the effect of finite sampling error on the SVM model produced, we employ Lemma \ref{lemma:k_k_tilde_main} to show that we can bound error in the model output -- proof in \ref{sec:error_bounds}.

\begin{lemma}\label{lemma:pred_error_fin_samp}
Let $f(\cdot)=\sum_i \alpha_i k(\cdot, x_i)$ be the ideal Nystr\"{o}m approximated model and $\widetilde{f}(\cdot)=\sum_i \alpha_i' \widetilde{k}(\cdot, x_i)$ be the equivalent perturbed solution as a result of additive finite sampling error on kernel estimations. With high probability, we can bound the error as,
\begin{equation}
    |f(x)-\widetilde{f}(x)| \leq \mathcal{O}_{M}\left( \frac{N^{4/3}\sqrt{L}}{\sqrt{M}} \right)
\end{equation}
\noindent where $\mathcal{O}_M$ expresses the term hardest to suppress by $M$.
\end{lemma}
\noindent Lemma \ref{lemma:pred_error_fin_samp} indicates that $M\sim \mathcal{O}(N^3L)$ will suffice to suppress errors from sampling. This is in comparison to $M\sim \mathcal{O}(N^4)$ for without the Nystr\"{o}m approximation \cite{10.1038/s41567-021-01287-z}. However, note that this is to approximate the linear function $f$ for which the SVM has some robustness with only  $\mathrm{sign} [f(\cdot)]$ \eqref{eq:split_function} required.
A caveat to also consider is that it is not necessarily true that a smaller $|| K - \widetilde{K} ||_2$ or $|f(x)-\widetilde{f}(x)| $ will give greater model performance  \cite{sharp_anal_kern_matr_approx}.

To address the complexity of the QRF model, for an error of $\epsilon = \mathcal{O}(1/\sqrt{M})$ on kernel estimations, we can show that training requires $\mathcal{O}(TL(d-1)N\epsilon^{-2})$ samples and single instance prediction complexity of $\mathcal{O}(TL(d-1)\epsilon^{-2})$ circuit samples -- discussed further in Appendix \ref{section:complexity_appendix}. The main result here is that we are no longer required to estimate $\mathcal{O}(N^2)$ elements, as is the case for QSVMs. Though this is a profound reduction for larger datasets with $L<<N$, it should be noted that datasets may require $L=\mathcal{O}(N)$. Nonetheless, the number of estimations will never be greater than $N^2$, on the basis that kernel estimations are stored in memory across trees.

Finally, to show that the QRF contains hypotheses unlearnable by both classical learners and linear quantum models, we extend the concept class generated with the \textit{discrete logarithm problem} (DLP) in \cite{10.1038/s41567-021-01287-z}, from a single dimensional clustering problem to one in two dimensions. We construct concepts that separate the $2$D log-space (torus) into four regions that can not be differentiated by a single hyperplane. Hence, the class is unlearnable by linear QNNs and quantum kernel machines -- even with an appropriate DLP quantum feature map. Classical learners are also unable to learn such concepts due to the assumed hardness of the DLP problem. Since QDTs essentially create multiple hyerplanes in feature space, there exists $f\in\mathbb{H}_{\mathrm{QDT}}$ to emulate a concept in this class. Further details are presented in Appendix \ref{section:qrf_advantage}.

\section{Numerical Results \& Discussion}\label{section:results}

\begin{figure*}
    \centering
    \includegraphics[width=500pt]{Figures/tree_params_figure.png}
    \vspace{-0.15cm}
    \caption{\textbf{(a)} Data pre-preprocessing stages preceding model analysis. See Appendix \ref{section:data_and_relabel_process} for further information on the relabelling process. All numerical results in this Figure use the $\mathcal{D}_{\text{FM}}$ dataset with a $3:2$ training/testing split over $300$ randomly sampled instances over a quantum kernel with $\Phi_{\text{IQP}}$ embedding (Appendix \ref{section:data_and_relabel_process}). The embedding is over $7$ qubits, unless otherwise stated. Figures \textbf{(b)} to \textbf{(d)} are the analysis of NQKE split function hyperparameters with each point an average over $50$ trials. \textbf{(b)} has a fixed $C=10$, while \textbf{(c)} has $(M, L)=(2048, 10)$, and \textbf{(c)} has $(M,C)=(2048, 10)$ fixed. The Figures \textbf{(e)} through \textbf{(i)} are the results over the entire QRF model with each point being averaged across 40 independent trials.
    Figure \textbf{(e)} with fixed parameters $(M,T, C)=(2048,1, 10)$, indicates with large enough $L$, the tree structure has marginal benefit. This is largely a consequence of relabelling instances to be optimal for the quantum kernel. Figure \textbf{(f)} illustrates the effect of increasing the number of trees in the ensemble, for a range of partition sizes -- i.e. the number of randomly sampled instances given to each tree. This is for fixed parameters $(M, L, d, C) = (2048, 5, 4, 10)$. \textbf{(g)} Performance of the QRF on multi-class problems with the original data labels of $\mathcal{D}_{\text{FM}}$. The exact classes introduced for each multi-class problem is given in Appendix \ref{section:data_and_relabel_process}. These results have the following fixed parameters: $(T, M, L, C)=(1, 2048, 5, 10)$. Finally we have the correlation results of Figures \textbf{(h)} and \textbf{(i)}, with \textbf{(h)} displaying a heatmap of the correlation between trees, and \textbf{(i)} showing the average correlation between trees in a QRF model plotted against its accuracy on test instances. These have fixed parameters, $(T, M, L, d, C) = (10, 2048, 5, 4, 10)$.}
    \label{fig:tree_params_figure}
\end{figure*}

The broad structure of the QRF allows for models to be designed specifically for certain problems. This includes the selection of a set of embeddings at each level of the tree that distinguish data points based on different aspects of the real-world dataset.
However, this brings about the requirement of selecting hyperparameters that are not necessarily trivial to determine as optimal.
There is often an interplay between hyperparameters that can mutually influence their effectiveness.

In this Section we provide numerical results (with all results simulated with shot noise) of the performance of the model when hyperparameters are varied. Importantly, this is preceded by the classical pre-processing of data -- illustrated in Figure \ref{fig:tree_params_figure}a -- which involves performing \textit{principal component analysis} (PCA) to adjust the dimensionality of a dataset to an embedding on a particular number of qubits. This is followed by normalisation $x_i\in[0,\pi]$ and a relabelling step. This work explores three types of relabelling strategies: (i) utilising the original labels of the dataset, $y_i=y_i^\mathrm{orig}$, (ii) relabelling dataset so that the geometric difference between the classical \textit{radial basis function} kernel and the quantum kernel with respect to a specific embedding are maximised \cite{10.1038/s41467-021-22539-9}, $(y_i)_i = \mathscr{R}_{\Phi}^{\text{QK}}(\{(x_i, y_i^\mathrm{orig})\}_i)$, and, (iii) relabelling the dataset so that they form discontinuous regions in quantum feature space, $(y_i)_i = \mathscr{R}_{\Phi}^{\text{QRF}}(\{(x_i, y_i^\mathrm{orig})\}_i)$, where the label $\text{QRF}$ indicates that such a relabelling process was constructed to highlight the performance of the QRF over other models. Further elaboration of relabelling can be found in Appendix \ref{section:data_and_relabel_process}.

The effect of varying hyperparameters on the split function are shown in Figures \ref{fig:tree_params_figure}b-d.
The result of increasing the number of landmark points are evident in \ref{fig:tree_params_figure}b and \ref{fig:tree_params_figure}d, with greater $L$ allowing for a better Nystr\"{o}m approximation. However, the increase in $L$ highlights the importance of bounding the first term of Eq. \eqref{eq:k_k_tilde_bound} as it then requires a larger number of circuit samples to limit the detrimental effects of shot noise. Hence, results exhibit behaviour that were identified theoretically in Lemma \ref{lemma:k_k_tilde_main}. Figure \ref{fig:tree_params_figure}c indicates the appropriate scale of the regularisation parameter, with signs of slight overfitting with larger $C$ favouring accuracy over wider margins.

The overall performance of the QRF with respect to the hyperparameters of maximum depth $d$, number of trees $T$, number of landmark points $L$, and partition size are shown in Figures \ref{fig:tree_params_figure}e and \ref{fig:tree_params_figure}f. Increasing the depth of the tree is seen to compensate the approximation error introduced by selecting a smaller subset of landmark points. Furthermore, increasing the number of trees in the forest is seen to have a significant impact on the learner's performance, with the number of instances randomly sampled for each tree (partition size) only modestly increasing model accuracy. This suggests no obvious advantage using all available training instances for all trees. Instead, we observe that a bagging approach, with trees receiving subsets of points, is quite effective.

One of the most distinguishing characteristics of the QRF model is its ability to accommodate multi-class problems. In Figure \ref{fig:tree_params_figure}g we see the performance of a single QDT on multi-class problems with the class split types: OAA and ES. The latter performs far better, with the QRF able to achieve similar accuracy on the four-class problem as with the binary problem. Crucially, this is with no change to the algorithm and further complexity. The leaf nodes simply output a probability distribution over $|\mathcal{C}|$ classes.
Though we naturally expect the required depth of the tree to increase with $|\mathcal{C}|$, considering the exponential growth of the number of leaf nodes -- with respect to the depth of the tree -- we conjecture a required sub-linear growth in depth, $d=\mathcal{O}(\mathrm{poly}(\log |\mathcal{C}|))$.

Ensemble machine learning methods are only useful in cases where its constituent weak classifiers are not highly correlated. Though we want individual learners to be correlated to the expected class labels, this must be balanced with the ineffectiveness of highly correlated classifiers. Figure \ref{fig:tree_params_figure}h and \ref{fig:tree_params_figure}i illustrate this by observing the \textit{Spearman's rank correlation} \cite{doi:https://doi.org/10.1002/9781118445112.stat05964} between trees.

The analysis of model hyperparameters now allows us to observe the QRF's performance in comparison to other machine learning models. Figure \ref{fig:model_comparison} shows accuracy of four models: the QRF, the SVM with a quantum kernel (QSVM), the classical \textit{random forest} (CRF) and the SVM with a classical RBF kernel. The comparison is given for embeddings $\Phi_{\text{IQP}}$ and $\Phi_{\text{Eff}}$ (that apply to the two quantum models, see Appendix \ref{section:data_and_relabel_process} for definitions) and the following datasets: $\mathcal{D}_{\text{Io}}$, $\mathcal{D}_{\text{H}}$, $\mathcal{D}_{\text{BC}}$ and $\mathcal{D}_{\text{FM}}$. Further information about these embeddings and datasets can be found in Appendix \ref{section:qrf_construction}. We see that the QRF is able to out-perform the QSVM in all cases where the dataset has been relabelled. This is surprising as $\mathscr{R}_\Phi^{\text{QK}}$ is specifically designed for a single quantum kernel. In these cases the QRF also out-performs the other classical learners. The results on the original class labels are not surprising. In most cases the classical learners have a greater performance, which conforms with results that address the requirement of an embedding with the correct inductive bias \cite{inductive_bias_of_qk}. In other words we require a kernel function that induces an \textit{effective} metric on the dataset. Hence we reach a limitation of quantum kernel methods: there exist no general $\Phi$ for which the kernel function produces a favourable higher dimensional map for classical datasets. We discuss this and other limitations further in Appendix \ref{sec:limitations_q_kernels}. However, it is for this reason we do not attempt to solve the predicament of obtaining useful embeddings, instead emphasising providing proof of concepts through relabeling datasets.

\begin{figure}
    \centering
    \includegraphics[width=240pt]{Figures/model_comparison.pdf}
    \caption{The comparison of quantum and classical classifiers on four datasets: $\mathcal{D}_{\text{Io}}$, $\mathcal{D}_{\text{H}}$, $\mathcal{D}_{\text{BC}}$ and $\mathcal{D}_{\text{FM}}$. The numeric simulations were carried out with two embeddings: $\Phi_{\text{IQP}}$ and $\Phi_{\text{Eff}}$; and the three relabelling strategies: original label, QK relabelled ($\mathscr{R}_\Phi^{\text{QK}}$) and QRF relabelled ($\mathscr{R}_\Phi^{\text{QRF}}$). Note: the two embeddings affect only the quantum models. QRF models have the following hyperparameters, $(T, M, L,d) = (5, 2048, 10, 4)$ for the first three datasets and we increase to $T=10$ for the larger $\mathcal{D}_{\text{FM}}$ dataset. The CRF has identical $(T,d)$ to the QRF. Due to the sensitivity of the models to $C$, the $\Phi_{\text{IQP}}$ model results have $C=10$ for QRF, QSVM and SVM classifiers, while $C=100$ is set for $\Phi_{\text{Eff}}$. Kernel estimations are stored throughout the depth and across the ensemble of trees. Hence, $\sigma$ states the ratio between kernel estimations required by QRF against QSVM.
    Note: $\mathscr{R}_\Phi^{\text{QRF}}$ generates new labels for each trial by randomly selecting new pivot points. This results in the larger error bars shown. The datasets, embeddings and relabelling process are elaborated further in Appendix \ref{section:data_and_relabel_process}.}
    \label{fig:model_comparison}
\end{figure}

\section{Conclusion}\label{section:conclusion}

The well-known \textit{no-free-lunch} theorem in ML \cite{shalev-shwartz_ben-david_2014} refers to the inconvenient fact that there exists no single learning model that is optimal for all datasets. The commonly explored QNN and QSVM methods are linear models in quantum feature space, with the form $f(x)=\Tr[\rho(x)O_\theta]$.
In this work, we show that the introduction of a learner with a decision tree structure generates a unique quantum learning model that is able to produce hypotheses that can provide a clear advantage due to its non-linearity. Furthermore, the tree structure of the QDT can compensate for a quantum embedding that is unable to perfectly separate instances in feature space. This is demonstrated both, theoretically with the construction of a DLP-based concept class, and numerically with the QRF that is observed to outperform QSVMs in most cases with identical embeddings.

In addition, unlike many other quantum models, the QRF can produce a probabilistic output of the estimated likelihood of class labels for a particular instance. This further allows for multi-class problems to be addressed with no structural change in the model. Employing a low-rank Nystr\"{o}m approximation (as we know generalisation suffers with large rank \cite{10.1038/s41467-021-22539-9}) to reduce the effectiveness of split functions, reduces the complexity from $\mathcal{O}(N^4)$ to $\mathcal{O}(N^3L)$ shots per estimation required to limit the errors introduced by finite sampling on the linear function representing the hyperplane. However, one should note that the split function takes the sign of the function and therefore such a large number of samples are not required in practice. Moreover, as QDTs are intended to be weak learners there are potential regularisation effects due to such noise \cite{noisy_quantum_kernel_machines}. The Nystr\"{o}m method also allows for the complexity of the number of kernel estimations of the entire model to not grow quadratically with the number of training instances -- which is a significant limitation of QSVMs. Possible improvements to NQKE can observe alternate methods to select landmark points. Examples include the selection of columns by incomplete Cholesky decomposition with column pivoting \cite{10.5555/645529.657980}.

The required model hyperparameters for any given data set and its relationship with the number of qubits used for the quantum split function, is not a trivial question and there are various avenues in which the QRF can be optimised. Though we have numerically explored varying a range of hyperparameters, the QRF model can easily be altered to respond to the specific problem at hand. This includes the use of multiple kernel embeddings down the QDT, or further reducing the complexity of the model for larger datasets by taking only a subset of instances for training at each node -- not only a subset of points to act as landmarks. Furthermore, there is the possibility of utilising techniques such as boosting \cite{boosting_qsvm} and other pruning methods, such as CART's cost complexity pruning algorithm \cite{li1984classification}, to improve the performance of the QRF model.

\textit{The QRF model can be accessed at the following code repository \cite{github}. Data is available upon reasonable request.}

\begin{acknowledgments}

MS is supported by the Australian Government Research Training Program (RTP) Scholarship. CDH is partially supported by the Laby Foundation research grant. The authors would like to thank Mario Kieburg for useful discussions. We acknowledge the support provided by the University of Melbourne through the establishment of an IBM Network Quantum Hub.
The large number of simulations required for this work were made feasible through access to the University of Melbourne's High Performance Computer, \textit{Spartan} \cite{spartan}.

MS conceived the QRF approach and carried out the mathematical development, computations and analysis under the supervision of CH and LH. MS wrote the paper with input from all authors.

\end{acknowledgments}