Learning the N-input Parity Function with a Single Qubit and Single Measurement Sampling

1. Introduction: The Parity Problem

Representing the XOR function is a fundamental problem for neural networks (NNs), as it involves solving a nonlinearly separable classification problem. Unlike other problems represented by logical functions with two-bit inputs, the XOR problem requires the use of two layers of neurons [1]. The parity function extends the XOR function from 2 to N inputs, forming a proportionally challenging binary classification problem for $2^N$ inputs, an illustration of which in the 3-dimensional space ($N=3$) is depicted in Figure 1. According to the universal approximation theorem, any function, including the parity function, can be represented by a deep NN. However, it is always preferable to employ an NN architecture with minimal cost in terms of the number of synapses, neurons, and layers. Another key consideration is the morphology of the solution landscape, as well as the representational capacity and extrapolation capabilities of the chosen NN model. An optimal NN architecture using perceptrons with one hidden layer of size N is described in [2,3], while a modular NN architecture is presented in [3], and a single product unit approach is discussed in [4].

Variational Quantum Circuits (VQCs) [5,6,7] combine the properties of quantum circuits and classical methods to solve classification tasks in a hybrid way. In particular, quantum circuits perform the data processing, which encompasses data encoding, label extraction, etc., whereas optimization is handled by classical –typically local optimization– algorithms. The data flow follows a sequential combination of the processes, meaning that data collected via measurements on the output quantum state of the VQC are subsequently exploited to perform the optimization task. The purpose of this hybrid formation is to optimally harvest the best of both worlds in an attempt to create new efficient learning methods. However, quantum measurement results are intrinsically stochastic, presenting a challenge in developing suitable optimization algorithms [8].

In this work, our aim is to employ an optimized VQC model in place of an NN and demonstrate that it can be efficiently trained to solve the N-bit parity problem. It is widely known all boolean logical functions can be represented by a quantum circuit[9]. Considering this, one can easily construct a quantum Boolean parity function using N qubits as input and an ancilla qubit to register the output, as shown in Figure 2. While the problem of computing the parity of an oracle function has been explored in literature, it has been shown [10] that at least $N/2$ oracle queries are required, offering no quantum advantage. Finally, another class of parity problems has been used to evaluate the effectiveness of different VQC models in [11], further highlighting the importance of finding an efficient way of solving the parity problem.

The remainder of the manuscript is divided so as to showcase two main concepts; the first concerns the presentation of the quantum variational model for solving the parity problem; the second follows the methodology for training this model. More specifically, in Section 2, we propose a variational Single-Qubit Classifier and investigate the solution landscape for the parity problem. In Section 3, we propose and test a new stochastic gradient method, to the best of our knowledge, for which batches are represented by density matrices. This method can be transformed into a doubly stochastic approach if single stochastic measurement outcomes are considered and averaged over the batches. We present results from numerical tests on both methods, demonstrating the effectiveness of the approach.

2. A Single-Qubit Classifier For The Parity Problem

Given a series of N bits, $\overrightarrow{x}=x_1,x_2,\dots ,x_N$, where $x_i\in \left\{0,1\right\}$, the parity function determines whether the value "1" appears an even or odd number of times in the series. The parity problem can, thus, be considered a classification problem, where the $n=2^N$ data are divided into two equally populated classes: class A (even parity) and class B (odd parity), as shown in Figure 1. This problem can also be extended to real numbers, i.e., noisy data $\overrightarrow{\tilde{x}}$, generated by adding Gaussian noise $\mathcal{N}(\mu ,{\sigma }^2)$ to the integer values of $\overrightarrow{x}$.

Following the logic of [12], where a general model of a variational qudit classifier was introduced, we specify the following model for the parity problem:

$$exp\left[-i\sum _{k=1}^3{s}_k\pi {\widehat{\sigma }}_k\right]exp\left[-i\left(\sum _{j=1}^N{w}_j{x}_j\right)\pi {\widehat{\sigma }}_1\right]\left|0\right.〉.$$

(1)

We will refer to this model as the Single-Qubit Classifier (SQC). In Equation (1), ${\widehat{\sigma }}_1,{\widehat{\sigma }}_2,{\widehat{\sigma }}_3$ represent the Pauli operators, $|0\rangle$ is the eigenstate of ${\widehat{\sigma }}_3$ with eigenvalue 1, and $\overrightarrow{s}=(s_1,s_2,s_3)$, $\overrightarrow{w}=(w_1,w_2,..,w_N)$, are the variational parameters, which are optimized during training. The training procedure updates the variational parameters given data points $\overrightarrow{x_i}=(x_{i1},x_{i2},...,x_{iN})$, with each data point labeled as $y_j=1$ for even or $y_j=-1$ for odd parity. The classes are, hence, represented by labels that can be directly matched to the mean measurement outcome of operator ${\widehat{\sigma }}_3$, ${\langle {\widehat{\sigma }}_3\rangle }_j$, where positive is classified as even parity and negative is classified as odd parity.

2.1. Landscape Of Solution

The SQC in Equation (1) is inspired by the quantum digital circuit shown in Figure 2, the relation between the two can be directly derived. For noiseless data, a periodic (global) solution emerges, allowing the SQC to classify all inputs with perfect accuracy when ${\overline{w}}_j=1/2+l$ and $s_k=0$ for all j and k, with $l\in \mathbb{Z}$. In this solution, the SQC reproduces the conditional evolution of the ancillary qubit in Figure 2: an X gate is applied whenever $x_j=1$. An even number of gate applications results in the final state $|0\rangle$, while an odd number of applications leads to the state $|1\rangle$. A measurement of the observable ${\widehat{\sigma }}_3$ at the end of the variational circuit can perfectly distinguish the two cases, resulting in the accurate classification of the input.

Even if the parameters converge to a solution with some deviation $\delta w$ around the global solution $\overline{w}$, the model classifies noiseless data successfully. Correspondingly, if the parameters are fixed at $\overline{w}$, the classifier can tolerate noise of the same order in the data. We observe that, for $N=2$, the classifier can tolerate a maximum deviation $\delta w_{max}\sim 1/8$, or equivalently ${\sigma }_{max}\sim 1/24$. For N inputs, the total noise in the exponent of Equation (1) increases by $\sqrt{N}\sigma$, meaning that the tolerable noise for each individual input scales with ${\sigma }_{max}\sim 1/24\sqrt{N}$. Finally, the second rotation of the SQC’s state – represented by the second term of Equation (1) – is added, so that the model can tolerate recurrent additive errors in the form of $\mathcal{N}(\mu \ne 0,{\sigma }^2)$ in the input data. In the remainder of this work, we will ignore this term and set $\overrightarrow{s}=\mathbf{0}$ in Equation (1).

We initially employ the square loss function $\mathcal{L}\left(\overrightarrow{w}\right)={\sum }_{j=1}^N{(y_j-{\langle {\widehat{\sigma }}_3\rangle }_j)}^2$ for evaluating the result of the learning procedure of the VQC, however, we observe that no local minima appear in this case. The solution landscape for $N=2$ can be seen in Figure 3 (a) over two periods. Since illustrations of the landscape’s morphology can only be extracted as cross-sections when $N>2$, we plot an example for $N=2,6,$ and 10, where the cross-section is applied along the principal hyper-diagonal, namely, where $x_1=x_2=\dots =x_N=z$, shown in Figure 3 (b). We notice that, as N increases, the number of local minima along the diagonal increases moderately. Nonetheless, the existence of local minima persists if $z\to 1-z$ for some inputs. This suggests that the number of minima within a period increases at least linearly with N, while the global solution remains unique.

On the other hand, our numerical investigation shows that input sets which include a subset of points whose coordinates coincide with that of the global solution and, simultaneously, the rest of the data points have coordinates of local minima, also correspond to local minima of the loss function. Reaching the aforementioned minima will also yield a perfect classification of the input. Based on the analysis of numerical results in the following section, we conclude that both local minima that achieve accuracy 1, meaning that all points are correctly classified, and those that do not, increase in number proportionally as N grows. A final observation on Figure 3 (b), supported by our numerical investigations, is that as N increases, the area in the parametric space covered by barren plateaus also increases. This makes the approach to a local minimum more time-consuming or, in some cases, impossible.

3. Learning Methods

The versatile landscape of the parametric space, combined with the fact that there is at least one known solution, renders the parity problem an attractive test-bed for investigating different optimization techniques for the training of SQV. We will, first, provide a description of the common steps involved in training a variational quantum classifier, followed by an introduction to the basic optimization technique for adjusting the weights, namely the Gradient Descent (GD) method. It is essential to gain a deeper understanding of the solution landscape through the application of GD before more advanced optimization techniques can be explored.

Let us restrict the discussion to the case of binary classification tasks, where the data ${\overrightarrow{x}}_j$ are to be classified into two classes, A and B, consisting of $N_A$ and $N_B$ points respectively, where $N_A+N_B=n$. We denote the output circuit state by $|\psi (\overrightarrow{w},{\overrightarrow{x}}_j)\rangle$ and assume that a measurement on a single observable $\widehat{G}$ (specifically, ${\widehat{\sigma }}_3$ in our case), which has two eigenvalues $\pm 1$, is performed on the quantum state. For fixed values of the parameters $\overrightarrow{w}$, the standard procedure involves running the quantum circuit multiple times to estimate

$$\begin{array}{cc}\hfill g(\overrightarrow{w},{\overrightarrow{x}}_j)=& \langle \psi (\overrightarrow{w},{\overrightarrow{x}}_j)|\widehat{G}|\psi (\overrightarrow{w},{\overrightarrow{x}}_j)\rangle \\ \hfill \approx {\tilde{g}}_M(\overrightarrow{w},{\overrightarrow{x}}_j)=& {\displaystyle \frac{1}{M}}\sum _{m=1}^M{G}_m(\overrightarrow{w},{\overrightarrow{x}}_j)\end{array}$$

(2)

where $G_m(\overrightarrow{w},{\overrightarrow{x}}_j)\in \left\{-1,+1\right\}$ represents the stochastic outcome of the m-th (independent) experiment. For $M\gg 1$, we have ${\tilde{g}}_M(\overrightarrow{w},{\overrightarrow{x}}_j)\to g(\overrightarrow{w},{\overrightarrow{x}}_j)$.

The step of estimating $g(\overrightarrow{w},{\overrightarrow{x}}_j)$ via the quantum circuit and the corresponding measurements is followed by a classical procedure to update the parameters $\overrightarrow{w}$ using a classical optimization method. The first step in this procedure is to construct a loss function $\mathcal{L}\left(\overrightarrow{w}\right)$ based on $g(\overrightarrow{w},{\overrightarrow{x}}_j)$ and the information about the classes. Here, we use the negative log-likelihood loss, which is commonly employed to learn the weights in NNs. If the classes A and B – corresponding to even and odd parity in our case study, respectively – are labeled as $y=\pm 1$, the probability that the training data ${\overrightarrow{x}}_i$ belongs to the class y can be evaluated as follows:

$$P\left(y\right|{\overrightarrow{x}}_j)={\displaystyle \frac{yg(\overrightarrow{w},{\overrightarrow{x}}_j)+1}{2}}.$$

(3)

By defining the class label $y_j$ for each input ${\overrightarrow{x}}_j$, the loss function can be constructed by summing the negative log-likelihood of the probabilities as follows:

$$\mathcal{L}\left(\overrightarrow{w}\right)=-{\displaystyle \frac{1}{n}}\sum _{j=1}^{n}log\left(P\left(y_j|{\overrightarrow{x}}_j\right)\right)$$

(4)

or,

$$\mathcal{L}\left(\overrightarrow{w}\right)=-{\displaystyle \frac{1}{N_A}}\sum _{j\in A}log\left(P\left(y=1|{\overrightarrow{x}}_j^A\right)\right)-{\displaystyle \frac{1}{N_B}}\sum _{j\in B}log\left(P\left(y=-1|{\overrightarrow{x}}_j^B\right)\right)$$

(5)

Using GD, the weights are updated according to

$${\overrightarrow{w}}_{t+1}={\overrightarrow{w}}_t-\alpha \overrightarrow{\nabla }\mathcal{L}\left({\overrightarrow{w}}_t\right)$$

(6)

where $\alpha$ is the learning rate. The numerical differentiation that produces the gradient $\mathcal{L}\left(\overrightarrow{w}\right)$ of the loss function requires at least $2N$ evaluations of g, defined in Equation (2), to be estimated by the quantum circuit, where N is the dimension of each input data vector. The analytical formulas [13,14] of the partial derivatives of the gradient vector can potentially relate these derivatives to the mean values of other observables in the quantum state, offering an improvement in the precision of the optimization procedure. In our case, where the model is relatively simple, the analytic formula is straightforward and involves the mean value of only one observable:

$${\displaystyle \frac{\partial g(\overrightarrow{w},{\overrightarrow{x}}_i)}{\partial w_j}}=2x_j^{\left(i\right)}\pi \langle \psi (\overrightarrow{w},{\overrightarrow{x}}_i)|{\widehat{\sigma }}_2|\psi (\overrightarrow{w},\overrightarrow{x})\rangle =2x_j^{\left(i\right)}\pi h(\overrightarrow{w},{\overrightarrow{x}}_i)$$

(7)

where $x_j^{\left(i\right)}$ denotes the jth component of the vector ${\overrightarrow{x}}_i$. Using the chain rule and Equation (7), one can accurately evaluate the elements of the gradient ${\displaystyle \frac{\partial \mathcal{L}\left(\overrightarrow{w}\right)}{\partial w_j}}$.

For the parity problem with $N=6,\dots ,10$-bit input, we apply the GD method to optimize the parameters of the SQC, as given by Equation (1). We use a constant learning rate $\alpha$ and a fixed number of epochs, and we employ the analytical formula in Equation (7) for the gradient. We report two important metrics concerning the classification results; the average accuracy $A=(\#\mathrm{inputs}\mathrm{correctly}\mathrm{classified})/\left(\#\mathrm{all}\mathrm{inputs}\right)$, listed in Table 1, as well as the percentage of times that accuracy reached 1, i.e., where all inputs were correctly classified after training. The metrics are reported over 50 runs where the weight vector is randomly initialized for each run.

The results of the GD optimization, shown in Table 1, as well as the evolution of the observed loss during optimization, lead to the following conclusions. The difficulty of the problem increases rapidly with the input dimension N, as expected, highlighting the need for more fine-tuned optimization methods. For $N=6,7,8$, the few occurrences of unsuccessful training are mainly due to the optimization method being trapped in the small neighborhood around local minima. However, for $N\ge 9$, the presence of barren plateaus becomes equally problematic and the loss function appears to oscillate randomly over these regions.

3.1. The ESGD Method: Towards An Optimization With Single Measurement Sampling

Inspired by recent studies [8,15,16] on exploiting the intrinsic stochasticity of quantum measurements in gradient optimization methods and, at the same time, motivated by the need to reduce quantum resources, we explore whether the SQC can be trained using just one measurement (one repetition of the quantum circuit) for a very rough estimate of $g(\overrightarrow{w},\overrightarrow{x})$ and ${\displaystyle \frac{\partial g(\overrightarrow{w},\overrightarrow{x})}{\partial w_j}}$. We begin by justifying this approach through the introduction of the ESGD method and then test the use of ${\tilde{g}}_{M=1}(\overrightarrow{w},{\overrightarrow{x}}_j)=\pm 1$ in the loss function. It is clear that some averaging must be involved to induce smoothness in the loss function, and intuition suggests that this average should be taken over different input data ${\overrightarrow{x}}_j$. In the following, we present a logical procedure to justify this averaging, while keeping the negative log-likelihood as the loss function, contrary to previous works.

We will now provide a brief description of the Stochastic Gradient Descent (SGD) optimization method. With SGD, n data points are randomly split into k batches, each batch containing $n_k=n/k$ points. During an epoch, the loss function, in our case given by Equation (8), and its derivatives are estimated k times.

$${\mathcal{L}}^{\left(k\right)}\left(\overrightarrow{w}\right)=-{\displaystyle \frac{1}{n_k}}\sum _{j=(k-1)n_k+1}^{k{n}_k+1}log\left(P\left(y_j|{\overrightarrow{x}}_j\right)\right)$$

(8)

Focusing, now, on the specific conditions of the parity problem, where $N_1=N_2=n/2$, we introduce a “balanced” cost functional with an equal number of points ($n_k/2$) from each class ${\overrightarrow{x}}_j^A$ and ${\overrightarrow{x}}_j^B$ in each batch:

$$\begin{array}{cc}\hfill {\mathcal{L}}^{\left(k\right)}\left(\overrightarrow{w}\right)& =-{\displaystyle \frac{2}{n_k}}\sum _{j=(k-1)n_k/2+1}^{k{n}_k/2}log\left(P\left(y=1|{\overrightarrow{x}}_j^A\right)\right)\\ & -{\displaystyle \frac{2}{n_k}}\sum _{j=(k-1)n_k/2+1}^{k{n}_k/2}log\left(P\left(y=-1|{\overrightarrow{x}}_j^B\right)\right).\end{array}$$

(9)

However, quantum mechanics offers an additional perspective, which has been partially explored in a different context in [17]. Let us construct $2k$ density matrices from the batches of states as follows:

$$\begin{array}{cc}\hfill {\widehat{\rho }}_{+1}^k& ={\displaystyle \frac{2}{n_k}}\sum _{j=(k-1)n_k/2+1}^{k{n}_k/2}|\psi (\overrightarrow{w},{\overrightarrow{x}}_j^A)\rangle \langle \psi (\overrightarrow{w},{\overrightarrow{x}}_j^A)|\end{array}$$

(10)

$$\begin{array}{cc}\hfill {\widehat{\rho }}_{-1}^k& ={\displaystyle \frac{2}{n_k}}\sum _{j=(k-1)n_k/2+1}^{k{n}_k/2}|\psi (\overrightarrow{w},{\overrightarrow{x}}_j^B)\rangle \langle \psi (\overrightarrow{w},{\overrightarrow{x}}_j^B)|.\end{array}$$

(11)

For the kth batch one can estimate

$$F(\overrightarrow{w},{\overrightarrow{x}}_{j\in k}^{A,B})=Tr\left({\widehat{\rho }}_{\pm 1}^k\widehat{G}\right)={\displaystyle \frac{2}{n_k}}\sum _{j=(k-1)n_k/2+1}^{k{n}_k/2}g(\overrightarrow{w},{\overrightarrow{x}}_j^{A,B}),$$

(12)

and then define $P_k(\pm 1|{\overrightarrow{x}}_j)={\displaystyle \frac{\pm 1F(\overrightarrow{w},{\overrightarrow{x}}_{j\in k}^{A,B})+1}{2}}$ and finally the loss function

$$\tilde{\mathcal{L}}{\left(\overrightarrow{w}\right)}^{\left(k\right)}=-log\left(P_k\left(y=1|{\overrightarrow{x}}_j^A\right)\right)-log\left(P_k\left(y=-1|{\overrightarrow{x}}_j^B\right)\right).$$

(13)

This method, which we shall refer to as Ensemble Stochastic Gradient Descent (ESGD), has a geometric interpretation. If a Bloch vector is assigned to each of the density matrices from Equations (10)-(), then this gradient descent method guides the ensemble of Bloch vectors associated with each class to point in opposite directions on the Bloch sphere ($\pm \widehat{z}$). The ESGD method is derived on solid foundations, since it is closely related to SGD, as described in the next. The first requires at least $n_k=2$ input points to populate each batch, since it includes the evaluation of Equations (10) and (). Considering, hence, the base case $n_k=2$, we have that

$$\begin{array}{cc}\hfill {\widehat{\rho }}_{+1}^k& =|\psi (\overrightarrow{w},{\overrightarrow{x}}_j^A)\rangle \langle \psi (\overrightarrow{w},{\overrightarrow{x}}_j^A)|\end{array}$$

(14)

$$\begin{array}{cc}\hfill {\widehat{\rho }}_{-1}^k& =|\psi (\overrightarrow{w},{\overrightarrow{x}}_j^B)\rangle \langle \psi (\overrightarrow{w},{\overrightarrow{x}}_j^B)|\end{array}$$

(15)

and the label assignment becomes

$$F(\overrightarrow{w},{\overrightarrow{x}}_{j\in k}^{A,B})=Tr\left({\widehat{\rho }}_{\pm 1}^k\widehat{G}\right)=g(\overrightarrow{w},{\overrightarrow{x}}_j^{A,B}),$$

(16)

reducing to the SGD base case. Finally, as with the standard GD, the components of the gradient can be calculated analytically using Equations (13) and (7).

Now, we will proceed to describe a modification of ESGD, according to which we consider a single measurement per input ${\overrightarrow{x}}_j$ for calculating either numerically, or analytically using the parametric shift rule [14], the gradient. We will refer to this as Doubly Stochastic Gradient Descent (DSGD), since it incorporates both batch processing and the stochasticity from single measurements. The convergence of DSGD to ESGD becomes clearer with the following formula, as the number of measurements M increases:

$$\begin{array}{cc}& F(\overrightarrow{w},{\overrightarrow{x}}_{j\in k}^{A,B})=Tr\left({\widehat{\rho }}_{\pm 1}^k\widehat{G}\right)\approx {\tilde{F}}_M(\overrightarrow{w},{\overrightarrow{x}}_{j\in k}^{A,B})\\ & ={\displaystyle \frac{2}{M{n}_k}}{\sum }_{j=(k-1)n_k/2+1}^{k{n}_k/2}{\sum }_{m=1}^M{G}_m(\overrightarrow{w},{\overrightarrow{x}}_{j\in k}^{A,B}).\end{array}$$

(17)

One may modify ESGD using any finite M to estimate ${\tilde{F}}_M(\overrightarrow{w},{\overrightarrow{x}}_{j\in k}^{A,B})$, then ${\tilde{P}}_k(\pm 1|{\overrightarrow{x}}_j)=\left(\pm 1{\tilde{F}}_M(\overrightarrow{w},{\overrightarrow{x}}_{j\in k}^{A,B})+1\right)/2$ and finally the loss function in Equation (13). In this work, however, we exhibit the power of the method for $M=1$.

We employ the SGD, ESGD, and DSGD optimization methods and conduct classification experiments on instances of the parity problem for an input dimension of $N=10$. We chose this dimension to illustrate the results, as it marks the point where GD starts to fail and because it allows for varying batch sizes. Our results, based on 25 experiments with random weight vector initialization, are presented in Table 2, Table 3 and Table 4. From these tables, it is evident that both ESGD and DSGD outperform SGD for all batch sizes. Specifically, the mean accuracy of ESGD and its stochastic counterpart, DSGD, increases as the batch size grows, whereas SGD experiences significant performance degradation, as shown in Figure 4. SGD eventually converges to the performance of GD (see Table 1 for $N=10$). The close relationship between SGD and ESGD for $n_k=2$ is also apparent, with both showing poor performance due to the continuous, uncorrelated treatment of the data, which fails to collectively guide the weight vector toward a solution. In contrast, ESGD consistently achieves high accuracy, close to 100% for any $n_k\ge 32$. DSGD, being stochastic, does not easily achieve 100% accuracy, as it tends to oscillate even once average convergence is reached.

Some important technical details about the implementation of the methods are as follows: The random weight vector initialization differs for each point in Figure 4. The performance of the methods is naturally influenced by the choice of hyperparameters, such as the learning rate and the number of epochs. For all three methods and each batch size, we optimized the learning rate to achieve the best outcomes. The number of function evaluations is kept constant at 1024, meaning the number of epochs increases with $n_k$. Note that the $N=2,4$ cases are not included in Table 4 due to the fact that the doubly stochastic derivation introduces excessive noise, preventing the model parameters from converging to a solution vector. Finally, in the DSGD method, to achieve smoother evolution, we average over the last four steps when evaluating the gradient vector.

4. Conclusions

This work presents two key contributions that advance both the theoretical understanding and practical application of quantum machine learning techniques. First, we demonstrate that the N-bit parity problem can be efficiently solved using a single qubit, highlighting the remarkable potential of qubits as powerful learning units in quantum circuits. This result underscores the fact that any Boolean logical function, including parity, can be implemented within a quantum framework. Our preliminary explorations of other logical functions, such as AND and OR, suggest that this approach may not be universally applicable to all N-input Boolean functions, thus motivating further exploration into the use of qudits. This opens an exciting avenue for future work, where the exploration of quantum systems beyond qubits can be further investigated.

The second significant contribution of this work is the introduction of a novel optimization method, ESGD, which utilizes batches represented by density matrices. Our numerical investigations show that ESGD outperforms traditional SGD in terms of accuracy and efficiency, particularly for the parity problem. Furthermore, ESGD maintains its effectiveness even when reduced to the extreme scenario of using a single measurement per input, a result that highlights its robustness in noisy environments. While our findings are promising, further validation across different problems and for larger quantum systems is required to fully assess the broader applicability and impact of ESGD.

In conclusion, this work makes substantial progress in applying quantum circuits to machine learning tasks, demonstrating the potential of quantum-based optimization methods and offering a new direction for future research in quantum-enhanced learning algorithms.