Quantum algorithms for feedforward neural networks

TL;DR

This paper introduces quantum algorithms for training and evaluating feedforward neural networks, achieving quadratic speedups over classical methods and potentially enhancing network regularization.

Contribution

It presents novel quantum algorithms for neural network training and evaluation, leveraging quantum subroutines and memory for efficiency and robustness.

Findings

Quadratic speedup in training and evaluation times for large networks.

Potential intrinsic resistance to overfitting due to quantum algorithm properties.

Quantum-inspired classical algorithms with similar dimension dependence but higher overhead.

Abstract

Quantum machine learning has the potential for broad industrial applications, and the development of quantum algorithms for improving the performance of neural networks is of particular interest given the central role they play in machine learning today. In this paper we present quantum algorithms for training and evaluating feedforward neural networks based on the canonical classical feedforward and backpropagation algorithms. Our algorithms rely on an efficient quantum subroutine for approximating the inner products between vectors in a robust way, and on implicitly storing large intermediate values in quantum random access memory for fast retrieval at later stages. The running times of our algorithms can be quadratically faster in the size of the network than their standard classical counterparts since they depend linearly on the number of neurons in the network, as opposed to the number of connections between neurons as in the classical case. This makes our algorithms suited for large-scale, highly-connected networks where the number of edges in the network dominates the classical algorithmic running time. Furthermore, networks trained by our quantum algorithm may have an intrinsic resilience to overfitting, as the algorithm naturally mimics the effects of classical techniques such as drop-out used to regularize networks. Our algorithms can also be used as the basis for new quantum-inspired classical algorithms which have the same dependence on the network dimensions as their quantum counterparts, but with quadratic overhead in other parameters that makes them relatively impractical.

Tables

Table 1. Table 1: ( ϵ , γ ) italic-ϵ 𝛾 (\epsilon,\gamma) -feedforward neural network accuracy (percentage of correctly labelled test points) on the MNIST handwriting data set, for various values of ϵ italic-ϵ \epsilon , and γ = 0.05 𝛾 0.05 \gamma=0.05 . A network with L = 4 𝐿 4 L=4 layers and dimensions [ n 1 , n 2 , n 3 , n 4 ] = [ 784 , 100 , 30 , 10 ] subscript 𝑛 1 subscript 𝑛 2 subscript 𝑛 3 subscript 𝑛 4 784 100 30 10 [n_{1},n_{2},n_{3},n_{4}]=[784,100,30,10] was used, with M = 100 , T = 7500 , η = 0.05 formulae-sequence 𝑀 100 formulae-sequence 𝑇 7500 𝜂 0.05 M=100,T=7500,\eta=0.05 , tanh activation function and mean squared error cost function.

$ϵ$	0	0.1	0.3	0.5
Standard	96.9%	97.1%	96.9%	96.2%
Low rank	-	96.8%	96.8%	96.3%
$\log \left(1/\gamma \right)/ϵ$	-	30	10	6
$R_a$	-	24.4	21.9	21.3
$R_{\delta }$	-	0.6	0.1	1.3
$R_W$	-	0.1	0.1	0.1

Table 2. Table 2: Average number of correctly labelled test points (out of a total of 30 30 30 ) for the Iris flower data set, for Gaussian distributed noise as well as noise from the RIPE distribution. For the case of RIPE, the R a , R δ subscript 𝑅 𝑎 subscript 𝑅 𝛿 R_{a},R_{\delta} and R W subscript 𝑅 𝑊 R_{W} terms presented are also averaged over 5 5 5 network evaluations. ( ϵ , γ ) italic-ϵ 𝛾 (\epsilon,\gamma) -feedforward neural network accuracy (percentage of correctly labelled test points) were trained on the Iris data set, for various values of ϵ italic-ϵ \epsilon , and γ = 0.05 𝛾 0.05 \gamma=0.05 . A network with L = 3 𝐿 3 L=3 layers and dimensions [ n 1 , n 2 , n 3 ] = [ 4 , 10 , 3 ] subscript 𝑛 1 subscript 𝑛 2 subscript 𝑛 3 4 10 3 [n_{1},n_{2},n_{3}]=[4,10,3] was used, with M = 10 , T = 1200 , η = 0.07 formulae-sequence 𝑀 10 formulae-sequence 𝑇 1200 𝜂 0.07 M=10,T=1200,\eta=0.07 , tanh activation function and mean squared error cost function.

$ϵ$	0	0.1	0.3	0.5
Gaussian	28.6	29.8	29.2	29.0
RIPE	-	28.6	28.4	28.0
$\log \left(1/\gamma \right)/ϵ$	-	30	10	6
$R_a$	-	6.0	4.4	4.0
$R_{\delta }$	-	0.7	0.6	0.6
$R_W$	-	0.4	0.4	0.4

Table 3. Table 3: Comparison of running times between the quantum, quantum-inspired and standard classical algorithms. Note that R a c ​ l ≥ R a 2 superscript subscript 𝑅 𝑎 𝑐 𝑙 superscript subscript 𝑅 𝑎 2 R_{a}^{cl}\geq R_{a}^{2} , R δ c ​ l ≥ R δ 2 superscript subscript 𝑅 𝛿 𝑐 𝑙 superscript subscript 𝑅 𝛿 2 R_{\delta}^{cl}\geq R_{\delta}^{2} and R e c ​ l ≥ R e 2 superscript subscript 𝑅 𝑒 𝑐 𝑙 superscript subscript 𝑅 𝑒 2 R_{e}^{cl}\geq R_{e}^{2} .

	Training	Evaluation
Quantum	$\tilde{O}\left({\left(TM\right)}^{1.5}N\frac{\log \left(1/\gamma \right)}{ϵ}\left(R_a+R_{\delta }+R_W\right)\right)$	$\tilde{O}\left(N\frac{\log \left(1/\gamma \right)}{ϵ}{R}_e\right)$
Quantum-inspired	$\tilde{O}\left({\left(TM\right)}^{2}N\frac{\log \left(1/\gamma \right)}{{ϵ}^2}\left(R_a^{cl}+R_{\delta }^{cl}\right)\right)$	$\tilde{O}\left(N\frac{\log \left(1/\gamma \right)}{{ϵ}^2}{R}_e^{cl}\right)$
Standard classical	$O\left(TME\right)$	$O\left(E\right)$

Equations

$$z_j^l=⟨W_j^l,a^{l-1}⟩+b_j^l,$$

$$W_{jk}^{t+1,l}$$

$$b_j^{t+1,l}$$

$$\left|s_j^l-⟨W_j^l,a^{l-1}⟩\right|\le \max \left\{ϵ\left|⟨W_j^l,a^{l-1}⟩\right|,ϵ\right\}\text{with probability }\ge 1-\gamma ,$$

$$\left|s-⟨x,y⟩\right|\le \{\begin{array}{ll}ϵ\left|⟨x,y⟩\right|& \text{in time }\tilde{O}\left(\frac{T_U\log (1/\gamma )}{ϵ}\frac{\Vert x\Vert \Vert y\Vert }{|⟨x,y⟩|}\right)\\ ϵ& \text{in time }\tilde{O}\left(\frac{T_U\log (1/\gamma )}{ϵ}\Vert x\Vert \Vert y\Vert \right)\end{array}$$

$$W_{jk}^{t,l}=\sum _{\tau =0}^{t-1}\sum _{m=1}^M\frac{-{\eta }^{\tau ,l}}{M}{\delta }_j^{\tau ,m,l}{a}_k^{\tau ,m,l-1}$$

$$\left|s_j^{t,m,l}-⟨W_j^{t,l},a^{t,m,l-1}⟩\right|\le \max \{ϵ|⟨W_j^{t,l},a^{t,m,l-1}⟩|,ϵ\}\text{mbox}withprobability\ge 1-\gamma$$

$$T_{RIPE}(W_j^{t,l},a^{t,m,l-1})=\tilde{O}\left(\frac{T_U\log (1/\gamma )}{ϵ}\frac{\Vert W_j^{t,l}\Vert \Vert a^{t,m,l-1}\Vert }{\max \left\{1,\left|⟨W_j^{t,l},a^{t,m,l-1}⟩\right|\right\}}\right),$$

$$\tilde{O}\left(\sqrt{TM}N\frac{\log (1/\gamma )}{ϵ}{R}_a^{t,m}\right),$$

$$\left|s_j^{t,m,l}-⟨{\left(W^{t,l+1}\right)}^j,{\delta }^{t,m,l+1}⟩\right|\le \max \{ϵ|⟨{\left(W^{t,l+1}\right)}^j,{\delta }^{t,m,l+1}⟩|,ϵ\}\text{mbox}withprobability\ge 1-\gamma$$

$$\tilde{O}\left(\sqrt{TM}N\frac{\log (1/\gamma )}{ϵ}{R}_{\delta }^{t,m}\right),$$

$$b_j^{t+1,l}$$

$$\left|s-⟨x,y⟩\right|\le \{\begin{array}{ll}ϵ\left|⟨x,y⟩\right|& \text{in time }\tilde{O}\left(\frac{\log (1/\gamma )}{{ϵ}^2}\frac{{\Vert x\Vert }^2{\Vert y\Vert }^2}{|⟨x,y⟩|}\right)\\ ϵ& \text{in time }\tilde{O}\left(\frac{\log (1/\gamma )}{{ϵ}^2}{\Vert x\Vert }^2{\Vert y\Vert }^2\right)\end{array}$$

$$\tilde{O}\left(N\frac{\log (1/\gamma )}{{ϵ}^2}{R}_e^{cl}\right),$$

$$\tilde{O}\left((TM{)}^{2}N\frac{\log (1/\gamma )}{{ϵ}^2}\left(R_a^{cl}+R_{\delta }^{cl}\right)\right)$$

$${\Vert X^{[t,l,j]}\Vert }_F=\sqrt{\sum _{\tau =0}^{t-1}\sum _{\mu =1}^M{\left(\frac{\eta }{M}{\delta }_j^{\tau ,\mu ,l}\Vert a^{\tau ,\mu ,l-1}\Vert \right)}^2}$$

$$R_a^{t,m}$$

$$\left|s^{\prime }-⟨x,y⟩\right|$$

$$\le$$

$$|X^{[t,l,j]}⟩=\frac{1}{{\Vert X^{[t,l,j]}\Vert }_F}\sum _{\tau =0}^{t-1}\sum _{\mu =1}^M-\frac{{\eta }^{\tau ,l}}{M}{\delta }_j^{\tau ,\mu ,l}\Vert a^{\tau ,\mu ,l-1}\Vert |\tau ⟩|\mu ⟩$$

$$U_X|t⟩|l⟩|j⟩|0⟩|0⟩$$

$$U_a|l⟩|t⟩|m⟩|0⟩$$

$$|t⟩|l⟩|j⟩\frac{1}{{\Vert X^{[t,l,j]}\Vert }_F}\sum _{\tau =0}^{t-1}\sum _{\mu =1}^M-\frac{{\eta }^{\tau ,l}}{M}{\delta }_j^{\tau ,\mu ,l}|\tau ⟩|\mu ⟩\sum _k{a}_k^{\tau ,\mu ,l-1}|k⟩.$$

$$|t⟩|l⟩|j⟩\frac{1}{{\Vert X^{[t,l,j]}\Vert }_F\sqrt{Mt}}\sum _x\sum _y\left(\sum _{k=1}^N\sum _{\tau =0}^{t-1}\sum _{\mu =1}^M-\frac{{\eta }^{\tau ,l}}{M}(-1{)}^{\tau \cdot x+\mu \cdot y}{\delta }_j^{\tau ,\mu ,l}{a}_k^{\tau ,\mu ,l-1}|k⟩\right)|x⟩|y⟩$$

$$=|t⟩|l⟩|j⟩\left(\sin \theta |W_j^{t,l}⟩|0⟩|0⟩+\cos \theta |\text{junk}⟩|00^{\perp }⟩\right),$$

$$|W_j^{t,l}⟩=\frac{1}{\Vert W_j^{t,l}\Vert }\sum _{k=1}^{n_l}{W}_{jk}^{t,l}|k⟩=\frac{1}{\Vert W_j^{t,l}\Vert }\sum _{k=1}^N\left(\sum _{\tau =0}^{t-1}\sum _{\mu =1}^M\frac{-{\eta }^{\tau ,l}}{M}{\delta }_j^{\tau ,\mu ,l}{a}_k^{\tau ,\mu ,l-1}\right)|k⟩,$$

$$T_W=O\left(\frac{{\Vert X^{[t,l,j]}\Vert }_F}{\Vert W_j^{t,l}\Vert }\sqrt{TM}\mathrm{polylog}(TMN)\right).$$

$$|{\tilde{X}}^{[t,l,j]}⟩=\frac{1}{{\Vert {\tilde{X}}^{[t,l,j]}\Vert }_F}\sum _{\tau =0}^{t-1}\sum _{\mu =1}^M-\frac{{\eta }^{\tau ,l}}{M}{a}_j^{\tau ,\mu ,l-1}\Vert {\delta }^{\tau ,\mu ,l}\Vert |\tau ⟩|\mu ⟩.$$

$$a^{\tau ,\mu ,l-1}$$

$$a^{t,m,l-1}$$

$$X^{[t,l,j]}$$

$$⟨Z⟩$$

$$=\sum _{k,\tau ,\mu }\frac{a_k^{\tau ,\mu ,l-1}}{\Vert a^{\tau ,\mu ,l-1}\Vert }{\left(X^{[t,l,j]}\right)}_{\tau ,\mu }{a}_k^{t,m,l-1}$$

Full text

I Introduction

Machine learning has been one of the great success stories of computation in recent times. From self-driving cars and speech recognition to cancer detection and product recommendation, the applications of machine learning have had a transformative impact on our lives. At the heart of many machine learning tasks are artificial neural networks: groups of interconnected computational nodes loosely modelled on the neurons in our brains. In a supervised learning scenario, a neural network is first trained to recognize a set of labelled data by learning a hierarchy of features that together capture the defining characteristics of each label. Once trained, the network can then be evaluated on previously unseen data to predict their corresponding labels. While their origins can be traced back to the 1940s, artificial neural networks enjoyed a resurgence of interest in the 1980s when a method for computing gradients known as backpropagation rumelhart1986learning was popularized, which dramatically increased the efficiency with which networks could be trained to recognize data. The following decades have seen researchers develop ever more sophisticated techniques to improve the performance of neural networks, and the best algorithms can now outperform humans on realistic image recognition tasks he2016deep . Regardless of the additional techniques used, backpropagation remains the key algorithmic component for neural network training.

In spite of the progress made over the years, significant resources are still required to train deep networks needed for industrial and commercial grade problems, and long training times on clusters of GPUs are often necessary. Furthermore, a network trained on a particular data set for a particular task may not perform well on another task. For this, an entirely different network may need to be trained from scratch. There is thus considerable benefit in any method for speeding up or otherwise improving the training and evaluation of neural networks.

The developments in machine learning have been accompanied by the rise of another potentially transformative technology: quantum computing. From the original proposal of Feynman feynman1982simulating in the 1980s for using the rules of quantum mechanics to carry out computation, quantum computing has developed into one of the most exciting fields of research today. The last several years have seen increased attention on applications of quantum computing for industry, of which quantum chemistry, material design, optimization and quantum machine learning are leading candidates. A natural and fundamentally important question at the intersection of machine learning and quantum computing is whether quantum algorithms can offer any improvement on the classical algorithms currently used for neural networks.

Several challenges present themselves immediately. Firstly, the process of training and evaluating neural networks is highly sequential. Quantum algorithms can be a natural fit for performing tasks in parallel by harnessing the phenomenon of quantum superposition. However, they are poorly suited to situations where data must be computed and stored at many intermediate steps, a process which destroys quantum coherence. Secondly, neural network performance relies critically on the ability to perform non-linear transformations of the data. Quantum mechanics, however, is linear, and effectively implementing non-linear transformations can be non-trivial cao2017quantum . Finally, if a quantum algorithm is to be used for training a classical neural network, the training data as well as the parameters which define the network must first be encoded in quantum states. Unless there is an efficient method for preparing these quantum states, the state-preparation procedure can be a time-consuming bottleneck preventing a quantum algorithm from running more quickly than classically. This is especially challenging for the parameters corresponding to the network weights which are so numerous as to rule out efficient methods for generating their corresponding quantum states directly.

There is some qualitative cause for optimism. The standard classical neural network algorithms typically employed make heavy use of linear algebra, for which quantum algorithms may have an advantage in certain cases biamonte2017quantum . Additionally, it is desirable for neural networks to be robust to noise and small errors, which can achieved classically by introducing perturbations during the training process srivastava2014dropout . In quantum computing, the randomness inherent in the outcome of measurements can have the effect of such introduced perturbations and, if harnessed correctly, a quantum algorithm may achieve a natural robustness to noise without additional cost. This is important as, in many cases, raw computational speed is not necessarily the main concern in practical machine learning. Other issues such as the size of the data sets, generalization errors, or how robust the algorithms are to noise and perturbations may indeed be of much greater practical concern. Nevertheless, there is an explicit connection between speed and performance. In many cases, the bottleneck is the fact that one can only spend, say, a day for training and not a year. With this restriction, one selects a neural network with a given size so that training can be completed in the allotted time. With faster training, larger neural networks can be trained, and more experiments can be conducted for different choices of network architecture, loss function and hyperparameters. Taken together, these can all lead to better eventual performance.

In this work we consider quantum algorithms for feedforward neural networks, in which the propagation of information proceeds from start to finish, layer by layer, without any loops or cycles. Feedforward neural networks not only constitute a canonical category of neural network of fundamental importance in their own right, they also form key components in practical machine learning architectures such as deep convolutional networks and autoencoders. We propose quantum algorithms for training and evaluating robust versions of classical feedforward neural networks—which we call $(\epsilon,\gamma)$-feedforward neural networks—based on the standard classical algorithms for feedforward and backpropagation, adapted to overcome the challenges and exploit the opportunities mentioned above. To achieve this, we use the fact that quantum computers can compute approximate inner products of vectors efficiently, and define a robust inner product estimation procedure that outputs such an estimate in a quantum data register. By having inner products stored in register, rather than in the phases of quantum states, we are able to directly implement nonlinear transformations on the data, circumventing the problem of non-linearity entirely. We address the state-preparation problem for short length vectors by storing them in a particular data-structure kerenidis2017quantum ; kerenidis2018gradient ; wossnig2018quantum in quantum random access memory (qRAM) giovannetti2008quantum ; prakash2014quantum , so that their elements may be queried in quantum superposition, and their corresponding states generated efficiently. For the network weight matrices which are too big to be stored efficiently in this manner, we reconstruct their corresponding quantum states indirectly by superposing histories of shorter length vectors stored in qRAM.

In a recent important theoretical work, Tang tang2018quantum showed that the data structure required for fast qRAM-based inner product estimation can also be used to classically estimate inner products. Quantum algorithms based on such data structures can thus give rise to quantum-inspired classical ones as well tang2018quantum02 ; gilyen2018quantum , i.e. classical algorithms based on their quantum counterparts can be defined with only a polynomial slow-down in running time. However, in practice the polynomial factors can make a difference, and analysis of a number of such algorithms arrazola2019quantum shows that care is needed when assessing their performance relative to the quantum algorithms from which they were inspired. The situation is the same with the quantum algorithms we propose here. New classical algorithms for training and evaluating $(\epsilon,\gamma)$-feedforward neural networks based on our quantum algorithms can indeed be defined, and in a later section we compare these with both their quantum versions as well as the standard classical algorithms. While interesting from a theoretical perspective, these quantum-inspired algorithms will always have a worse performance than their quantum versions, and it is unclear whether they will ever perform well enough in practice to replace existing methods.

Quantum machine learning schuld2015introduction ; biamonte2017quantum ; dunjko2018machine in general and quantum supervised learning lloyd2013quantum ; rebentrost2014quantum ; wiebe2015quantum ; liu2015fast ; schuld2016prediction ; qSFA in particular are fast growing areas of research. However, while quantum generalizations of neural networks have been proposed for feedforward networks cao2017quantum ; wan2017quantum ; romero2017quantum ; farhi2018classification ; liao2018quantum ; schuld2019quantum , Boltzmann Machines wiebe2016quantum ; verdon2017quantum ; kieferova2017tomography ; wiebe2019generative and Hopfield networks rebentrost2018quantum , the current work is, to our knowledge, the first proposed quantum algorithms for training and evaluating feedforward neural networks that can, in principle, offer an advantage over the canonical classical feedforward and backpropagation algorithms. While it remains to be seen whether our algorithms can be used in the future to obtain a practical advantage over classical methods, we believe, backed with our theoretical analysis and initial simulations, that our results indicate one very promising way that quantum techniques can be applied to neural networks, and that in the future more sophisticated practical techniques may be built on these or other ideas.

II Results

We will discuss $(\epsilon,\gamma)$-feedforward neural networks in more detail later. For now, it is enough to think of these as network in which, with probability at least $1-\gamma$, inner products between vectors can be computed to within error $\epsilon$ of the true value. Our main results are the following.

Quantum Training: [See Theorem 1] There exists a quantum algorithm for training $(\epsilon,\gamma)$-feedforward neural networks in time $\tilde{O}\left((TM)^{1.5}N\frac{\log(1/\gamma)}{\epsilon}R\right)$, where $T$ is the number of update iterations, $M$ the number of input samples in each mini-batch, $N$ is the total number of neurons in the network, and $R$ is a factor that depends on the network and training samples, and which numerical evidence suggests is small for many practical parameter regimes.

Here and in what follows, $\tilde{O}$ hides polylogarithmic factors in $TMN$. Once the neural network is trained it can be used to label new input data. This evaluation is essentially the same as the feedforward algorithm that is used in the training, and we obtain the following result.

Quantum Evaluation: [See Theorem 2] For a neural network whose weights are explicitly stored in a qRAM data structure, there exists a quantum algorithm for evaluating an $(\epsilon,\gamma)$-feedforward neural network in time $\tilde{O}\left(N\frac{\log(1/\gamma)}{\epsilon}R_{e}\right)$, where $N$ is the total number of neurons in the network, and $R_{e}$ is a factor that depends on the network and training samples, and is expected to be small for practical parameter regimes.

In contrast, the classical training algorithm has running time $O(TME)$ while the evaluation algorithm takes time $O(E)$, where $E$ is the total number of edges in the neural network. It is interesting to note that the complexity order reduces from the number of edges (neural connections) classically to the number of vertices (neurons) in the quantum algorithm. This gap can be very large; indeed, an average neuron in the human brain has 7000 synaptic connections to other neurons drachman2005we .

The running time of our quantum algorithms also depends on $\log(1/\gamma)/\epsilon$ and on terms $R$ and $R_{e}$ respectively, which depend on the values of certain vector norms which appear during the network evaluation and will be explicated in later sections. While the combined effect of these terms must be taken into account, we will give arguments and numerical evidence that indicate that, in practice, these terms do not contribute significantly to the running time.

II.1 Feedforward neural networks

Feedforward neural networks are covered by a vast amount of literature, and good introductory references are available goodfellow2016deep ; nielsen2017deep . Here we simply give a brief summary of the concepts required for the current work.

A feedforward neural network consists of $L$ layers, with the $l$-th layer containing $n_{l}$ neurons. A weight matrix $W^{l}\in\mathbb{R}^{n_{l}\times n_{l-1}}$ is associated between layers $l-1$ and $l$, and a bias vector $b^{l}\in\mathbb{R}^{n_{l}}$ is associated to each layer $l$ (Fig.1). The total number of neurons is $N=\sum_{l=1}^{L}n_{l}$ and the total number of edges in the network is $E=\sum_{l=2}^{L}n_{l}\cdot n_{l-1}$. For each level $l$, let $a^{l}\in\mathbb{R}^{n_{l}}$ to be the vector of outputs (activations) of the neurons. Given a non-linear activation function $f$, the network feedforward rule is given by $a^{l}_{j}=f(z^{l}_{j})$, where $z^{l}_{j}=\sum_{k}W^{l}_{jk}a^{l-1}_{l}+b^{l}_{j}$ (Fig 2). It will be convenient to express this as

[TABLE]

where $\langle x,y\rangle$ is the Euclidean inner product between vectors $x$ and $y$, and $W^{l}_{j}$ is the $j$-th row of $W^{l}$, i.e. $\left(W^{l}_{j}\right)_{k}=W^{l}_{jk}$. For simplicity, we assume that the activation function is the same across all layers in the network, although our algorithm works more generally when the activation function is allowed to vary across layers. Common activation functions are the sigmoid $\sigma(z)=\frac{1}{1+e^{-z}}$ (when neurons take values in $[0,1]$), the hyperbolic function $\tanh(z)=2\sigma(2z)-1$ (when neurons take values in $[-1,1]$), and the Rectified Linear Unit (ReLU function) $f(x)=\max(0,x)$. A training set $\mathcal{T}=\{(x^{1},y^{1}),(x^{2},y^{2}),\ldots\}$ for a learning task consists of vectors $x^{i}\in\mathbb{R}^{n_{1}}$ and corresponding labels $y^{i}\in\mathbb{R}^{n_{L}}$. This set is used, via a training algorithm, to adjust the network weights and biases so as to minimize a chosen cost function $C:\mathbb{R}^{n_{L}}\rightarrow\mathbb{R}$, which quantifies the network performance. The goal is to obtain parameters such that when the network evaluates a new input which is not part of the training set, it outputs the correct label with a high degree of accuracy.

Classically, network training consists of the following steps:

Initialization of weights and biases. Various techniques are used, but a common choice is to set the biases to a small constant, and to draw the weight matrices randomly from a normal distribution according to $W^{l}_{jk}\sim\mathcal{N}\left(0,\frac{1}{\sqrt{n_{l-1}}}\right)$ glorot2010understanding . 2. 2.

Feedforward. Select a pair $(x,y)$ from the training set. Assign the neurons in the first layer of the network to have activations $a^{1}_{j}=x_{j}$. Pass through the network layer-by-layer, at each layer $l$ computing and storing the vectors $z^{l}=\left\langle W^{l}_{j},a^{l-1}\right\rangle+b^{l}_{j}$ and $a^{l}=f\left(z_{j}^{l}\right)$. The running time of this procedure is $O\left(\sum_{l=2}^{L}n_{l}\cdot n_{l-1}\right)=O\left(E\right)$: at each level $l$ one must evaluate $n_{l}$ activations, and each activation involves calculating an inner product of dimension $n_{l-1}$ which takes times $n_{l-1}$. 3. 3.

Backpropagation. Given $z^{L}$ and $a^{L}$ at the end of the feedforward process, define vector $\delta^{L}$ with components $\delta^{L}_{j}=f^{\prime}\left(z^{L}_{j}\right)\frac{\partial C}{\partial a^{L}_{j}}$, where $C$ is the chosen cost function. Then, proceeding backwards through the network, compute and store the vectors $\delta^{l}_{j}=f^{\prime}\left(z^{l}_{j}\right)\left\langle\left(W^{l+1}\right)^{T}_{j},\delta^{l+1}\right\rangle$, where $\left(W^{l+1}\right)^{T}$ is the matrix transpose of $W^{l+1}$. The running time of the backpropagation algorithm is again $O(E)$, by the same reasoning as for the feedforward algorithm. 4. 4.

Update weights and biases. Repeat the feedforward and backpropagation steps for a mini-batch of inputs of size $M$, and then update the weights and biases by stochastic gradient descent, according to:

[TABLE]

where the superscipts $t\in[T]$ and $m\in[M]$ denote the iteration number and mini-batch element respectively, and $\eta^{t,l}>0$ are update step sizes. 5. 5.

Iterate. Repeat the previous weight and bias update step $T$ times, each time with a different mini-batch.

Note that the method for choosing each mini-batch of inputs is left to the user. A common choice in practice, and the one that used in our numerical simulations, is to divide the training procedure into a number of periods known as epochs. At the beginning of each epoch, the training set is randomly shuffled and successive mini-batches of size $M$ are chosen without replacement. The epoch ends when all training examples have been selected, after which the next epoch begins. Alternatively, one may consider selecting the mini-batches by randomly sampling-with-replacement from the full training set.

The overall running time of the classical training algorithm is $O(TME)$ since there are $TM$ steps and each step requires performing algorithms to carry out the feedforward and backpropagation procedures. Note that the product $TM$ is equal to the number of epochs multiplied by the size $\left|\mathcal{T}\right|$ of the training set. The fact that the training of general feedforward neural networks depends linearly on the number of edges makes large-size fully-connected feedforward networks computationally expensive to train. Our quantum algorithm reduces the dependence on the network size from the number of edges to the number of vertices.

Once the network training is complete, the classical evaluation algorithm consists of running the feedforward step on a previously unseen input, i.e. on data that lies outside the training set. The success of deep learning tells us that the network training produces weights and biases that perform well in classifying new data.

II.2 Robust feedforward neural networks

A perennial concern in neural network training is the problem of overfitting. Since a large network may have many more free parameters than training data points, it is easy to train a network that recognizes the training data accurately, yet performs poorly on data that lies outside the training set. Various techniques are used in practice to reduce overfitting. For instance, one may consider adding a term to the cost function which penalizes a network with too many large parameters. Alternatively, one can train the data on an ensemble of different networks, and average the performance across the ensemble. While this may be prohibitively costly to perform exactly, the effect can be approximated by randomly deactivating certain neurons in each layer, a technique known as dropout. A similar result can be achieved by adding random multiplicative Gaussian noise to the activation of each neuron during training. These techniques are both seen to be effective at preventing overfitting srivastava2014dropout , and can also improve network performance by avoiding local minima in parameter space.

Motivated by the benefits of such random network perturbations during training, and anticipating the requirement of quantum processes to generate randomized outcomes, we consider a generalization of the training and evaluation algorithms where the inner products in the feedforward and backpropagation steps may not be evaluated exactly. Instead, with probability $1-\gamma$, they are estimated to within some error tolerance $\epsilon$, either relative or absolute, depending on whether the inner product is large or small. That is, the feedforward step computes values $s^{l}_{j}$ satisfying

[TABLE]

and similarly for the inner product calculation between $\left(W^{l+1}\right)_{j}^{T}$ and $\delta^{l+1}$ in the backpropgation step. Note that we do not specify how the $s^{l}_{j}$ are generated above. The way this is realized will be left to specific implementations of the algorithm. For instance, a simple classical implementation is to first compute the inner product $\left\langle W_{j}^{l},a^{l-1}\right\rangle$ and then add independent Gaussian noise bounded by the maximum of $\epsilon\left|\left\langle W_{j}^{l},a^{l-1}\right\rangle\right|$ and $\epsilon$. In the quantum case, the $s^{l}_{j}$ will be generated by a quantum inner product procedure which is not perfect, but rather outputs an estimate of the true inner product satisfying the $(\epsilon,\gamma)$ conditions required. The reason we allow for either a relative or absolute error is to ensure that the quantum procedure can be carried out efficiently regardless of the magnitude of the inner products involved.

We refer to such networks as $(\epsilon,\gamma)$-feedforward neural networks, of which the standard classical feedforward neural network is a special case. Our simulation results show that, for reasonable tolerance parameters, the generalization to $(\epsilon,\gamma)$ estimates of inner products does not hurt network performance.

For small enough $\epsilon$, the running time of classically computing the feedforward and backpropagation steps remains $O(E)$, since in general one needs to look at a large fraction of the coordinates of two vectors to obtain an $\epsilon$-error approximation to their inner product. The classical training time for $(\epsilon,\gamma)$-feedforward neural networks is thus $O(TME)$, as for the original case where inner products are evaluated exactly.

II.3 Quantum Training

While it is clearly desirable to improve on the $O(TME)$ classical result using a quantum algorithm, there are several obstacles to achieving this. As mentioned in the introduction, feedforward neural network training and evaluation is a highly sequential procedure, where at each point one needs to know the results of previously computed steps. Quantum algorithms, by contrast, are typically well suited to performing tasks in parallel, but not for performing tasks which require sequential measurements to be performed, a process which destroys quantum coherence. In addition, a critical step classically is the application of a non-linear activation function to each neuron. Given that quantum mechanics is inherently linear, applying non-linearity to quantum states is non-trivial. Finally, the size of each weight matrix is $n_{l}\times n_{l-1}$, so even explicitly writing down these matrices for every step of the algorithm takes time $O(E)$.

We address these challenges by using a hybrid quantum-classical procedure which follows the classical training algorithm closely. In our algorithm all sequential steps are taken classically. At each step, quantum operations are only invoked for estimating the inner products of vectors, and for reading and writing data to and from qRAM.

Given a vector $x\in\mathbb{R}^{n}$, define the corresponding normalized quantum state $\left|x\right\rangle=\frac{1}{\left\|x\right\|}\sum_{j=1}^{n}x_{j}\left|j\right\rangle$, where $\left\|\cdot\right\|$ is the $\ell_{2}$ norm. The inner product of two quantum states therefore satisfies $\langle x|y\rangle:=\frac{\left\langle x,y\right\rangle}{\left\|x\right\|\left\|y\right\|}$. Two key ingredients are the following:

Robust Inner Product Estimation (RIPE). If quantum states $\left|x\right\rangle$ and $\left|y\right\rangle$ can each be created in time $T_{U}$, and if estimates of the norms $\left\|x\right\|$ and $\left\|y\right\|$ are known to within $\epsilon/3$ multiplicative error, then a generalization of the inner product estimation subroutine of kerenidis2018qmeans allows one to perform the mapping $\left|x\right\rangle\left|y\right\rangle\left|0\right\rangle\rightarrow\left|x\right\rangle\left|y\right\rangle\left|s\right\rangle$ where, with probability at least $1-\gamma$:

[TABLE]

The inner product estimates $s$ above are computed in a quantum register and not in the phase of a quantum state, enabling non-linear activation functions to be applied to it. When the data required is stored in qRAM (see below), this inner product calculation is efficient and provides roughly a factor $O(n_{l})$ saving in time per layer of the network.

Quantum Random Access Memory. For the RIPE algorithm to efficiently compute an approximation of the inner product $\left\langle x,y\right\rangle$, we need an efficient way to prepare the states $\left|x\right\rangle$ and $\left|y\right\rangle$. A qRAM is a device that allows for classical data be queried in superposition. That is, if the classical vector $x\in\mathbb{R}^{N}$ is stored in qRAM, then a query to the qRAM implements the unitary $\sum_{j}\alpha_{j}\left|j\right\rangle\left|0\right\rangle\rightarrow\sum_{j}\alpha_{j}\left|j\right\rangle\left|x_{j}\right\rangle$. Importantly, if the elements $x_{j}$ of $x$ arrive as a stream of entries $(j,x_{j})$ in some arbitrary order, then $x$ can be stored in a particular data structure kerenidis2017quantum —which we will refer to as an $\ell_{2}$-Binary Search Tree ($\ell_{2}$-BST)—in a time linear in $N$ (up to logarithmic factors) and, once stored, $\left|x\right\rangle=\frac{1}{\left\|x\right\|}\sum_{j}x_{j}\left|j\right\rangle$ can be created in time polylogarithmic in $N$ (Fig.3).

Conceptually, we would like to follow the classical training algorithm, modified so that the network biases $b^{t,l}$, weights $W^{t,l}$, pre-activations $z^{t,m,l}$, activations $a^{t,m,l}$ and backpropagation vectors $\delta^{t,m,l}$ are stored in a qRAM $\ell_{2}$-BST at every step. Their corresponding quantum states can then be efficiently created, and the RIPE algorithm used to estimate the inner products required to update the network. However, there is a problem. The weight matrices $W^{t,1},W^{t,2},\ldots,W^{t,L}$ have a combined total of $E$ entries. Thus it will take time $\tilde{O}(E)$ just to write all the matrix elements into qRAM. This is true even for the initial weight matrices if they are generated by choosing independent random variables for each element, as is common classically. Each time the weights are updated, it will take an additional $\tilde{O}(E)$ time to write the new values to qRAM. Furthermore, the RIPE algorithm requires estimates of the norms of the vectors involved. This is not a problem for the $a^{t,m,l}$ and $\delta^{t,m,l}$ vectors, as their storage in the $\ell_{2}$-BST allows their norms to be accessed. However, if the weight matrices are not stored in this way, then the norms of their rows and columns are also not available to use explicitly.

We circumvent these issues in two ways.

Low rank initialization. We set the initial weights $W^{1,l}$ to be low rank by selecting a small number $r$ of pairs of random vectors $a^{0,m,l}$ and $\delta^{0,m,l}$ ($m=1,\ldots,r$) and taking the sums of their outer-products. If we define $\eta^{0,1}=-1$ then equation (2) for the weights can be expressed as

[TABLE]

where, for a fixed $t$ and $l$, only $r$ of the $M$ possible $\delta^{0,m,l}$ and $\delta^{0,m,l}$ are non-zero. We shall see that this low-rank initialization does not affect the network performance in practice, and in section III we provide a classical justification for this initialization as well as numerical evidence supporting the use of $r=O(\log n_{l})$ random pairs.

Implicit storage of weight matrices. For each $t\in[T],l\in\{2,\ldots L\},j\in[n_{L}]$, define the matrix $X^{[t,l,j]}\in\mathbb{R}^{t\times M}$ with matrix elements $\left(X^{[t,l,j]}\right)_{\tau m}=\frac{-\eta^{\tau,l}}{M}\delta_{j}^{\tau,m,l}\left\|a^{\tau,m,l-1}\right\|$, with $\tau\in\{0,\ldots t-1\}$, $m\in[M]$ . We will store all the matrix elements of each $X^{[t,l,j]}$ in qRAM (Fig.4), which allows for efficient computation of the states $\left|W^{t,l}_{j}\right\rangle$ on the fly, as opposed to explictly storing all the values $W^{t,l}_{j,k}$ which is prohibitively expensive. More specifically, in the Methods section we show that by querying the rows of each matrix in superposition over iterations $\tau<t$, it is possible to generate the weight states $\left|W^{t,l}_{j}\right\rangle$ in time $T_{W}=\tilde{O}\left(\frac{\left\|X^{[t,l,j]}\right\|_{F}}{\left\|W^{t,l}_{j}\right\|}\sqrt{TM}\right)$, and estimate $\left\|W^{t,l}_{j}\right\|$ to multiplicative error $\xi$ in time $O(T_{W}/\xi)$ respectively, where $\left\|X^{[t,l,j]}_{F}\right\|$ is the Frobenius norm $X^{[t,l,j]}$ . Similar results apply to generating the states $\left|\left(W^{l}\right)^{T}_{j}\right\rangle$ corresponding to the columns of the weight matrices, and estimating their norms.

With these ideas, one can define a quantum $(\epsilon,\gamma)$-feedforward algorithm which adapts the classical one to make use of RIPE, qRAM access and the implicit storage of the weight matrices.

Subroutine 1**.**

(Quantum $(\epsilon,\gamma)$-Feedforward)

Inputs: indices $t\in[T],m\in[M]$; input pair $(x^{t,m},y^{t,m})$ in qRAM; unitaries $U_{W_{j}^{t,l},a^{t,m,l-1}}$ for creating $\left|W^{t,l}_{j}\right\rangle$ and $\left|a^{t,m,l-1}\right\rangle$ in time $T_{U}$, estimates of their norms $\overline{\left\|W^{t,l}_{j}\right\|}$ and $\overline{{\left\|a^{t,m,l-1}\right\|}}$ to relative error at most $\xi=\epsilon/3$, and vectors $b^{t,l}$ in qRAM for $l\in[L]$; activation function $f$; accuracy parameters $\epsilon,\gamma>0$.

For $j=1$ to $n_{1}$ do: 2. 2.

$a^{t,m,1}_{j}=x^{t,m}_{j}$ 3. 3.

For $l=2$ to $L$ do: 4. 4.

For $j=1$ to $n_{L}$ do: 5. 5.

Use the RIPE algorithm with unitary $U_{W_{j}^{t,l},a^{t,m,l-1}}$ to compute $s_{j}^{t,m,l}$, such that

[TABLE] 6. 6.

Compute $z_{j}^{t,m,l}=s_{j}^{t,m,l}+b_{j}^{t,l}$ and store $z_{j}^{t,m,l}$ in qRAM 7. 7.

Compute $a_{j}^{t,m,l}=f(z_{j}^{t,m,l})$ and store $a_{j}^{t,m,l}$ in qRAM

The cost of performing this feedforward procedure is $O\left(\sum_{l=2}^{L}\sum_{j=1}^{n_{l}}T_{RIPE}(W^{t,l}_{j},a^{t,m,l-1})\right)=O\left(N\overline{T}_{RIPE}\right)$, where $T_{RIPE}(W^{t,l}_{j},a^{t,m,l-1})$ is the time required to perform Robust Inner Product Estimation between vectors $W^{t,l}_{j}$ and $a^{t,m,l-1}$, and $\overline{T}_{RIPE}$ denotes the average time (over all layers and neurons) to perform this inner product estimation. Using the running time of RIPE that achieves the larger of the absolute or relative estimation errors, one obtains

[TABLE]

where $T_{U}$ is the time required to prepare state $\left|W^{t,l}_{j}\right\rangle$ and $\left|a^{t,m,l-1}\right\rangle$. $\left|a^{t,m,l-1}\right\rangle$ can be created in time polylogarithmic in $N$ since the required data is, by assumption, stored in qRAM. By the implicit storage of weight matrices, $\left|W^{t,l}_{j}\right\rangle$ can be created in time $\tilde{O}\left(\frac{\left\|X^{[t,l,j]}\right\|_{F}}{\left\|W^{t,l}_{j}\right\|}\sqrt{TM}\right)$. The overall running time of the quantum feedforward algorithm is therefore

[TABLE]

where $R^{t,m}_{a}=\frac{1}{N-n_{1}}\sum_{l=2}^{L}\sum_{j=1}^{n_{l}}\frac{\left\|X^{[t,l,j]}\right\|_{F}\left\|a^{t,m,l-1}\right\|}{\max\left\{1,\left|\left\langle W^{t,l}_{j},a^{t,m,l-1}\right\rangle\right|\right\}}$.

The factor $R^{t,m}_{a}$ does not appear in the classical algorithms, and while it is a priori not clear what impact this will have on the running time, we give evidence in the discussion section that this does not impact the running time significantly in practice. The upside is we save a factor of $O(N)$ compared with the classical case, since our algorithm depends linearly on $N$ and not on the number of edges $E$. Last, there is an overhead of $\sqrt{TM}$ which is a consequence of only saving the weight matrices implicitly. For large neural networks one expects that $N>>\sqrt{TM}$.

We can similarly define a quantum backpropagation algorithm:

Subroutine 2**.**

(Quantum $(\epsilon,\gamma)$-Backpropagation)

Inputs: indices $t\in[T],m\in[M]$; input pair $(x^{t,m},y^{t,m})$ in qRAM; vectors $a^{t,m,l},z^{t,m,l},l\in[L]$ in qRAM; unitaries $U_{(W^{t,l+1})^{j},\delta^{t,m,l+1}}$ for creating $\left|(W^{t,l+1})^{j}\right\rangle$ and $\left|\delta^{t,m,l+1}\right\rangle$ in time $T_{U}$, estimates of their norms $\overline{\left\|(W^{t,l+1})^{j}\right\|}$ and $\overline{{\left\|\delta^{t,m,l+1}\right\|}}$ to relative error at most $\xi=\epsilon/3$, and vectors $b^{t,l}$ in qRAM for $l\in[L]$; derivative activation function $f^{\prime}$; parameters $\eta^{t,l}$ for $l\in[L]$; accuracy parameters $\epsilon,\gamma>0$.

For $j=1$ to $n_{L}$ do: 2. 2.

$\delta^{t,m,L}_{j}=f^{\prime}(z^{t,m,L}_{j})\frac{\partial C}{\partial a^{L}_{j}}$. 3. 3.

For $l=L-1$ to $1$ do: 4. 4.

For $j=1$ to $n_{l}$ do: 5. 5.

Use the RIPE algorithm with unitary $U_{(W^{t,l})^{j},\delta^{t,m,l+1}}$ to compute $s_{j}^{t,m,l}$, such that

[TABLE] 6. 6.

Compute $\delta_{j}^{t,m,l}=f^{\prime}(z^{t,m,l}_{j})s_{j}^{t,m,l}$ and store $\delta_{j}^{t,m,l}$, $-\frac{\eta^{t,l}}{M}\delta_{j}^{t,m,l}\left\|a^{t,m,l-1}\right\|$, and $-\frac{\eta^{\tau,l}}{M}\left\|\delta^{t,m,l}\right\|a_{j}^{t,m,l-1}$ in qRAM.

Similar to the feedforward case, the running time of the quantum backpropagation algorithm is

[TABLE]

where $R_{\delta}^{t,m}=\frac{1}{N-n_{l}}\sum_{l=1}^{L-1}\sum_{j=1}^{n_{l}}\frac{\left\|\tilde{X}^{[t,l+1,j]}\right\|_{F}\left\|\delta^{t,m,l+1}\right\|}{\max\left\{1,\left|\left\langle\left(W^{t,l+1}\right)^{T}_{j},\delta^{t,m,l+1}\right\rangle\right|\right\}}$.

The quantum training algorithm consists of running the feedforward and the backpropagation algorithms for all inputs in a mini-batch of size $M$, and iterating this procedure $T$ times. After each mini batch has been processed, we explicitly update the biases $b$ in qRAM, but we do not explicitly update the weights $W$, since this would take time $\tilde{O}(E)$. Instead, we compute an estimate of the norm of the rows and columns of the weight matrices and keep a history of the $a$ and $\delta$ vectors in memory so that we can create the quantum states corresponding to the weights on the fly.

Subroutine 3**.**

(Quantum $(\epsilon,\gamma)$-Training)

Inputs: input pairs $(x^{t,m},y^{t,m})$ for all $t\in T,m\in[M]$, parameters $\eta^{t,l}$, for $t\in[T],l\in[L]$, and $\epsilon,\gamma>0$.

Initialise the weights and biases $W^{1,l},b^{1,l}$ for $l\in[L]$ with a low-rank initialization. 2. 2.

For $t=1$ to $T$ do: 3. 3.

For $m=1$ to $M$ do: 4. 4.

Run the quantum $(\epsilon,\gamma)$-feedfoward algorithm. 5. 5.

Run the quantum $(\epsilon,\gamma)$-backpropagation algorithm. 6. 6.

Compute the biases and update the qRAM with

[TABLE] 7. 7.

Compute the estimates of the norms $\overline{\left\|W^{t+1,l}_{j}\right\|}$ and $\overline{\left\|(W^{t+1,l})^{j}\right\|}$ with relative error $\xi=\epsilon/3$.

Theorem 1**.**

The running time of the quantum training algorithm is $\tilde{O}\left((TM)^{1.5}N\frac{\log(1/\gamma)}{\epsilon}R\right)$, where $R=R_{a}+R_{\delta}+R_{W}$, $R_{a}=\frac{1}{TM}\sum_{t,m}R^{t,m}_{a}$, $R_{\delta}=\frac{1}{TM}\sum_{t,m}R^{t,m}_{\delta}$, and $R_{W}=\frac{1}{T}\sum_{t}(R^{t}_{W_{r}}+R^{t}_{W_{c}})$, with $R^{t}_{W_{r}}=\frac{1}{M}\frac{1}{N-n_{1}}\sum_{l=2}^{L}\sum_{j=1}^{n_{l}}\frac{\left\|X^{[t,l,j]}\right\|_{F}}{\left\|W^{t,l}_{j}\right\|}$ and $R^{t}_{W_{c}}=\frac{1}{M}\frac{1}{N-n_{1}}\sum_{l=2}^{L}\sum_{j=1}^{n_{l}}\frac{\left\|\tilde{X}^{[t,l,j]}\right\|_{F}}{\left\|\left(W^{t,l}\right)^{T}_{j}\right\|}$.

Proof.

The terms $R_{a}$ and $R_{\delta}$ come from the feedforward and backpropagation subroutines. The $R_{W}$ term comes from the estimation of the norms $\left\|W^{t,l}_{j}\right\|$ and $\left\|\left(W^{t,l}\right)^{T}_{j}\right\|$ which only happens once for each mini-batch. For a given $t,l,j$ the estimation of the norm $\overline{\left\|W^{t+1,l}_{j}\right\|}$ takes time $O(T_{W}/\xi)$, with $T_{W}=\tilde{O}\left(\frac{\left\|X^{[t,l,j]}\right\|_{F}}{\left\|W^{t,l}_{j}\right\|}\sqrt{TM}\right)$, and we take $\xi=\epsilon/3$. Hence, we get the ratio $R^{t}_{W_{r}}=\frac{1}{M}\frac{1}{N-n_{1}}\sum_{l=2}^{L}\sum_{j=1}^{n_{l}}\frac{\left\|X^{[t,l,j]}\right\|_{F}}{\left\|W^{t,l}_{j}\right\|}$ and similarly for the estimation of the columns $R^{t}_{W_{c}}=\frac{1}{M}\frac{1}{N-n_{1}}\sum_{l=2}^{L}\sum_{j=1}^{n_{l}}\frac{\left\|\tilde{X}^{[t,l,j]}\right\|_{F}}{\left\|\left(W^{t,l}\right)^{T}_{j}\right\|}$. Then we can define $R_{W}=\frac{1}{T}\sum_{t}(R^{t}_{W_{r}}+R^{t}_{W_{c}})$ as needed.

∎

The trade-off in avoiding an $O\left(E\right)$ scaling is an extra cost factor of $O(\sqrt{TM})$, which makes the overall running time of our quantum training algorithm essentially $\tilde{O}\left((TM)^{1.5}N\right)$ and, while the exact running time also depends on various other factors, loosely speaking has an advantage over the classical training when $\sqrt{TM}\ll N$.

II.4 Quantum Evaluation

Once a neural network has been trained, it can be evaluated on new data to output a predicted label. While the initial training may only occur once, evaluation of new data labels may occur thousands or millions of times thereafter. There are thus large gains to be had from even small improvements to the efficiency of network evaluation, and we realize such an improvement with our second quantum algorithm.

The quantum procedure for neural network evaluation is essentially the same as the quantum feedforward algorithm. Assume there is a new input pair $(x,y)$ that we want to evaluate, and that we have unitaries $U_{W_{j}^{l},a^{l-1}}$ for creating $\left|W^{l}_{j}\right\rangle$ and $\left|a^{l-1}\right\rangle$ in time $T_{U}$, estimates of their norms $\overline{\left\|W^{l}_{j}\right\|}$ and $\overline{{\left\|a^{l-1}\right\|}}$ to relative error at most $\xi=\epsilon/3$, and vectors $b^{l}$ in qRAM for $l\in[L]$. The running time of the quantum feedforward subroutine implies the following.

Theorem 2**.**

There exists a quantum algorithm for evaluating an $(\epsilon,\gamma)$-feedforward neural network in time $\tilde{O}\left(T_{U}N\frac{\log(1/\gamma)}{\epsilon}R_{e}\right)$, where $R_{e}=\frac{1}{N-n_{1}}\sum_{l=2}^{L}\sum_{j=1}^{n_{l}}\frac{\left\|W^{l}_{j}\right\|\left\|a^{l-1}\right\|}{\max\left\{1,\left|\left\langle W^{l}_{j},a^{l-1}\right\rangle\right|\right\}}$ and $T_{U}$ is the time required to prepare any of the states $\left|W^{l}_{j}\right\rangle$ and $\left|a^{l-1}\right\rangle$. For a neural network whose weights are already explictly stored in an $\ell_{2}$-BST, $T_{U}$ is polylogarithmic in $E$, and the total running time is $\tilde{O}\left(N\frac{\log(1/\gamma)}{\epsilon}R_{e}\right)$.

In contrast, the classical evaluation algorithm requires running time $O(E)$. For a network with parameters trained via our quantum training algorithm, we store the network weight matrices $W$ in memory only implicitly, which leads to a time $T_{U}$ scaling as $O(\sqrt{TM})$. In this case the overall running time equates to $\tilde{O}\left(\sqrt{TM}N\frac{\log(1/\gamma)}{\epsilon}R_{e}\right)$.

II.5 Simulation

In this section we show that $(\epsilon,\gamma)$-networks can achieve comparable results to standard neural networks, and give numerical evidence that the hard-to-estimate parameters $R_{a}$, $R_{\delta}$ and $R_{W}$ may not be large for problems of practical interest. We give further arguments to support this in the Discussion section.

We classically simulate the training and evaluation of $(\epsilon,\gamma)$-feedforward neural networks on the MNIST handwritten digits data set consisting of 60,000 training examples and 10,000 test examples. A network with $L=4$ layers and dimensions $[n_{1},n_{2},n_{3},n_{4}]=[784,100,30,10]$ was used, with $M=100,T=7500,\eta=0.05$, tanh activation function and mean squared error cost function $C=\frac{1}{2M}\sum_{m}\left\|y^{t,m}-a^{t,m,L}\right\|^{2}$. For this network size we have $N=924$, $E=81,700$. Our results are summarized in Table 1. ‘Standard’ weight initialization refers to drawing the initial weight matrix elements from appropriately normalized Gaussian distributions $W^{l}_{jk}\sim\mathcal{N}\left(0,\frac{1}{\sqrt{n_{l-1}}}\right)$ glorot2010understanding , and ‘Low rank’ refer to weight initialization as described in Section IV.2, with rank $r=6$. In all cases, Gaussian noise drawn from $N(\epsilon/2,0)$ was added to each inner product evaluated throughout the network. Note that no cost function regularization was used to improve the network performance, and as $\epsilon$ was varied the other network hyper-parameters were not re-optimized. We find that for $\gamma=0.05$ and various values of $\epsilon$, the network achieves high accuracy whilst incurring only modest contributions to the running time from the quantum-related terms $R_{a},R_{\delta}$ and $R_{W}$. For the $\epsilon=0.3$ case we calculate the values of $R_{a}$, $R_{\delta}$ and the two components of $R_{W}$, as a function of the number of gradient update steps $t$. These results appear in Fig. 5, where we see that these values quickly stabilize to small constants during the training procedure.

While the MNIST data set is evaluated using Gaussian distributed noise, we also numerically evaluate the performance of $(\epsilon,\gamma)$-feedforward networks on the Iris flower data set using both Gaussian noise and, for comparison, noise corresponding to the quantum RIPE subroutine where the norms of all vectors are known exactly (see Supplemental Material). This data set—consisting of 120 training examples and 30 test examples, with each data point corresponding to a length $4$ real vector and having one of three labels—is small enough that classically sampling from the RIPE distribution (which corresponds to simulating the output of a quantum circuit) during the network training can be performed efficiently on a standard desktop computer. A network with $L=3$ layers, dimension $[n_{1},n_{2},n_{3}]=[4,10,3]$, tanh activation function and mean squared error cost function was chosen, with $M=10$, $T=1200$ and $\eta=0.07$. For each $(\epsilon,\gamma)$ pair we repeated the network training and evaluation 5 times and show the average number of correctly labelled test points in Table 2. For $\gamma=0.05$ and the same choices of $\epsilon$ as in the MNIST case, we find that the choice of noise distribution, Gaussian or RIPE, does not lead to significant difference in the network performance.

II.6 Quantum-inspired classical algorithms

Our quantum training and evaluation algorithms rely on the $\ell_{2}$-BST data structure in order to efficiently estimate inner products via the RIPE procedure. By a result of Tang, such a data structure in fact allows inner products to be efficiently estimated classically as well:

Classical $\ell_{2}$-BST-based Inner Product Estimation tang2018quantum . If $x,y\in\mathbb{R}^{n}$ are stored in $\ell_{2}$-BSTs then, with probability at least $1-\gamma$, a value $s$ can be computed satisfying

[TABLE]

We can use this concept to derive quantum-inspired classical analogues of our algorithms. If the network weights are explicitly stored in $\ell_{2}$-BSTs, then the quantum RIPE procedure can be directly replaced by this classical inner product estimation routine to give a classical algorithm for network evaluation which has running time

[TABLE]

where $R^{cl}_{e}=\frac{1}{N-n_{1}}\sum_{l=2}^{L}\sum_{j=1}^{n_{l}}\frac{\left\|W^{l}_{j}\right\|^{2}\left\|a^{l-1}\right\|^{2}}{\max\left\{1,\left|\left\langle W^{l}_{j},a^{l-1}\right\rangle\right|^{2}\right\}}\geq R_{e}^{2}$.

Dequantizing the quantum training algorithm is slightly more complicated as the network weights are only stored implicitly. Nonetheless, the required inner products can still be estimated efficiently, and in the Methods section we show a quantum-inspired algorithm can be given which runs in time

[TABLE]

where $R_{a}^{cl}\geq R_{a}^{2}$ and $R_{\delta}^{cl}\geq R_{\delta}^{2}$.

The quantum-inspired algorithm is thus slower than the quantum one by a factor of $\frac{\sqrt{TM}}{\epsilon}\frac{R_{a}^{2}+R_{\delta}^{2}}{\left(R_{a}+R_{\delta}+R_{W}\right)}$. We analyze their relative performance further in the Discussion section.

III Discussion

We have presented quantum algorithms for training and evaluating $(\epsilon,\gamma)$-feedforward neural networks which take time $\tilde{O}\left((TM)^{1.5}N\frac{\log(1/\gamma)}{\epsilon}R\right)$ and $\tilde{O}(NT_{U}\frac{\log(1/\gamma)}{\epsilon}R_{e})$ respectively. These are the first algorithms based on the canonical classical feedforward and backpropagation procedures with running times better than the number of inter-neuron connections, opening the way to utilizing larger sized neural networks. Furthermore, our algorithms can be used as the basis for quantum-inspired classical algorithms which have the same dependence on the network dimensions, but a quadratic penalty in other parameters. Let us now make a number of remarks on our design choices and the performance of our algorithms.

Performance of $(\epsilon,\gamma)$-feedforward neural networks. While the notion of introducing noise in the inner product calculations of an $(\epsilon,\gamma)$-feedforward neural network is, as previously mentioned, similar in spirit to known classical techniques for improving neural network performance, there are certain differences. Classically, dropout and multiplicative Gaussian noise srivastava2014dropout involve introducing perturbations post-activation: $a^{l}_{j}\rightarrow a^{l}_{j}r_{j}$, where $r_{j}\sim N(1,1)$ (multiplicative Guassian noise) or $r_{j}\sim Bernoulli(p)$ (dropout), whereas in our case the noise due to inner product estimation occurs pre-activation. Furthermore, classical methods are typically employed during the training phase but, once the network parameters have been trained, new points are evaluated without the introduction of noise. This is in contrast to our quantum algorithm where the evaluation of new points inherently also involves errors in inner product estimation. However, numerical simulation (see Tables 1 and 2) of the noise model present in our quantum algorithm shows that $(\epsilon,\gamma)$-networks may tolerate modest values of noise, and values of $\epsilon$ and $\gamma$ can be chosen for which the network performance does not suffer significantly, whilst at the same time do not contribute greatly to the factor of $\log(1/\gamma)/\epsilon$ that appears in the quantum running time.

Low rank initialization. One assumption we make in our quantum algorithm is a low-rank initialization of the network weight matrices, compared with freedom to choose full-rank weights classically. This assumption is made in order to avoid a time of $\tilde{O}(E)$ to input the initial weight values into qRAM. Low rank approximations to network weights have found applications classically in both speeding up testing of trained networks sainath2013low ; denton2014exploiting ; yu2017compressing as well as in network training tai2015convolutional , in some cases delivering significant speedups without sacrificing much accuracy. We find numerically that the low-rank initialization we require for our quantum algorithm works as well as full rank for a range of $\epsilon$ values (see Table 1).

Quantum training running time. Compared with the classical running time of $O(TME)$, our quantum algorithm scales with the number of neurons in the network as opposed to the number of edges. However, this comes at the cost of a square root penalty in the number of iterations and mini-batch size, and it is an open question to see how to remove this term. The quantum algorithm also has additional factors of $\log(1/\gamma)/\epsilon$ and $R=R_{a}+R_{\delta}+R_{W}$ which do not feature in the classical algorithm. While the $\epsilon$ and $\gamma$ can be viewed as hyperparameters, which can be freely chosen, the impact of the $R$ terms warrant further discussion.

The ratio $\left\|X^{[t,l,j]}\right\|/\left\|W^{t,l}_{j}\right\|$ appears in both $R_{a}$ and $R_{W}$ (in the $R^{t}_{W_{r}}$ contribution) and similarly the ratio $\left\|\tilde{X}^{[t,l,j]}\right\|/\left\|\left(W^{t,l}\right)^{T}_{j}\right\|$ appears in $R_{\delta}$ and the $R^{t}_{W_{r}}$ contribution to $R_{W}$. While exact values may be difficult to predict, one can expect the following large $t$ behaviour: Classically, initial weight matrices are typically chosen so that the entries are drawn from a normal distribution with standard deviation $1/\sqrt{\text{row length}}$, which would give $\left\|W^{1,l}_{j}\right\|\approx 1$, and a similar scenario can hold with low rank initialization. As the weights are updated according to equation (2), and since changes in individual matrix elements may be positive or negative, for a constant step size $\eta$ one expects $\left\|W^{t,l}_{j}\right\|$ to roughly grow proportionally to $\sqrt{t\eta}$. The norm $\left\|X^{[t,l,j]}\right\|$ has value

[TABLE]

and, as $t$ increases and the network becomes close to well trained, we expect $\delta^{\tau,\mu,l}_{j}\rightarrow 0$, whereas for activations bounded in the range $[-1,1]$, as is the case for the tanh function, $\left\|a^{t,l-1}\right\|\leq\sqrt{n_{l-1}}$. We thus expect $\left\|X^{[t,l,j]}\right\|$ to saturate for large $t$ and not grow in an unbounded fashion. Fig. 5(a) showing the time averaged values of $R^{t}_{W_{r}}$ and $R^{t}_{W_{c}}$ for an $(\epsilon,\gamma)$-network trained on the MNIST data set is consistent with these ratios saturating to very small values over time, in fact values less than $0.014$. A similar result can be seen in Fig. 6(a) for the IRIS data set.

The term $R_{a}=\frac{1}{TM}\sum_{t,m}R_{a}^{t,m}$ is an average over iterations and mini-batch elements of terms $R_{a}^{t,m}$, which themselves are averages over neurons in the network of ratios and products of matrix and vector norms:

[TABLE]

As discussed, we expect $\left\|X^{[t,l,j]}\right\|$ to saturate for large $t$, and for activations in $[-1,1]$ we have $\left\|a^{t,m,l-1}\right\|\leq\sqrt{n_{l-1}}$. However, as the network becomes well trained, one expects the inner products $\left\langle W^{t,l}_{j},a^{t,m,l-1}\right\rangle$ to become large in magnitude so that the neurons have post-activation values close to $\pm 1$. It is thus reasonable to expect the $R_{a}$ to saturate or even decline for large $t$. This is consistent with our results in Fig. 5(b) and 6(b). One expects similar large $t$ behaviour for $R_{\delta}$, except intuitively $R_{\delta}$ should be much smaller than $R_{a}$ since $\left\|\delta^{t,m,l}\right\|$ should become very small as the network becomes well trained. Figs. 5(c) and 6(c) display this expected behaviour. While these simulation results are already promising, we expect the quantum advantage in the running time to become more prominent for larger size neural networks.

Quantum evaluation running time. The classical running time of the evaluation algorithm is $O(E)$. In contrast, our quantum evaluation algorithm runs in time $\tilde{O}\left(N\frac{\log(1/\gamma)}{\epsilon}R_{e}\right)$ if the entries of the weight matrices are explicitly stored in qRAM. By the same arguments given above for the training algorithm, we expect the penalty term $R_{e}$ to be small for well-trained networks, and in this case we expect the quantum training algorithm to be able to provide a significant speed-up over its classical counterpart.

Quantum vs. quantum-inspired running times. While both the quantum and quantum-inspired classical algorithms presented here have a running time linear in $N$, the quantum-inspired algorithms come with a quadratic overhead in other parameters. A comparison of their running times is given in Table 3. The standard classical training time of $TME$ is outperformed by our quantum training algorithm when $\sqrt{TM}\frac{N}{E}\frac{\log(1/\gamma)}{\epsilon}\left(R_{a}+R_{\delta}+R_{W}\right)\ll 1$, whereas the quantum-inspired classical training algorithms can only outperform the standard classical algorithm when $TM\frac{N}{E}\frac{\log(1/\gamma)}{\epsilon^{2}}\left(R_{a}^{cl}+R_{\delta}^{cl}\right)\ll 1$. To understand the significance of the quadratic overheads required by the quantum-inspired algorithms, consider a back-of-the-envelope calculation based on the landmark neural network of Krizhevsky, Sutskever and Hinton krizhevsky2012imagenet . The convolutional neural network they use to recognize the ImageNet LSVRC data set has $TM\approx 10^{8}$, and fully connected final layers corresponding to $N\approx 2\times 10^{4},E\approx 6\times 10^{7}$. Taking $\epsilon=0.1,\gamma=0.05$ and assuming values of $R_{a},R_{\delta}$ and $R_{W}$ the same order of magnitude as we numerically evaluated for the MNIST case gives $\sqrt{TM}\frac{N}{E}\frac{\log(1/\gamma)}{\epsilon}\left(R_{a}+R_{\delta}+R_{W}\right)\approx 10^{3}$ and $TM\frac{N}{E}\frac{\log(1/\gamma)}{\epsilon^{2}}\left(R_{a}^{cl}+R_{\delta}^{cl}\right)\approx 1.4\times 10^{9}$, a factor of $10^{6}$ in favor of the quantum algorithm. While neither the quantum nor the quantum-inspired classical algorithms can compete with the standard classical algorithm in this particular case, plausible changes to the network and data set parameters can be chosen where the quantum training algorithm has an advantage. The quantum-inspired algorithm though would require network parameters many orders of magnitude different to the practical example considered here. In fact, it remains to be seen if there are any real cases where quantum-inspired algorithms can be better in practice than the standard classical ones, and the evidence so far is negative arrazola2019quantum .

Let us add a final remark. In our training algorithm we used classical inputs and showed that the number of iterations required for convergence is similar to the case of classical robust training. One can also consider using superpositions of classical inputs for the training, which could conceivably reduce the number of iterations or size of mini-batch required. We leave this as an interesting open direction for future work.

IV Methods

IV.1 Robust Inner Product Estimation running time.

In kerenidis2018qmeans , the authors give a quantum inner product estimation (IPE) algorithm which allows one to perform the mapping $\left|x\right\rangle\left|y\right\rangle\left|0\right\rangle\rightarrow\left|x\right\rangle\left|y\right\rangle\left|s\right\rangle$, where $s$ satisfies $\left|s-\langle x,y\rangle\right|\leq\epsilon$ with probability at least $1-\gamma$. The running time is $\widetilde{O}\left(\frac{T_{U}\log(1/\gamma)}{\epsilon}\left\|x\right\|\left\|y\right\|\right)$, where $T_{U}$ is the time require to implement unitary operations for creating states $\left|x\right\rangle,\left|y\right\rangle,\left|\left\|x\right\|\right\rangle,\left|\left\|y\right\|\right\rangle$.

The RIPE algorithm is a generalization to the case where one only has access to estimates $\overline{\left\|x\right\|},\overline{\left\|y\right\|}$ of the vector norms, satisfying $\left|\overline{\left\|x\right\|}-\left\|x\right\|\right|\leq\frac{\epsilon}{3}\left\|x\right\|$ and $\left|\overline{\left\|y\right\|}-\left\|y\right\|\right|\leq\frac{\epsilon}{3}\left\|y\right\|$ and where one would like to obtain inner product estimates to either additive or multiplicative error. For normalized vectors $\left|x\right\rangle,\left|y\right\rangle$, the IPE algorithm runs in time $\widetilde{O}\left(\frac{T_{U}\log(1/\gamma)}{\epsilon^{\prime}}\right)$ and outputs $\left|s-\langle x|y\rangle\right|\leq\epsilon^{\prime}$. By taking $\epsilon^{\prime}=\frac{\epsilon}{4}\left|\langle x|y\rangle\right|$ we obtain a relative error algorithm in time $\widetilde{O}\left(\frac{T_{U}\log(1/\gamma)}{\epsilon}\frac{1}{\left|\langle x|y\rangle\right|}\right)$. Outputting the estimator $s^{\prime}=\overline{\left\|x\right\|}\,\overline{\left\|y\right\|}s$ then satisfies

[TABLE]

for small enough $\epsilon$. An absolute error estimate can similarly be obtained by taking replacing $\epsilon$ above with $\epsilon/|\langle x,y\rangle|$. We thus have an algorithm that runs in time $T_{\text{IPE}}(x,y)=\widetilde{O}\left(\frac{T_{U}\log(1/\gamma)}{\epsilon}\frac{\left\|x\right\|\left\|y\right\|}{\max\{1,\left|\langle x,y\rangle\right|\}}\right)$ and achieves an error of $|s-\langle x,y\rangle|\leq\max\left\{\epsilon\left|\left\langle x,y\right\rangle\right|,\epsilon\right\}$.

IV.2 Constructing the weight matrix states and estimating their norms

If, at each iteration $t$, the values $a^{t,m,l}_{j}$ and $-\frac{\eta^{t,l}}{M}\delta_{j}^{t,m,l}\left\|a^{t,m,l-1}\right\|$ are stored in qRAM for all $m\in[M]$, then the states $\left|a^{t,m,l-1}\right\rangle$ and

[TABLE]

can be created coherently in time polylogarithmic in $TMN$. That is, unitary operators $U_{X}$ and $U_{a}$ can be implemented in this time that effect the transformations:

[TABLE]

Application of $U_{X}$ on the first five registers of state $\left|t\right\rangle\left|l\right\rangle\left|j\right\rangle\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle$, followed by application of $U_{a}$ on registers $\{2,4,5,6\}$ produces the state

[TABLE]

Applying the Hadamard transformations $\left|\tau\right\rangle\rightarrow\frac{1}{\sqrt{t}}\sum_{x}(-1)^{\tau\cdot x}\left|x\right\rangle$ and $\left|\mu\right\rangle\rightarrow\frac{1}{\sqrt{M}}\sum_{y}(-1)^{\mu\cdot y}\left|y\right\rangle$ leads to

[TABLE]

where

[TABLE]

$\sin^{2}\theta=\frac{\left\|W^{t,l}_{j}\right\|^{2}}{\left\|X^{[t,l,j]}\right\|_{F}^{2}Mt}$ and $\left|00^{\perp}\right\rangle$ is a state orthogonal to $\left|0\right\rangle\left|0\right\rangle$. By the well-known quantum procedures of amplitude amplification and amplitude estimation brassard2002quantum , given access to a unitary operator $U$ acting on $k$ qubits such that $U\left|0\right\rangle^{\otimes k}=\sin(\theta)\left|x,0\right\rangle+\cos(\theta)\left|G,1\right\rangle$ (where $\left|G\right\rangle$ is arbitrary), $\sin^{2}(\theta)$ can be estimated to additive error $\epsilon\sin^{2}(\theta)$ in time $O\left(\frac{T(U)}{\epsilon\sin(\theta)}\right)$ and $\left|x\right\rangle$ can be generated in expected time $O\left(\frac{T(U)}{\sin(\theta)}\right)$,where $T(U)$ is the time required to implement $U$. Amplitude amplification applied to the unitary preparing the state in (8) allows one to generate $\left|W^{t,l}_{j}\right\rangle$ in time

[TABLE]

Similarly, amplitude estimation can be used to find an $s$ satisfying $\left|s-\frac{\left\|W^{t,l}_{j}\right\|^{2}}{\left\|X^{[t,l,j]}\right\|_{F}^{2}Mt}\right|\leq\xi\frac{\left\|W^{t,l}_{j}\right\|^{2}}{\left\|X^{[t,l,j]}\right\|_{F}^{2}Mt}$ in time $O\left(T_{W}/\xi\right)$. Outputting $\overline{\left\|W_{j}^{t,l}\right\|}=\left\|X^{[t,l,j]}\right\|_{F}\sqrt{Mts}$ then satisfies $\left|\overline{\left\|W^{t,l}_{j}\right\|}-\left\|W^{t,l}_{j}\right\|\right|\leq\xi\left\|W^{t,l}_{j}\right\|$, since $|\sqrt{s}\left\|X^{[t,l,j]}\right\|_{F}\sqrt{Mt}-\left\|W^{t,l}_{j}\right\||\leq\xi\left\|W^{t,l}_{j}\right\|$.

If the values $\delta^{t,m,l}_{j}$ and $\frac{-\eta^{t,l}}{M}a^{t,m,l-1}_{j}\left\|\delta^{t,m,l}\right\|$ are also stored at every iteration, analogous results also hold for creating quantum states $\left|\left(W^{t,l}\right)^{T}_{j}\right\rangle$ corresponding to the columns of $W^{t,l}$, except in this case the key ratio $\left\|W^{t,l}_{j}\right\|/\left\|X^{[t,l,j]}\right\|_{F}$ that appeared in the previous proof is replaced by $\left\|\left(W^{t,l}\right)^{j}\right\|/\left\|\tilde{X}^{[t,l,j]}\right\|_{F}$, where $\left\|\tilde{X}^{[t,l,j]}\right\|_{F}$ is the norm of the quantum state

[TABLE]

IV.3 Quantum-inspired classical algorithm for $(\epsilon,\gamma)$-feedforward network training.

Suppose that for a given $l\in\{2,\ldots,L\}$ and $j\in[n_{l}]$ the following are stored in an $\ell_{2}$-BST:

[TABLE]

It is then possible to sample $(\tau,\mu)$ with probability $P(\tau,\mu)=\frac{\left(X^{[t,l,j]}\right)^{2}_{\tau,\mu}}{\left\|X^{[t,l,j]}\right\|^{2}_{F}}$, and $k$ with conditional probability $P(k\mid\tau,\mu)=\left(\frac{a^{\tau,\mu,l-1}_{k}}{\left\|a^{\tau,\mu,l-1}\right\|}\right)^{2}$ in time $O(\operatorname{polylog}(TMN))$. The random variable $Z(k,\tau,\mu)\stackrel{{\scriptstyle\mathsf{def}}}{{=}}a_{k}^{t,m,l-1}\frac{\left\|a^{\tau,\mu,l-1}\right\|}{a_{k}^{\tau,\mu,l-1}}\frac{\left\|X^{[t,l,j]}\right\|_{F}^{2}}{\left(X^{[t,l,j]}\right)_{\tau,\mu}}$ can be computed in time $O(\operatorname{polylog}(TMN))$, and has expectation

[TABLE]

and variance

[TABLE]

By the Chebyshev and Chernoff-Hoeffding inequalities, taking the median of $O\left(\log(1/\gamma)\right)$ averages, each an average of $O\left(\frac{1}{\epsilon^{\prime 2}}\right)$ independent copies of $Z$, produces an estimate $s$ within $\epsilon^{\prime}\sigma$ of $\langle Z\rangle$. Taking $\epsilon^{\prime}=\frac{\epsilon\max\left\{\left|\langle W^{t,l}_{j},a^{t,m,l-1}\rangle\right|,1\right\}}{\sqrt{MT}\left\|a^{t,m,l-1}\right\|\left\|X^{[t,l,j]}\right\|_{F}}$ then allows an $s^{t,m,l}_{j}$ to be computed satisfying $\left|s^{t,m,l}_{j}-\left\langle W^{t,l}_{j},a^{t,m,l-1}\right\rangle\right|\leq\epsilon\max\left\{\left|\left\langle W^{t,l}_{j},a^{t,m,l-1}\right\rangle\right|,1\right\}$ with probability at least $1-\gamma$.

Replacing the RIPE procedure in Subroutine 1 with this method for computing $s^{t,m,l}_{j}$ gives a quantum-inspired classical $(\epsilon,\gamma)$-feedforward subroutine which runs in time

[TABLE]

where ${R^{cl}_{a}}^{t,m}=\frac{1}{N-n_{1}}\sum_{l=2}^{L}\sum_{j=1}^{n_{l}}\left(\frac{\left\|X^{[t,l,j]}\right\|_{F}\left\|a^{t,m,l-1}\right\|}{\max\left\{1,\left|\langle W^{t,l}_{j},a^{t,ml-1}\rangle\right|\right\}}\right)^{2}$. Using $\left(\sum_{i=1}^{n}x_{i}\right)^{2}\leq n\sum_{i=1}^{n}x_{i}^{2}$ for $x_{i}\geq 0$, it follows that ${R^{cl}_{a}}^{t,m}\geq\left(R_{a}^{t,m}\right)^{2}$.

The RIPE procedure in Subroutine 2 can similarly be replaced to obtain a quantum-inspired classical $(\epsilon,\gamma)$-backpropagation subroutine which runs in time

[TABLE]

where ${R_{\delta}^{cl}}^{t,m}\geq\left(R_{\delta}^{t,m}\right)^{2}$.

Finally, by substituting the quantum $(\epsilon,\gamma)$-feedforward and quantum $(\epsilon,\gamma)$-backpropagation subroutines with their quantum-inspired classical counterparts in Subroutine 3, one obtains a quantum-inspired classical $(\epsilon,\gamma)$-training algorithm which runs in time

[TABLE]

with $R_{a}^{cl}\geq R_{a}^{2}$ and $R_{\delta}^{cl}\geq R_{\delta}^{2}$.

V Supplemental Material

V.1 Classically sampling from the RIPE distribution

The quantum $(\epsilon,\gamma)$-feedforward and quantum $(\epsilon,\gamma)$-backpropagation subroutines (subroutines 1 and 2 respectively) of the main text make use of the Robust Inner Produce Estimation (RIPE) procedure to compute values $s$ that satisfy

[TABLE]

with probability at least $1-\gamma$, for given input vectors $x,y\in\mathbb{R}^{n}$. In this section we give a classical subroutine for generating samples from RIPE distribution, for the special case of RIPE where one has exact knowledge of the the norms $\left\|x\right\|$ and $\left\|y\right\|$ (the generalization to the case where we only have estimates of the norms is straightforward). At a high level, this procedure is based on the following ideas (see [38] for more details):

• •

The inner product estimation procedure of [38] on which RIPE is based implicitly assumes access to a unitary operator for efficiently creating the state $\left|\phi\right\rangle=\frac{\left\|x\right\|\left|0\right\rangle\left|x\right\rangle+\left\|y\right\|\left|1\right\rangle\left|y\right\rangle}{\sqrt{\left\|x\right\|^{2}+\left\|y\right\|^{2}}}$. Applying a Hadamard operator to the first register produces the state

[TABLE]

where $a=\frac{\left\|x\right\|^{2}+\left\|y\right\|^{2}-2\langle x,y\rangle}{2\left(\left\|x\right\|^{2}+\left\|y\right\|^{2}\right)}$. Obtaining an estimate $\bar{a}$ to $a$ satisfying $\left|\bar{a}-a\right|\leq\epsilon_{a}=\frac{\epsilon\max\{1,\left|\langle x,y\rangle\right|\}}{\left\|x\right\|^{2}+\left\|y\right\|^{2}}$ will therefore allow us to compute an estimate $s$ to $\langle x,y\rangle$ satisfying $\left|s-\langle x,y\rangle\right|\leq\epsilon\max\{1,\left|\langle x,y\rangle\right|\}$, by taking $s=\left(\left\|x\right\|^{2}+\left\|y\right\|^{2}\right)\left(1-2\bar{a}\right)/2$.

• •

Performing amplitude estimation [44] on $\left|\Psi\right\rangle$ with $\log M$ ancilla qubits returns a value $\tilde{a}$ satisfying $\left|\tilde{a}-a\right|\leq\frac{\pi}{M}+\left(\frac{\pi}{M}\right)^{2}$ with probability at least $8/\pi^{2}$. This process requires time $\tilde{O}\left(MT\right)$, where $T$ is the time required to implement the unitary required for the creation of state $\left|\phi\right\rangle$. Taking $M=\left\lceil\frac{\pi}{2\epsilon_{a}}\left(1+\sqrt{1+4\epsilon_{a}}\right)\right\rceil$ suffices to ensure that $\left|\tilde{a}-a\right|\leq\epsilon_{a}$.

• •

By the Hoeffding bound, repeating the above procedure $Q=\left\lceil\frac{\log(1/\gamma)}{2(8/\pi^{2}-1/2)^{2}}\right\rceil_{odd}$ times and taking the median of the results gives a value $\bar{a}$ satisfying $\left|\bar{a}-a\right|\leq\epsilon_{a}$ with probability at least $1-\gamma$. The notation $\left\lceil z\right\rceil_{odd}$ denotes the smallest odd integer greater than or equal to $z$.

The amplitude amplification part of the sampling subroutine makes use of a distance function $d:\mathbb{R}^{2}\rightarrow\mathbb{Z}$ given by $d(\omega_{0},\omega_{1})=\min_{z\in\mathbb{Z}}\{z+\omega_{1}-\omega_{0}\}$.

Subroutine 4**.**

(Classically sampling from the RIPE distribution)

Inputs: $x,y\in\mathbb{R}^{n}$, $\epsilon,\gamma>0$, $Q=\left\lceil\frac{\log(1/\gamma)}{2(8/\pi^{2}-1/2)^{2}}\right\rceil_{odd}$.

Compute $a=\frac{\left\|x\right\|^{2}+\left\|y\right\|^{2}-2\langle x,y\rangle}{2\left(\left\|x\right\|^{2}+\left\|y\right\|^{2}\right)}$, $\theta_{a}=\sin^{-1}\left(\sqrt{a}\right)$, $\epsilon_{a}=\frac{\epsilon\max\{1,\left|\langle x,y\rangle\right|\}}{\left\|x\right\|^{2}+\left\|y\right\|^{2}}$, $M=\left\lceil\frac{\pi}{2\epsilon_{a}}\left(1+\sqrt{1+4\epsilon_{a}}\right)\right\rceil$ 2. 2.

For $j=1$ to $M$ do: 3. 3.

$a_{j}=\sin^{2}\left(\pi j/M\right)$ 4. 4.

$p(a_{j})=\left|\frac{\sin\left(Md\left(j/M,\theta_{a}/\pi\right)\right)}{M\sin\left(d\left(j/M,\theta_{a}/\pi\right)\right)}\right|^{2}$ 5. 5.

For $q=1$ to $Q$ do: 6. 6.

Sample $\tilde{a}^{(q)}\sim p$ 7. 7.

Compute $\bar{a}=\mathsf{median}(\tilde{a}^{(1)},\ldots\tilde{a}^{(Q)})$ 8. 8.

Return $s=\left(\left\|x\right\|^{2}+\left\|y\right\|^{2}\right)(1-2\bar{a})/2$

Examples of such samples are shown in Fig 7 for randomly chosen vectors $x$ and $y$ given by

[TABLE]

$\epsilon=0.3$ and various values of $\gamma$. These vectors have norms $\left\|x\right\|=10.34$ and $\left\|y\right\|=14.35$, and inner product $\langle x,y\rangle=95.38$.