Title: QUANTUM ALGORITHMS FOR DEEP CONVOLUTIONAL NEURAL NETWORK

Abstract: Quantum computing is a powerful computational paradigm with applications in several fields, including machine learning. In the last decade, deep learning, and in particular Convolutional Neural Networks (CNN), have become essential for applications in signal processing and image recognition. Quantum deep learning, however, remains a challenging problem, as it is difficult to implement non linearities with quantum unitaries Schuld et al. (2014) . In this paper we propose a quantum algorithm for evaluating and training deep convolutional neural networks with potential speedups over classical CNNs for both the forward and backward passes. The quantum CNN (QCNN) reproduces completely the outputs of the classical CNN and allows for non linearities and pooling operations. The QCNN is in particular interesting for deep networks and could allow new frontiers in the image recognition domain, by allowing for many more convolution kernels, larger kernels, high dimensional inputs and high depth input channels. We also present numerical simulations for the classification of the MNIST dataset to provide practical evidence for the efficiency of the QCNN.
In this paper, we design a quantum convolutional neural network (QCNN) algorithm with a modular architecture that allows for any number of layers, any number and size of kernels, and that can support a large variety of non linearity and pooling methods. Our main technical contributions include a new notion of a quantum convolution product, the development of a quantum sampling technique well suited for information recovery in the context of CNNs and a proposal for a quantum backpropagation algorithm for efficient training of the QCNN.

1. INTRODUCTION:

The growing importance of deep learning in research, in industry and in our society will require extreme computational power as the dataset sizes and the complexity of these algorithms is expected to increase. Quantum computers are a good candidate to answer this challenge. The recent progress in the physical realization of quantum processors and the advances in quantum algorithms increases the importance of understanding their capabilities and limitations. In particular, the field of quantum machine learning has witnessed many innovative algorithms that offer speedups over their classical counterparts Kerenidis et al. (2019) ; Lloyd et al. (2013; 2014) ; Kerenidis & Prakash (2017b) ; Wiebe et al. (2014a) .

Quantum deep learning refers to the problem of creating quantum circuits that mimic and enhance the operations of neural networks. It has been studied in several works Allcock et al. (2018) ; Rebentrost et al. (2018) ; Wiebe et al. (2014b) but remains challenging as it is difficult to implement non linearities with quantum unitaries Schuld et al. (2014) . In this work we propose a quantum algorithm for convolutional neural networks (CNN), a type of deep learning designed for visual recognition, signal processing and time series. We also provide results of numerical simulations to evaluate the running time and accuracy of the quantum convolutional neural network (QCNN). Note that our algorithm is theoretical and could be compiled on any type of quantum computers (trapped ions, superconducting qubits, cold atoms, photons, etc.)

The CNN was originally developed by LeCun et al. (1998) in the 1980's. They have achieved great practical success over the last decade Krizhevsky et al. (2012) and have been used in cutting-edge domains like autonomous cars Bojarski et al. (2016) and gravitational wave detection George & Huerta (2018) . Despite these successes, CNNs suffer from computational bottlenecks due to the size of the optimization space and the complexity of the inner operations, these bottlenecks make deep CNNs resource expensive.

The growing interest in quantum machine learning has led researchers to develop different variants of Quantum Neural Networks (QNN). The quest for designing quantum analogs of neural networks is challenging due to the modular layer architecture of the neural networks and the presence of non linearities, pooling, and other non unitary operations, as explained in Schuld et al. (2014) . Several strategies have been tried in order to to implement some features of neural networks Allcock et al. (2018) ; Wiebe et al. (2014b) ; Beer et al. (2019) in the quantum setting.

Variational quantum circuits provide another path to the design of QNNs, this approach has been developed in Farhi & Neven (2018) ; Henderson et al. (2019) ; Killoran et al. (2018) . A quantum convolutional neural network architecture using variational circuits was recently proposed Cong et al. (2018) . However further work is required to provide evidence that such techniques can outperform classical neural networks in machine learning settings.

2. PRELIMINARIES 2.1 CONVOLUTION PRODUCT AS MATRIX MULTIPLICATION:

We briefly introduce the formalism and notation concerning classical convolution product and its equivalence with matrix multiplication. More details can be found in Appendix (Section C). A single layer of the classical CNN does the following operations: from an input image X ∈ R H ×W ×D seen as a 3D tensor, and a kernel K ∈ R H×W ×D ×D +1 seen as a 4D tensor, it performs a convolution product and outputs X +1 = X * K , with X +1 ∈ R H +1 ×W +1 ×D +1 . This convolution operation is equivalent to the matrix multiplication A F = Y +1 where A , F and Y +1 are suitably vectorized versions of X , K and X +1 respectively. The output of the layer of the CNN is f (X +1 ) where f is a non linear function.

2.2. QUANTUM COMPUTING:

For a detailed introduction to quantum computing and its applications to machine learning in the context of this work, we invite the reader to look at Appendix F. We also refer to Nielsen & Chuang (2002b) for a more complete overview of quantum computing.

In this part we will discuss only briefly the core notions of quantum computing. Like a classical bit, a quantum bit (qubit) can be |0 , |1 , but can also be in a superposition state α |0 + β |1 with amplitudes (α, β) ∈ C such that |α| 2 + |β| 2 = 1. With n qubits it is then possible to construct a superposition of the 2 n binary combinations possible, each with a specific amplitude. We will note the i th combination (e.g. |01 • • • 110 ) as |i . A vector v ∈ R d can be encoded in a quantum state made of log(d) qubits. This encoding is a quantum superposition, where the components (v 1 , • • • , v d ) of v are used as the amplitudes of the d binary combinations. We note this state |v := 1 v i∈ [d] v i |i , where |i is a register representing the i th vector in the standard basis. Quantum computation proceeds by applying quantum gates which are defined to be unitary matrices acting on 1 or 2 qubits, for example the Hadamard gate that maps

|0 → 1 √ 2 (|0 + |1 ) and |1 → 1 √ 2 (|0 -|1

). The output of the computation is a quantum state that can be measured to obtain classical information. The measurement of a qubit α |0 + β |1 yields either 0 or 1, with probability equal to the square of the respective amplitude. A detailed discussion of the results from quantum machine learning and linear algebra used in this work can be found in Appendix (Section F).

The QCNN can be directly compared to the classical CNN as it has the same inputs and outputs. We show that it offers a speedup compared to certain cases of classical CNN for both the forward pass and for training using backpropagation. For each layer, on the forward pass (Algorithm 1), the speedup is exponential in the size of the layer (number of kernels) and almost quadratic on the spatial dimension of the input. We next state informally the speedup for the forward pass, the formal version appears as Theorem D.1.

Result 1 (Quantum Convolution Layer) Let X be the input and K be the kernel for layer of a convolutional neural network, and f : R → [0, C] with C > 0 be a non linear function so that f (X +1 ) := f (X * K ) is the output for layer . Given X and K stored in QRAM (Quantum Random Access Memory), there is a quantum algorithm that, for precision parameters > 0 and η > 0, creates quantum state |f (X +1 ) such that f (X +1 ) -f (X +1 ) ∞ ≤ 2 and retrieves classical tensor X +1 such that for each pixel j,

|X +1 j -f (X +1 j )| ≤ 2 if f (X +1 j ) ≥ η X +1 j = 0 if f (X +1 j ) < η (1)

The running time of the algorithm is

O 1 η 2 • M √ C √ E(f (X +1 ))

where E(•) represents the average value, O hides factors poly-logarithmic in the size of X and K and the parameter M = max p,q A p F q is the maximum product of norms from subregions of X and K .

We see that the number of elements in the input and the kernels appear only with a poly-logarithmic contribution in the running time. This is one of the main advantages of our algorithm and it allows us to use for larger and even exponentially deeper kernels. For the number of elements in the input, their number is hidden in the precision parameter η in the running time. Indeed, a sufficiently large fraction of pixels must be sampled from the output of the quantum convolution to retrieve the meaningful information. In the Numerical Simulations (Section 6) we provide empirical estimates for η. For details about the QRAM, see Appendix F.2.

Following the forward pass, a loss function L is computed for the output of a classical CNN. The backpropagation algorithm is then used to calculate, layer by layer, the gradient of this loss with respect to the elements of the kernels K , in order to update them through gradient descent. We state our quantum backpropagation algorithm next, the formal version of this result appears as Theorem E.1

Result 2 (Quantum Backpropagation for Quantum CNN) Given the forward pass quantum algorithm in Result 1, and given the kernel matrix F , input matrices A and Y , stored in the QRAM for each layer , and a loss function L, there is a quantum backpropagation algorithm that estimates each element of the gradient tensor ∂L ∂F within additive error δ ∂L ∂F , and updates F according to a gradient descent update rule. The running time of a single layer for quantum backpropagation is given by

O µ(A ) + µ( ∂L ∂Y +1 ) κ( ∂L ∂F ) + µ( ∂L ∂Y +1 ) + µ(F ) κ( ∂L ∂Y ) log 1/δ δ 2 (2)

where for a matrix V ∈ R n×n , κ(V ) is the condition number and µ(V ) ≤ √ n is a matrix dependent parameter defined in Equation (5).

For the quantum back-propagation algorithm, we introduce a quantum tomography algorithm with ∞ norm guarantees, that could be of independent interest. It is exponentially faster than tomography with 2 norm guarantees and is given as Theorem G.1 in Section G. Numerical simulations on classifying the MNIST dataset show that our quantum CNN achieves a similar classification accuracy as the classical CNN.

4. FORWARD PASS FOR QCNN:

The forward pass algorithm for the QCNN implements the quantum analog of a single quantum convolutional layer. It includes a convolution product between an input and a kernel, followed by the application of a non linear function and pooling operations to prepare the next layer's input. We provide an overview of the main ideas of the algorithm here, the complete technical details are given in the Appendix (Section D).

Algorithm 1 QCNN Layer:

Require: Matrix A representing input to layer and kernel matrix F stored in QRAM. Precision parameters and η, a boolean circuit for a non linear function f : R → [0, C]. Ensure: Outputs the data matrix A +1 for the next layer which is the result of the convolution between the input and the kernel, followed by a non linearity and pooling.

1:

Step 1: Quantum Convolution 1.1: Inner product estimation Perform the following mapping, using QRAM queries on rows A p and columns F q , followed by quantum inner product estimation (Theorems F.2 and F.4) to implement the mapping

1 K p,q |p |q → 1 K p,q |p |q |P pq |g pq Where P pq is -close to P pq = 1+ A p |F q 2 , K = √ H +1 W +1 D +1

is a normalisation factor and |g pq is some garbage quantum state that can be ignored. 1.2: Non linearity Use an arithmetic circuit and two QRAM queries to obtain Y +1 , an -approximation of the convolution output Y +1 p,q = (A p , F q ) and apply the non linear function f as a boolean circuit to obtain Update the QRAM for the next layer A +1 while sampling. The implementation of pooling (Max, Average, etc.) can be done by a specific update to the QRAM data structure described in Section D.2.2.

1 K p,q |p |q |f (Y +1 p,q ) |g pq

In this algorithm, we propose the first quantum algorithm for performing the convolution product.

Our algorithm is based on the observation that the convolution product can be regarded as a matrix product between reshaped matrices. The reshaped input's rows A p and the reshaped kernel's columns F q are loaded as quantum states, in superposition. Then the entries of the convolution A p |F q are estimated using a simple quantum circuit for inner product estimation and stored in an auxiliary register as in Step 1.1 of Algorithm 1.

One of the difficulties in the design of quantum neural networks is that non linear functions are hard to implement as unitary operations. We get around this difficulty by applying the non-linear function f as a boolean circuit to the output of the quantum inner product estimation circuit in Step 1.2 of Algorithm 1. Most of the non linear functions in the machine learning literature can be implemented using small sized boolean circuits, our algorithm thus allows for many possible choices of the nonlinear function f (see Appendix F.1 for details on non linear boolean circuits in quantum circuits).

Step 2 of Algorithm 1 develops a quantum importance sampling procedure wherein the pixels with high values of f (Y +1 pq ) are read out with higher probability. This is done by encoding these values into the amplitudes of the quantum state using the well known Amplitude Amplification algorithm Brassard et al. (2002) . This kind of importance sampling is a task that can be performed easily in the quantum setting and has no direct classical analog. Although it does not lead to asymptotic improvements for the algorithms running time, it could lead to improvements that are significant in practice.

More precisely, during the measurement of a quantum register in superposition, only one of its values appears, with a probability corresponding the the square of its amplitude. It implies that the output's pixels measured with more probability are the ones with the highest value f (Y +1 p,q ). Once measured, we read directly from the registers the position p, q and the value itself. Thus we claim that we measure only a fraction of the quantum convolution product output, and that the set of pixels measured collect most of the meaningful information for the CNN, the other pixels being set to 0.

After being measured, each pixel's value and position are stored in a QRAM to be used as quantum state for next layer's input. During this phase, it is possible to discard or aggregate some values to perform pooling operations as described in Step 3 of Algorithm 1. The forward pass for the QCNN thus includes the the convolution product, the non linearity f and pooling operation, in time polylogarithmic in the kernel's dimensions. In comparison, the classical CNN layer in linear in both kernel and input dimensions.

Note finally that quantum importance sampling in Step 2 implies that the non linear function f be bounded by a parameter C > 1. In our experiments we use the capReLu function, which is a modified ReLu function that becomes constant above C.

5. QUANTUM BACKPROPAGATION ALGORITHM:

The second algorithm required for the QCNN is the quantum backpropagation algorithm given as Algorithm 2. Like the classical backpropagation algorithm, it updates all kernels weights according to the derivatives of a given loss function L. In this sectiion, we explain the main ideas and compare it to the classical backpropagation algorithm, the complete details are given in Appendix (Section E).

Algorithm 2 Quantum Backpropagation:

Require: Precision parameter δ. Data matrices A and kernel matrices F stored in QRAM for each layer . Ensure: Outputs gradient matrices ∂L ∂F and ∂L ∂Y for each layer .

1: Calculate the gradient for the last layer L using the outputs and the true labels:

∂L ∂Y L 2: for = L -1, • • • , 0 do 3:

Step 1 : Modify the gradient With ∂L ∂Y +1 stored in QRAM, set to 0 some of its values to take into account pooling, tomography and non linearity that occurred in the forward pass of layer . These values correspond to positions that haven't been sampled nor pooled, since they have no impact on the final loss.

4::

Step 2 : Matrix-matrix multiplications With the modified values of ∂L ∂Y +1 , use quantum linear algebra algorithm (Theorem F.7) to perform the matrix-matrix multiplications (A ) T • ∂L ∂Y +1 and ∂L ∂Y +1 • (F ) T , allowing to obtain quantum states corresponding to ∂L ∂F and ∂L ∂Y .

5::

Step 3 : ∞ tomography Measure the previous outputs, as in Algorithm 3. This allows to estimate each entry of Step 4 : Gradient descent From the previous tomography, perform the gradient descent to update the values of F in QRAM:

F s,q ← F s,q -λ ∂L ∂F s,q ± 2δ ∂L ∂F 2 . 7: end for

We describe briefly detail the implementation of quantum backpropagation at layer . The algorithm assumes that ∂L ∂Y +1 is known. First, the backpropagation of the quantum convolution product is equivalent to the classical one, and we use the matrix-matrix multiplication formulation to obtain the derivatives ∂L ∂F and ∂L ∂Y . The first one is the result wanted and the second one is needed for layer -1. This matrix-matrix multiplication can be implemented as a quantum circuit, by decomposing into several matrix-vector multiplications, known to be efficient, with a running time depending on the ranks and Frobenius norm of the matrices. We obtain a quantum state corresponding to a superposition of all derivatives. We use again the ∞ tomography to retrieve each derivative with precision δ > 0 such that, for all kernel's weight F s,q we have approximated it's loss derivative with ∂L ∂F s,q

, with an error bounded by ∂L ∂F s,q -∂L ∂F s,q ≤ 2δ ∂L ∂F 2 . This implies that the gradient descent rule is perturbed by 2δ ∂L ∂F 2 at most, see Appendix (Section E.4). We also take into account the effects of quantum non linearity, quantum measurement and pooling. The quantum pooling operation is equivalent to the classical one, where pixels that were not selected during pooling see their derivative set to 0. Quantum measurement is similar, since pixels that haven't been measured don't contribute to the gradient. For the non linearity, as in the classical case, pixels with negative values were set to zero, hence should have no contribution to the gradient. Additionally, because we used the capReLu function, pixels bigger than the threshold C must also have null derivatives. This two rules can be implemented by combining them with measurement rules compared to classical backpropagation, see Appendix (Section E.2.2) for details.

6. NUMERICAL SIMULATIONS:

As described above, the adaptation of the CNNs to the quantum setting implies some modifications that could alter the efficiency of the learning or classifying phases. We now present some experiments to show that such modified CNNs can converge correctly, as the original ones.

The experiment, using the PyTorch library developed by Paszke et al. (2017) , consists of training classically a small convolutional neural network for which we have added a "quantum" sampling after each convolution. Instead of parametrising it with the precision η, we have choosed to use the sampling ratio σ that represents the fraction of pixels drawn during tomography. This two definitions are equivalent, as shown in Appendix (Section D.1.5), but the second one is more intuitive regarding the running time and the simulations.

We also add a noise simulating the amplitude estimation (parameter ), followed by a capReLu instead of the usual ReLu (parameter C), and a noise during the backpropagation (parameter δ). In the following results, we observe that our quantum CNN is able to learn and classify visual data from the widely used MNIST dataset. This dataset is made of 60.000 training images and 10.000 testing images of handwritten digits. Each image is a 28x28 grayscale pixels between 0 and 255 (8 bits encoding), before normalization.

Let's first observe the "quantum" effects on an image of the dataset. In particular, the effect of the capped non linearity, the introduction of noise and the quantum sampling.

We now present the full simulation of our quantum CNN. In the following, we use a simple network made of 2 convolution layers, and compare our quantum CNN to the classical one. The first and second layers are respectively made of 5 and 10 kernels, both of size 7x7. A three-layer fully connected network is applied at the end and a softmax activation function is applied on the last layer to detect the predicted outcome over 10 classes (the ten possible digits). Note that we didn't introduce pooling, being equivalent between quantum and classical algorithms and not improving the results on our CNN. The objective of the learning phase is to minimize the loss function, defined by the negative log likelihood of the classification on the training set. The optimizer used was a built-in Stochastic Gradient Descent.

Using PyTorch, we have been able to implement the following quantum effects (the first three points are shown in Figure 1 ):

-The addition of a noise, to simulate the approximation of amplitude estimation during the forward quantum convolution layer, by adding gaussian noise centered on 0 and with standard deviation 2M , with M = max p,q A p F q .

-A modification of the non linearity: a ReLu function which is constant above the value T (the cap).

-A sampling procedure to apply on a tensor with a probability distribution proportional to the tensor itself, reproducing the quantum sampling with ratio σ.

-The addition of a noise during the gradient descent, to simulate the quantum backpropagation, by adding a gaussian noise centered on 0 with standard deviation δ, multiplied by the norm of the gradient, as given by Equation (28).

Figure 1 : Effects of the QCNN on a 28x28 input image. From left to right: original image, image after applying a capReLu activation function with a cap C at 2.0, introduction of a strong noise during amplitude estimation with = 0.5, quantum sampling with ratio σ = 0.4 that samples the highest values in priority. The useful information tends to be conserved in this example. The side gray scale indicates the value of each pixel. Note that during the QCNN layer, a convolution is supposed to happen before the last image but we chose not to perform it for visualisation matter.

The CNN used for this simulation may seem "small" compared to the standards AlexNet developed by Krizhevsky et al. (2012) or VGG-16 by Simonyan & Zisserman (2014) , or those used in industry. However simulating this small QCNN on a classical computer was already very computationally intensive and time consuming, due to the"quantum" sampling task, apparently not optimized for a classical implementation in PyTorch. Every single training curve showed in Figure 9 could last for 4 to 8 hours. Hence adding more convolutional layers wasn't convenient. Similarly, we didn't compute the loss on the whole testing set (10.000 images) during the training to plot the testing curve. However we have computed the test losses and accuracies once the model trained (see Table 4 ), in order to detect potential overfitting cases.

We now present the result of the training phase for a quantum version of this CNN, where partial quantum sampling is applied, for different sampling ratio (number of samples taken from the resulting convolution). Since the quantum sampling gives more probability to observe high value pixels, we expect to be able to learn correctly even with small ratio (σ ≤ 0.5). We compare these training curve to the classical one. The learning has been done on two epochs, meaning that the whole dataset is used twice. The following plots show the evolution of the loss L during the iterations on batches. This is the standard indicator of the good convergence of a neural network learning phase. We can compare the evolution of the loss between a classical CNN and our QCNN for different parameters. Most results are presented in Appendix (Section H).

Our simulations show that the QCNN is able to learn despite the introduction of noise, tensor sampling and other modifications. In particular it shows that only a fraction of the information is meaningful for the neural network, and that the quantum algorithm captures this information in priority. This learning can be more or less efficient depending on the choice of the key parameters. For decent values of these parameters, the QCNN is able to converge during the training phase. It can then classify correctly on both training and testing set, indicating neither overfitting nor underfitting.

We notice that the learning curves sometimes present a late start before the convergence initializes, in particular for small sampling ratio. This late start can be due to the random initialization of the kernel weights, that performs a meaningless convolution, a case where the quantum sampling of the output is of no interest. However it is very interesting to see that despite this late start, the kernel start converging once they have found a good combination.

Overall, it is possible that the QCNN presents some behaviors that have no classical equivalence. Understanding their potential effects, positive or negative, is an open question, all the more so as the effects of the classical CNN's hyperparameters are already a topic of active research, see the work of Samek et al. (2017) for details. Note also that the neural network used in this simulation is Figure 2 : Training curves comparison between the classical CNN and the Quantum CNN (QCNN) for = 0.01, C = 10, δ = 0.01 and the sampling ratio σ from 0.1 to 0.5. We can observe a learning phase similar to the classical one, even for a weak sampling of 20% or 30% of each convolution output, which tends to show that the meaningful information is distributed only at certain location of the images, coherently with the purpose of the convolution layer. Even for a very low sampling ratio of 10%, we observe a convergence despite a late start. rather small. A following experiment would be to simulate a quantum version of a standard deeper CNN (AlexNet or VGG-16), eventually on more complex dataset, such as CIFAR-10 developed by Krizhevsky & Hinton (2009) or Fashion MNIST by Xiao et al. (2017) .

7. CONCLUSIONS:

We have presented a quantum algorithm for evaluating and training convolutional neural networks (CNN). At the core of this algorithm, we have developed a novel quantum algorithm for computing a convolution product between two tensors, with a substantial speed up. This technique could be reused in other signal processing tasks that could benefit from an enhancement by a quantum computer. Layer by layer, convolutional neural networks process and extract meaningful information. Following this idea of learning foremost important features, we have proposed a new approach of quantum tomography where the most meaningful information is sampled with higher probability, hence reducing the complexity of our algorithm.

Our QCNN is complete in the sense that almost all classical architectures can be implemented in a quantum fashion: any (non negative and upper bounded) non linearity, pooling, number of layers and size of kernels are available. Our circuit is shallow and could be run on relatively small quantum computers. One could repeat the main loop many times on the same shallow circuit, since performing the convolution product is simple, and is similar for all layer. The pooling and non linearity are included in the loop. Our building block approach, layer by layer, allows high modularity, and can be combined with work on quantum feedforward neural network developed by Allcock et al. (2018) .

The running time presents a speedup compared to the classical algorithm, due to fast linear algebra when computing the convolution product, and by only sampling the important values from the resulting quantum state. This speedup can be highly significant in cases where the number of channels D in the input tensor is high (high dimensional time series, videos sequences, games play) or when the number of kernels D +1 is big, allowing deep architectures for CNN, which was the case in the recent breakthrough of DeepMind AlphaGo algorithm of Silver et al. (2016) . The Quantum CNN also allows larger kernels, that could be used for larger input images, since the size the kernels must be a contant fraction of the input in order to recognize patterns. However, despite our new techniques to reduce the complexity, applying a non linearity and reusing the result of a layer for the next layer make register encoding and state tomography mandatory, hence preventing from having an exponential speedup on the number of input parameters.

Finally we have presented a backpropagation algorithm that can also be implemented as a quantum circuit. The numerical simulations on a small CNN show that despite the introduction of noise and sampling, the QCNN can efficiently learn to classify visual data from the MNIST dataset, performing a similar accuracy than the classical CNN.

APPENDIX A VARIABLE SUMMARY:

We recall the most important variables for layer . They represent tensors, their approximations, and their reshaped versions.

Data:

Variable Dimensions Indices

Input

X H × W × D (i , j , d ) Y (H W ) × D - A (H +1 W +1 ) × (HW D ) (p, r) Kernel K H × W × D × D +1 (i, j, d, d ) F (HW D ) × D +1 (s, q)

Table 1: Summary of input variables for the th layer, along with their meaning, dimensions and corresponding notations. These variables are common for both quantum and classical algorithms.

We have omitted indices for Y which don't appear in our work.

Data Variable Dimensions Indices Output of Quantum Convolution f (Y +1 ) (H +1 W +1 ) × D +1 (p, q) f (X +1 ) H +1 × W +1 × D +1 (i +1 , j +1 , d +1 ) Output of Quantum Tomography X +1 H +1 × W +1 × D +1 (i +1 , j +1 , d +1 ) Output of Quantum Pooling X +1 H +1 P × W +1 P × D +1 ( ĩ +1 , j +1 , d +1 )

Table 2 : Summary of variables describing outputs of the layer , with the quantum algorithm.

Data Variable Dimensions Indices:

Output of Classical Convolution

f (Y +1 ) (H +1 W +1 ) × D +1 (p, q) f (X +1 ) H +1 × W +1 × D +1 (i +1 , j +1 , d +1 ) Output of Classical Pooling X +1 H +1 P × W +1 P × D +1 ( ĩ +1 , j +1 , d +1 )

Table 3 : Summary of variables describing outputs of the layer , with the classical algorithm.

Classical and quantum algorithms can be compared with these two diagrams:

Quantum convolution layer :

EQUATION

We finally provide some remarks that could clarify some notations ambiguity:

-Formally, the output of the quantum algorithm is X +1 . It is used as input for the next layer + 1.

But we consider that all variables' names are reset when starting a new layer: X +1 ← X +1 .

-For simplicity, we have sometimes replaced the indices (i +1 , j +1 , d +1 ) by n to index the elements of the output.

-In Section D.2.2, the input for layer + 1 is stored as A +1 , for which the elements are indexed by (p , r ).

APPENDIX B PRELIMINARIES IN QUANTUM INFORMATION:

We introduce a basic and broad-audience quantum information background necessary for this work.

For a more detailed introduction we recommend Nielsen & Chuang (2002a).

B.1 QUANTUM INFORMATION:

Quantum Bits and Quantum Registers: The bit is the most basic unit of classical information.

It can be either in state 0 or 1. Similarly a quantum bit or qubit, is a quantum system that can be is state |0 , |1 (the braket notation |• is a reminder that the bit considered is a quantum system) or in superposition of both states α |0 + β |1 with coefficients α, β ∈ C such that |α| 2 + |β| 2 = 1. The amplitudes α and β are linked to the probabilities of observing either 0 or 1 when measuring the qubit, since P (0) = |α| 2 and P (1) = |β| 2 .

Before the measurement, any superposition is possible, which gives quantum information special abilities in terms of computation. With n qubits, the 2 n possible binary combinations can exist simultaneously, each with a specific amplitude. For instance we can consider an uniform distribution

1 √ n 2 n -1

i=0 |i where |i represents the i th binary combination (e.g. |01 • • • 1001 ). Multiple qubits together are often called a quantum register.

In its most general formulation, a quantum state with n qubits can be seen as vector in a complex Hilbert space of dimension 2 n . This vector must be normalized under 2 -norm, to guarantee that the squared amplitudes sum to 1.

Quantum Computation: To process qubits and therefore quantum registers, we use quantum gates. These gates are unitary operators in the Hilbert space as they should map unit-norm vectors to unit-norm vectors. Formally, we can see a quantum gate acting on n qubits as a matrix U ∈ C 2 n such that U U † = U † U = I, where U † is the conjugate transpose of U . Some basic single qubit gates includes the NOT gate 0 1 1 0 that inverts |0 and |1 , or the Hadamard gate 1

√ 2 1 1 1 -1 that maps |0 → 1 √ 2 (|0 + |1 ) and |1 → 1 √ 2 (|0 -|1

), creating the quantum superposition. Finally, multiple qubits gates exist, such as the Controlled-NOT that applies a NOT gate on a target qubit conditioned on the state of a control qubit.

The main advantage of quantum gates is their ability to be applied to a superposition of inputs. Indeed, given a gate U such that U |x → |f (x) , we can apply it to all possible combinations of x at once U ( 1

C x |x ) → 1 C x |f (x) .

We now state some primitive quantum circuits, which we will use in our algorithm: For two integers i and j, we can check their equality with the mapping |i |j |0 → |i |j |[i = j] . For two real value numbers a > 0 and δ > 0, we can compare them using |a |δ |0 → |a |δ |[a ≤ δ] . Finally, for a real value numbers a > 0, we can obtain its square |a |0 → |a |a 2 . Note that these circuits are basically a reversible version of the classical ones and are linear in the number of qubits used to encode the input values.

Any classical boolean function can be implemented in a quantum unitary, even though this seems at first contradictory with the requirements of unitaries (reversibility, linearity). Let σ : R → R be a classical function, we define U σ the unitary that acts as U σ |x |0 → |x |σ(x) . Using a second quantum register to encode the result of the function, the properties of quantum unitaries are respected.

B.2 QUANTUM SUBROUTINES FOR DATA ENCODING:

Knowing some basic principles of quantum information, the next step is to understand how data can be efficiently encoded using quantum states. While several approaches could exist, we present the most common one called amplitude encoding, which leads to interesting and efficient applications.

Let x ∈ R d be a vector with components (x 1 , • • • , x d ). Using only log(d) qubits, we can form |x , the quantum state encoding x, given by |x = 1 x d-1 j=0 x j |j . We see that the j th component x j becomes the amplitude of |j , the j th binary combination (or equivalently the j th vector in the standard basis). Each amplitude must be divided by x to preserve the unit 2 -norm of |x .

Similarly, for a matrix A ∈ R n×d or equivalently for n vectors A i for i ∈ [n], we can express each row of A as

|A i = 1 Ai d-1 i=0 A ij |j .

We can now explain an important definition, the ability to have quantum access to a matrix. This will be a requirements for many algorithms.

Definition 1 [Quantum Access to Data] We say that we have quantum access to a matrix A ∈ R n×d if there exist a procedure to perform the following mapping, for i ∈ [n], in time T :

• |i |0 → |i |A i • |0 → 1 A F i A i |i

By using appropriate data structures the first mapping can be reduced to the ability to perform a mapping of the form |i |j |0 → |i |j |A ij . The second requirement can be replaced by the ability of performing |i |0 → |i | A i or to just have the knowledge of each norm. Therefore, using matrices such that all rows A i have the same norm makes it simpler to obtain the quantum access.

The time or complexity T necessary for the quantum access can be reduced to polylogarithmic dependence in n and d if we consider the access to a Quantum Memory or QRAM. The QRAM Kerenidis & Prakash (2017a) is a specific data structure from which a quantum circuit can allow quantum access to data in time O(log (nd)).

Theorem B.1 (QRAM data structure, see Kerenidis & Prakash (2017a) ) Let A ∈ R n×d , there is a data structure to store the rows of A such that, 1. The time to insert, update or delete a single entry

A ij is O(log 2 (n)).

2. A quantum algorithm with access to the data structure can perform the following unitaries in time

T = O(log 2 n). (a) |i |0 → |i |A i for i ∈ [n]. (b) |0 → i∈[n] A i |i .

We now state important methods for processing the quantum information. Their goal is to store some information alternatively in the quantum state's amplitude or in the quantum register as a bitstring.

Theorem B.2 [Amplitude Amplification and Estimation Brassard et al. (2002) ] Given a unitary operator U such that U : |0 → √ p |y |0 + √ 1 -p |y ⊥ |1 in time T , where p > 0 is the probability of measuring "0", it is possible to obtain the state |y |0 using O( T √ p ) queries to U , or to estimate p with relative error δ using O( T δ √ p ) queries to U .

Theorem B.3 [Conditional Rotation] Given the quantum state |a , with a ∈ [-1, 1], it is possible to perform |a |0 → |a (a |0 + √ 1 -a |1 ) with complexity O(1).

Using Theorem F.3 followed by Theorem F.2, it then possible to transform the state 1

√ d d-1 j=0 |x j into 1 x d-1 j=0 x j |x j .

In addition to amplitude estimation, we will make use of a tool developed in Wiebe et al. (2014a) to boost the probability of getting a good estimate for the inner product required for the quantum convolution algorithm. In high level, we take multiple copies of the estimator from the amplitude estimation procedure, compute the median, and reverse the circuit to get rid of the garbage. Here we provide a theorem with respect to time and not query complexity.

Theorem B.4 (Median Evaluation, see Wiebe et al. (2014a) ) Let U be a unitary operation that maps

U : |0 ⊗n → √ a |x, 1 + √ 1 -a |G, 0 for some 1/2 < a ≤ 1 in time T .

Then there exists a quantum algorithm that, for any ∆ > 0 and for any 1/2 < a 0 ≤ a, produces a state |Ψ such that |Ψ -|0 ⊗nL |x ≤ √ 2∆ for some integer

L, in time 2T ln(1/∆) 2 |a 0 | -1 2 2 . B.3 QUANTUM SUBROUTINES FOR LINEAR ALGEBRA

In the recent years, as the field of quantum machine learning grew, its "toolkit" for linear algebra algorithms has become important enough to allow the development of many quantum machine learning algorithms. We introduce here the important subroutines for this work, without detailing the circuits or the algorithms.

Definition 2 For a matrix A, the parameter µ(A) is defined by

µ(A) = min p∈[0,1] A F , s 2p (A)s 2(1-p) (A T ) where s p (A) = max i ( A i p p ).

The next theorems allow to compute the distance between vectors encoded as quantum states, and use this idea to perform the k-means algorithm. Given quantum access in time T to two matrices U and V with rows u i and v j of dimension d, there is a quantum algorithm that, for any pair (i, j), performs the following mapping |i |j |0 → |i |j |d 2 (u i , v j ) , estimating the euclidean distance between u i and v j with precision

|d 2 (u i , v j ) -d 2 (u i , v j )| ≤ for any > 0.

The algorithm has a running time given by O(T η/ ),

where η = max ij ( u i v j ), assuming that min i ( u i ) = min i ( v i ) = 1. Theorem B.6 [Quantum k-means clustering Kerenidis et al. (2019)]

Given quantum access in time T to a dataset V ∈ R n×d , there is a quantum algorithm that outputs with high probability k centroids c 1 , • • • , c k that are consistent with the output of the k-means algorithm with noise δ > 0, in time

O(T × (kd η(V ) δ 2 κ(V )(µ(V ) + k η(V ) δ ) + k 2 η(V ) 1.5 δ 2 κ(V )µ(V ))) per iteration.

Definition 3 For a matrix V ∈ R n×d , its parameter η(V ) is defined as as maxi( vi 2 ) mini( vi 2 ) , or as

max i ( v i 2 ) assuming min i ( v i ) = 1.

In theorem F.6, the other parameters in the running time can be interpreted as follows : δ is the precision in the estimation of the distances, but also in the estimation of the position of the centroids. κ(V ) is the condition number of V and µ(V ) is defined above (Definition 5). Finally, in the case of well clusterable datasets, which should be the case when we will apply k-means during spectral clustering, the running simplifies to

O(T × (k 2 d η(V ) 2.5 δ 3 + k 2.5 η(V ) 2 δ 3 )).

Note that the dependence in n is hidden in the time T to load the data. This dependence becomes polylogarithmic in n if we assume access to a QRAM.

Theorem B.7 (Quantum Matrix Operations, Chakraborty et al. (2018) ) Let M ∈ R d×d and x ∈ R d . Let δ 1 , δ 2 > 0.

If M is stored in appropriate QRAM data structures and the time to prepare |x is T x , then there exist quantum algorithms that with probability at least 1 -1/poly(d) return

1. A state |z such that |z -|M x 2 ≤ δ 1 in time O((κ(M )µ(M ) + T x κ(M )) log(1/δ 1 )).

Note that this also implies

|z -|M x ∞ ≤ δ 1 2. Norm estimate z ∈ (1 ± δ 2 ) M x 2 , with relative error δ 2 , in time O(T x κ(M )µ(M ) δ2

log(1/δ 1 )).

The linear algebra procedures above can also be applied to any rectangular matrix V ∈ R n×d by considering instead the symmetric matrix V = 0 V V T 0 .

APPENDIX C CLASSICAL CONVOLUTIONAL NEURAL NETWORK (CNN):

CNN is a specific type of neural network, designed in particular for image processing or time series. It uses the Convolution Product as a main procedure for each layer. We will focus on image processing with a tensor framework for all elements of the network. Our goal is to explicitly describe the CNN procedures in a form that can be translated in the context of quantum algorithms.

As a regular neural network, a CNN should learn how to classify any input, in our case images. The training consists of optimizing a series of parameters, learned on the inputs and their corresponding labels.

C.1 TENSOR REPRESENTATION Images, or more generally layers of the network, can be seen as tensors. A tensor is a generalization of a matrix to higher dimensions. For instance an image of height H and width W can be seen as a matrix in R H×W , where every pixel is a greyscale value between 0 ans 255 (8 bit). However the three channels of color (RGB: Red Green Blue) must be taken into account, by stacking three times the matrix for each color. The whole image is then seen as a 3 dimensional tensor in R H×W ×D where D is the number of channels. We will see that the Convolution Product in the CNN can be expressed between 3-tensors (input) and 4-tensors (convolution filters or kernels), the output being a 3-tensor of different dimensions (spatial size and number of channels).

Figure 3 : RGB decomposition, a colored image is a 3-tensor.

C.2 ARCHITECTURE:

A CNN is composed of 4 main procedures, compiled and repeated in any order : Convolution layers, most often followed by an Activation Function, Pooling Layers and some Fully Connected layers at the end. We will note the current layer.

Convolution Layer : The th layer is convolved by a set of filters called kernels. The output of this operation is the ( + 1) th layer. A convolution by a single kernel can be seen as a feature detector, that will screen over all regions of the input. If the feature represented by the kernel, for instance a vertical edge, is present in some part of the input, there will be a high value at the corresponding position of the output. The output is commonly called the feature map of this convolution.

Activation Function : As in regular neural network, we insert some non linearities also called activation functions. These are mandatory for a neural network to be able to learn any function. In the case of a CNN, each convolution is often followed by a Rectified Linear Unit function, or ReLu. This is a simple function that puts all negative values of the output to zero, and lets the positive values as they are.

Pooling Layer : This downsampling technique reduces the dimensionality of the layer, in order to improve the computation. Moreover, it gives to the CNN the ability to learn a representation invariant to small translations. Most of the time, we apply a Maximum Pooling or an Average Pooling. The first one consists of replacing a subregion of P × P elements only by the one with the maximum value. The second does the same by averaging all values. Recall that the value of a pixel corresponds to how much a particular feature was present in the previous convolution layer.

Fully Connected Layer : After a certain number of convolution layers, the input has been sufficiently processed so that we can apply a fully connected network. Weights connect each input to each output, where inputs are all element of the previous layer. The last layer should have one node per possible label. Each node value can be interpreted as the probability of the initial image to belong to the corresponding class.

C.3 CONVOLUTION PRODUCT AS A TENSOR OPERATION:

Most of the following mathematical formulations have been very well detailed by Wu (2017) .

At layer , we consider the convolution of a multiple channels image, seen as a 3-tensor X ∈ R H ×W ×D . Let's consider a single kernel in R H×W ×D . Note that its third dimension must match the number of channels of the input, as in Figure 4 . The kernel passes over all possible regions of the input and outputs a value for each region, stored in the corresponding element of the output. Therefore the output is 2 dimensional, in R H +1 ×W +1 In a CNN, the most general case is to apply several convolution products to the input, each one with a different 3-tensor kernel. Let's consider an input convolved by D +1 kernels. We can globally see this process as a whole, represented by one 4-tensor kernel K ∈ R H×W ×D ×D +1 . As D +1 convolutions are applied, there are D +1 outputs of 2 dimensions, equivalent to a 3-tensor

X +1 ∈ R H +1 ×W +1 ×D +1

We can see on Figure 5 that the output's dimensions are modified given the following rule:

H +1 = H -H + 1 W +1 = W -W + 1 (4)

We omit to detail the use of Padding and Stride, two parameters that control how the kernel moves through the input, but these can easily be incorporated in the algorithms.

An element of X is determined by 3 indices (i , j , d ), while an element of the kernel K is determined by 4 indices (i, j, d, d ). For an element of X +1 we use 3 indices (i +1 , j +1 , d +1 ).

We can express the value of each element of the output X +1 with the relation The convolution operation X * K = X +1 is equivalent to the following matrix multiplication

EQUATION

EQUATION

where each column of Y +1 ∈ R (H +1 W +1 )×D +1 is a column first vectorized form of one of the D +1 channels of X +1 . Note that an element Y +1 p,q is the inner product between the p th row of A and the q th column of F . It is then simple to convert Y +1 into X +1 The indices relation between the elements Y +1 p,q and X +1 i +1 ,j +1 ,d +1 is given by:   

d +1 = q j +1 = p H +1 i +1 = p -H +1 p H +1 (7)

A summary of all variables along with their meaning and dimensions is given in Section A.

APPENDIX D QUANTUM CONVOLUTIONAL NEURAL NETWORK:

In this section we will design quantum procedures for the usual operations in a CNN layer. We start by describing the main ideas before providing the details. Steps are gathered in Algorithm 1.

First, to perform a convolution product between an input and a kernel, we will use the mapping between convolution of tensors and matrix multiplication from Section C.3, that can further be reduced to inner product estimation between vectors, in order to use quantum linear algebra procedures to perform these computations faster. The output will be a quantum state representing the result of the convolution product, from which we can sample to retrieve classical information to feed the next layer. This is stated by the following Theorem: Theorem D.1 (Quantum Convolution Layer) Given 3D tensor input X ∈ R H ×W ×D and 4D tensor kernel K ∈ R H×W ×D ×D +1 stored in QRAM, there is a quantum algorithm that computes a quantum states ∆-close to |f (X +1 ) with arbitrary small parameter ∆ > 0. |f (X +1 ) is close to the result of the convolution product X +1 = X * K followed by any non linear function f : R → R + , with an error bounded by f (X +1 ) -f (X +1 ) ∞ ≤ 2M for any precision > 0, where M is the maximum norm of a product between one of the D +1 kernels, and one of the regions of X of size HW D . The time complexity of this procedure is given by O (1/ ), where O hides factors poly-logarithmic in ∆ and in the size of X and K .

In a second step, we efficiently retrieve classical information from the output. Recall that a convolution can be seen as a pattern detection on the input image, where the pattern is the kernel. The output values correspond to "how much" the pattern was present in the corresponding region of the input. Low value pixels in the output indicate the absence of the pattern in the input at the corresponding regions. Therefore, by sampling according to these output values, where the high value pixels are sampled with more probability, we could retrieve less but only meaningful information for the neural network to learn. While sampling, we update the QRAM data structure with the new information (see Section D.2). We also perform the Pooling operation during this phase (see Section D.2.2). It is an interesting use case where amplitudes of a quantum state are proportional to the importance of the information they carry, giving a new utility to the probabilistic nature of quantum sampling. Numerical simulations are presented in Section 6 to have an empirical estimation of how many samples from the output state are necessary.

D.1 SINGLE QUANTUM CONVOLUTION LAYER:

In order to develop a quantum algorithm to perform the convolution as described above, we will make use of quantum linear algebra procedures. We will use quantum states proportional to the rows of A , noted |A p , and the columns of F , noted |F q (we omit the exponent in the quantum states to simplify the notation). These states are given by |A p = 1

Ap HW D -1 r=0:

A pr |r and |F q = 1 Fq D +1 -1 s=0

F sq |s . We suppose we can load these vectors in quantum states by performing the following queries:

|p |0 → |p |A p |q |0 → |q |F q (8)

Such queries, in time poly-logarithmic in the dimension of the vector, can be implemented with a Quantum Random Access Memory (QRAM). See Section D.2 for more details on the QRAM update rules and its integration layer by layer.

D.1.1 INNER PRODUCT ESTIMATION:

The following method to estimates inner product is derived from previous work by Kerenidis et al. (2019) . With the initial state |p |q 1 √ 2 (|0 + |1 ) |0 we apply the queries detailed above in a controlled fashion, followed simply by a Hadamard gate to extract the inner product A p |F q in an amplitude. 

0 (p, q) = Ppq H +1 W +1 D +1 = 1+ Ap|Fq 2H +1 W +1 D +1

Now we can relate to the Convolution product. Indeed, the triplets (p, q, 0) that are the most probable to be measured are the ones for which the value A p |F q is the highest. Recall that each element of Y +1 is given by Y +1 pq = (A p , F q ), where "(•, •)" denotes the inner product. We see here that we will sample most probably the positions (p, q) for the highest values of Y +1 , that corresponds to the most important points of X +1 , by the Equation ( 7). Note that the the values of Y +1 can be either positive of negative, which is not an issue thanks to the positiveness of P 0 (p, q).

A first approach could be to measure indices (p, q) and rely on the fact that pixels with high values, hence a high amplitude, would have a higher probability to be measured. However we have not exactly the final result, since A p |F q = (A p , F q ) = A p F q A p |F q . Most importantly we then want to apply a non linearity f (Y +1 pq ) to each pixel, for instance the ReLu function, which seems not possible with unitary quantum gates if the data is encoded in the amplitudes only. Morever, due to normalization of the quantum amplitudes and the high dimension of the Hilbert space of the input, the probability of measuring each pixel is roughly the same, making the sampling inefficient. Given these facts, we have added steps to the circuit, in order to measure (p, q, f (Y +1 pq )), therefore know the value of a pixel when measuring it, while still measuring the most important points in priority.

D.1.2 ENCODING THE AMPLITUDE IN A REGISTER:

Let U be the unitary that map |0 to |u

: |u = 1 √ H +1 W +1 D +1

p,q |p |q P pq |0, y pq + 1 -P pq |1, y pq . The amplitude P pq can be encoded in an ancillary register by using Amplitude Estimation (Theorem F.2) followed by a Median Evaluation (Theorem F.4). For any ∆ > 0 and > 0, we can have a state ∆-close to |u =

1 √ H +1 W +1 D +1

p,q |p |q |0 |P pq |g pq with probability at least 1 -2∆, where |P pq -P pq | ≤ and |g pq is a garbage state. This requires O( ln(1/∆) ) queries of U. In the following we discard the third register |0 for simplicity.

The benefit of having P pq in a register is to be able to perform operations on it (arithmetic or even non linear). Therefore we can simply obtain a state corresponding to the exact value of the the convolution product. Since we've built a circuit such that P pq = 1+ Ap|Fq 2

, with two QRAM calls, we can retrieve the norm of the vectors by applying the following unitary |p |q |P pq |g pq |0 |0 → |p |q |P pq |g pq | A p | F q . On the fourth register, we can then write Y +1 pq = A p F q A p |F q using some arithmetic circuits (addition, multiplication by a scalar, multiplication between registers). We then apply a boolean circuit that implements the ReLu function on the same register, in order to obtain an estimate of f (Y +1 pq ) in the fourth register. We finish by inverting the previous computations and obtain the final state:

|f (Y +1 ) = 1 √ H +1 W +1 D +1 p,q |p |q |f (Y +1 pq ) |g pq (9)

Because of the precision on |P pq , our estimation

Y +1 pq = (2P pq -1) A p F q , is obtained with error such that |Y +1 pq -Y +1 pq | ≤ 2 A p F q .

In superposition, we can bound this error by

|Y +1 pq -Y +1 pq | ≤ 2M

where we define

M = max p,q A p F q (10)

M is the maximum product between norms of one of the D +1 kernels, and one of the regions of X of size HW D . Finally, since the previous error estimation is valid for all pairs (p, q), the overall error committed on the convolution product can be bounded by Y

+1 -Y +1 ∞ ≤ 2M , where

. ∞ denotes the ∞ norm. Recall that Y +1 is just a reshaped version of X +1 . Since the non linearity adds no approximation, we can conclude on the final error committed for a layer of our QCNN

EQUATION

At this point, we have established Theorem D.1 as we have created the quantum state ( 9), with given precision guarantees, in time poly-logarithmic in ∆ and in the size of X and K .

We know aim to retrieve classical information from this quantum state. Note that |Y +1 pq is representing a scalar encoded in as many qubits as needed for the precision, whereas |A p was representing a vector as a quantum state in superposition, where each element A p,r is encoded in one amplitude (See Section F). The next step can be seen as a way to retrieve both encoding at the same time, that will allow an efficient tomography focus on the values of high magnitude.

D.1.3 CONDITIONAL ROTATION:

In the following sections, we omit the + 1 exponent for simplicity. Garbage states are removed as they will not perturb the final measurement. We now aim to modify the amplitudes, such that the highest values of |f (Y ) are measured with higher probability. A way to do so consists in applying a conditional rotation on an ancillary qubit, proportionally to f (Y pq ). We will detail the calculation since in the general case f (Y pq ) can be greater than 1. To simplify the notation, we note x = f (Y pq ). This step consists in applying the following rotation on a ancillary qubit:

|x |0 → |x x max x |0 + β |1 , where max x = max p,q f (Y pq ) and β = 1 -( x max x ) 2 .

Note that in practice it is not possible to have access to |max x from the state ( 9), but we will present a method to know a priori this value or an upper bound in section D.1.6. Let's note

α pq = f (Y pq ) maxp,q(f (Y pq ))

. The ouput of this conditional rotation in superposition on state ( 9) is then

1 √ HW D p,q |p |q |f (Y pq ) (α pq |0 + 1 -α 2 pq |1 ).

D.1.4 AMPLITUDE AMPLIFICATION:

In order to measure (p, q, f (Y pq )) with higher probability where f (Y pq ) has high value, we could post select on the measurement of |0 on the last register. Otherwise, we can perform an amplitude amplification on this ancillary qubit. Let's rewrite the previous state as ), as shown by Brassard et al. (2002) . Since f (Y pq ) ∈ R + , we have

1 √ HW D p,q α pq |p |q |f (Y pq ) |0 + 1 -α 2 pq |g pq |1 ,

α 2 pq = f (Y pq ) maxp,q(f (Y pq ))

. Therefore the number of queries is O max p,q (f (Y pq ))

1 1 HW D p,q f (Y pq ) = O √ maxp,q(f (Y pq )) √ Ep,q(f (Y pq ))

, where the notation E p,q (f (Y pq )) represents the average value of the matrix f (Y ). It can also be written E(f (X)) as in Result 1: E p,q (f (Y pq )) = 1 HW D p,q f (Y pq ). At the end of these iterations, we have modified with high probability the state to the following:

EQUATION

Where, to respect the normalization of the quantum state, α pq = αpq p,q α 2 pq HW D

. Eventually, the probability of measuring (p, q, f (Y pq )) is given by p(p, q, f (Y pq )) =

(α pq ) 2 HW D = f (Y pq ) p,q f (Y pq )

. Note that we have used the same type of name |f (Y ) for both state (9) and state (12) . For now on, this state name will refer only to the latter (12).

D.1.5 ∞ TOMOGRAPHY AND PROBABILISTIC SAMPLING:

We can rewrite the final quantum state obtained in (12) as

EQUATION

We see here that f (Y +1 pq ), the values of each pixel, are encoded in both the last register and in the amplitude. We will use this property to extract efficiently the exact values of high magnitude pixels. For simplicity, we will use instead the notation f (X The pixels with high values will have more probability of being sampled. Specifically, we perform a tomography with ∞ guarantee and precision parameter η > 0. See Theorem G.1 and Section G for details. The ∞ guarantee allows to obtain each pixel with error at most η, and require O(1/η 2 ) samples from the state (13). Pixels with low values f (X +1 n ) < η will probably not be sampled due to their low amplitude. Therefore the error committed will be significative and we adopt the rule of setting them to 0. Pixels with higher values f (X +1 n ) ≥ η, will be sample with high probability, and only one appearance is enough to get the exact register value f (X +1 n ) of the pixel, as is it also written in the last register.

To conclude, let's note X +1 n the resulting pixel values after the tomography, and compare it to the real classical outputs f (X +1 n ). Recall that the measured values f (X +1 n ) are approximated with error at most 2M with M = max p,q A p F q . The algorithm described above implements the following rules:

|X +1 n -f (X +1 n )| ≤ 2M if f (X +1 n ) ≥ η X +1 n = 0 if f (X +1 n ) < η (14)

Concerning the running time, one could ask what values of η are sufficient to obtain enough meaningful pixels. Obviously this highly depends on the output's size H +1 W +1 D +1 and on the output's content itself. But we can view this question from an other perspective, by considering that we sample a constant fraction of pixels given by σ • (H +1 W +1 D +1 ) where σ ∈ [0, 1] is a sampling ratio. Because of the particular amplitudes of state ( 13), the high value pixels will be measured and known with higher probability. The points that are not sampled are being set to 0. We see that this approach is equivalent to the ∞ tomography, therefore we have

1 η 2 = σ • H +1 W +1 D +1 .

We will use this analogy in the numerical simulations (Section 6) to estimate, for a particular QCNN architecture and a particular dataset of images, which values of σ are enough to allow the neural network to learn.

D.1.6 REGULARIZATION OF THE NON LINEARITY:

In the previous steps, we see several appearances of the parameter max p,q (f (Y +1 pq )). First for the conditional rotation preprocessing, we need to know this value or an upper bound. Then for the running time, we would like to bound this parameter. Both problems can be solved by replacing the usual ReLu non linearity by a particular activation function, that we note capReLu. This function is simply a parametrized ReLu function with an upper threshold, the cap C, after which the function remain constant. The choice of C will be tuned for each particular QCNN, as a tradeoff between accuracy and speed. Otherwise, the only other requirement of the QCNN activation function would be not to allow negative values. This is already often the case for most of the classical CNN. In practice, we expect the capReLu to be as good as a usual ReLu, for convenient values of the cap C (≤ 10). We performed numerical simulations to compare the learning curve of the same CNN with several values of C. See the numerical experiments presented in Section 6 for more details. 

D.2 QRAM UPDATE:

We wish to detail the use of the QRAM between each quantum convolution layer, and present how the pooling operation can happen during this phase. General results about the QRAM is given as Theorem F.1. Implementation details can be found in the work of Kerenidis & Prakash (2017a) . In this section, we will show how to store samples from the output of the layer , to create the input of layer + 1.

D.2.1 STORING THE OUTPUT VALUES DURING THE SAMPLING:

At the beginning of layer +1, the QRAM must store A +1 , a matrix where each elements is indexed by (p , r ), and perform |p |0 → |p |A +1 p . The data is stored in the QRAM as a tree structure described by Kerenidis & Prakash (2017b) . Each row A +1 When a point is sampled from the final state of the quantum convolution, at layer , as described in Section D.1.4, we obtain a triplet corresponding to the two positions and the value of a point in the matrix f (Y +1 ). We can know where this point belong in the input of layer + 1, the tensor X +1 , by Equation ( 7), since Y is a reshaped version of X .

The position in X +1 , noted (i +1 , j +1 , d +1 ), is then matched to several positions (p , r ) in A +1 .

For each p , we write in the tree T +1 p the sampled value at leaf r and update its parent nodes, as required in the work of Kerenidis & Prakash (2017b) . Note that leaves that weren't updated will be considered as zeros, corresponding to pixels with too low values, or not selected during pooling (see next section).

Having stored pixels in this way, we can then query |p |0 → |p |A p , using the quantum circuit developed by Kerenidis & Prakash (2017b) , where we correctly have

|A +1 p = 1 A +1 p r A +1

p r |r . Note that each tree has a logarithmic depth in the number of leaves, hence the running time of writing the output of the quantum convolution layer in the QRAM gives a marginal multiplicative increase, poly-logarithmic in the number of points sampled from |f (Y +1 ) , namely

O(log(1/η 2 )).

D.2.2 QUANTUM POOLING:

As for the classical CNN, a QCNN should be able to perform pooling operations. We first detail the notations for classical pooling. At the end of layer , we wish to apply a pooling operation of size P on the output f (X +1 ). We note X +1 the tensor after the pooling operation. For a point in f (X +1 ) at position (i +1 , j +1 , d +1 ), we know to which pooling region it belongs, corresponding to a position ( ĩ +1 , j +1 , d +1 ) in X +1 :

EQUATION

)

Figure 8 : A 2×2 tensor pooling. A point in f (X +1 ) (left) is given by its position (i +1 , j +1 , d +1 ).

A point in X +1 (right) is given by its position ( ĩ +1 , j +1 , d +1 ). Different pooling regions in f (X +1 ) have separate colours, and each one corresponds to a unique point in X +1 .

We now show how any kind of pooling can be efficiently integrated to our QCNN structure. Indeed the pooling operation will occur during the QRAM update described above, at the end of a convolution layer. At this moment we will store sampled values according to the pooling rules.

In the quantum setting, the output of layer after tomography is noted X +1 . After pooling, we will describe it by X +1 , which has dimensions

H +1 P × W +1 P × D +1

. X +1 will be effectively used as input for layer + 1 and its values should be stored in the QRAM to form the trees T +1 p , related to the matrix expansion Ã +1 . However X +1 is not known before the tomography is over. Therefore we have to modify the update rule of the QRAM to implement the pooling in an online fashion, each time a sample from |f (X +1 ) is drawn. Since several sampled values of |f (X +1 ) can correspond to the same leaf Ã +1 p r (points in the same pooling region), we need an overwrite rule, that will depend on the type of pooling. In the case of Maximum Pooling, we simply update the leaf and the parent nodes if the new sampled value is higher that the one already written. In the case of Average Polling, we replace the actual value by the new averaged value.

In the end, any pooling can be included in the already existing QRAM update. In the worst case, the running time is increased by O(P/η 2 ), an overhead corresponding to the number of times we need to overwrite existing leaves, with P being a small constant in most cases.

As we will see in Section E, the final positions (p, q) that were sampled from |f (X +1 ) and selected after pooling must be stored for further use during the backpropagation phase.

D.3 RUNNING TIME:

We will now summarise the running time for one forward pass of convolution layer . With Õ we hide the polylogaryhtmic factors. We first write the running time of the classical CNN layer, which is given by O H +1 W +1 D +1 • HW D . For the QCNN, the previous steps prove Result 1 and can be implemented in time

O 1 η 2 • M √ C √ E(f (X +1 ))

. Note that, as explain in Section D.1.5, the quantum running time can also be written

O σH +1 W +1 D +1 • M √ C √ E(f (X +1 )) , with σ ∈ [0, 1]

being the fraction of sampled elements among

H +1 W +1 D +1 of them.

It is interesting to notice that the one quantum convolution layer can also include the ReLu operation and the Pooling operation in the same circuit, for no significant increase in the running time, whereas in the classical CNN each operation must be done on the whole data again.

APPENDIX E QUANTUM BACKPROGATION:

The entire QCNN is made of multiple layers. For the last layer's output, we expect only one possible outcome, or a few in the case of a classification task, which means that the dimension of the quantum output is very small. A full tomography can be performed on the last layer's output in order to calculate the outcome. The loss L is then calculated, as a measure of correctness of the predictions compared to the ground truth. As the classical CNN, our QCNN should be able to perform the optimization of its weights (elements of the kernels) to minimize the loss by an iterative method.

Theorem E.1 (Quantum Backpropagation for Quantum CNN) Given the forward pass quantum algorithm in Algorithm 1, the input matrix A and the kernel matrix F stored in the QRAM for each layer , and a loss function L, there is a quantum backpropagation algorithm that estimates, for any precision δ > 0, the gradient tensor ∂L ∂F and update each element to perform gradient descent such that ∀(s, q), ∂L ∂F s,q -∂L ∂F s,q ≤ 2δ ∂L ∂F 2 . Let ∂L ∂Y be the gradient with respect to the th layer. The running time of a single layer for quantum backpropagation is given by

EQUATION

where for a matrix V , κ(V ) is the condition number and µ(V ) is defined in Equation ( 5).

E.1 CLASSICAL BACKPROPAGATION:

After each forward pass, the outcome is compared to the true labels and define a loss. We can update our weights by gradient descent to minimize this loss, and iterate. The main idea behind the backpropagation is to compute the derivatives of the loss L, layer by layer, starting from the last one.

At layer , the derivatives needed to perform the gradient descent are ∂L ∂F and ∂L ∂Y . The first one represents the gradient of the final loss L with respect to each kernel element, a matrix of values that we will use to update the kernel weights F s,q . The second one is the gradient of L with respect to the layer itself and is only needed to calculate the gradient ∂L ∂F -1 at layer -1.

E.1.1 CONVOLUTION PRODUCT:

We first consider a classical convolution layer without non linearity or pooling. Thus the output of layer is the same tensor as the input of layer + 1, namely X +1 or equivalently Y +1 . Assuming we know ∂L ∂X +1 or equivalently ∂L ∂Y +1 , both corresponding to the derivatives of the ( + 1) th layer's input, we will show how to calculate ∂L ∂F , the matrix of derivatives with respect to the elements of the previous kernel matrix F . This is the main goal in order to optimize the kernel's weights.

The details of the following calculations can be found in the work of Wu (2017) . We will use the notation vec(X) to represents the vectorized form of any tensor X.

Recall that A is the matrix expansion of the tensor X , whereas Y is a matrix reshaping of X . By applying the chain rule

∂L ∂vec(F ) T = ∂L ∂vec(X +1 ) T ∂vec(X +1 )

∂vec(F ) T , we can obtain:

EQUATION

See calculations details in the work of Wu (2017) . Equation ( 17) shows that, to obtain the desired gradient, we can just perform a matrix-matrix multiplication between the transposed layer itself (A ) and the gradient with respect to the previous layer ( ∂L ∂Y +1 ). Equation ( 17) explains also why we will need to calculate ∂L ∂Y in order to backpropagate through layer -1. To calculate it, we use the chain rule again for

∂L ∂vec(X ) T = ∂L ∂vec(X +1 ) T ∂vec(X +1 ) ∂vec(X ) T .

Recall that a point in A , indexed by the pair (p, r), can correspond to several triplets (i , j , d ) in X . We will use the notation (p, r) ↔ (i , j , d ) to express formally this relation. One can show that ∂L ∂Y +1 (F ) T is a matrix of same shape as A , and that the chain rule leads to a simple relation to calculate ∂L ∂Y :

EQUATION

We have shown how to obtain the gradients with respect to the kernels F and to the layer itself Y (or equivalently X ).

E.1.2 NON LINEARITY:

The activation function has also an impact on the gradient. In the case of the ReLu, we should only cancel gradient for points with negative values. For points with positive value, the derivatives remain the same since the function is the identity. A formal relation can be given by

∂L ∂X +1 i +1 ,j +1 ,d +1 = ∂L ∂f (X +1 ) i +1 ,j +1 ,d +1 if X +1 i +1 ,j +1 ,d +1 ≥ 0 0 otherwise (19) E.1.3 POOLING

If we take into account the pooling operation, we must change some of the gradients. Indeed, a pixel that hasn't been selected during pooling has no impact on the final loss, thus should have a gradient equal to 0. We will focus on the case of Max Pooling (Average Pooling relies on similar idea). To state a formal relation, we will use the notations of Section D.2.2: an element in the output of the layer, the tensor f (X +1 ), is located by the triplet (i +1 , j +1 , d +1 ). The tensor after pooling is noted X +1 and its points are located by the triplet ( ĩ +1 , j +1 , d +1 ). During backpropagation, after the calculation of ∂L ∂ X +1 , some of the derivatives of f (X +1 ) should be set to zero with the following rule:

∂L ∂f (X +1 ) i +1 ,j +1 ,d +1 = ∂L ∂ X +1 ĩ +1 , j +1 , d +1 if (i +1 , j +1 , d +1 ) was selected during pooling 0 otherwise (20) E.2 QUANTUM ALGORITHM FOR BACKPROPAGATION

In this section, we want to give a quantum algorithm to perform backrpopagation on a layer , and detail the impact on the derivatives, given by the following diagram:

EQUATION

We assume that backpropagation has been done on layer +1. This means in particular that ∂L ∂X +1 is stored in QRAM. However, as shown on Diagram (21), ∂L ∂X +1 corresponds formally to ∂L ∂ X +1 , and not ∂L ∂X +1 . Therefore, we will have to modify the values stored in QRAM to take into account non linearity, tomography and pooling. We will first consider how to implement ∂L ∂X and ∂L ∂F through backpropagation, considering only convolution product, as if ∂L ∂X +1 and ∂L ∂X +1 where the same. Then we will detail how to simply modify ∂L ∂X +1 a priori, by setting some of its values to 0.

E.2.1 QUANTUM CONVOLUTION PRODUCT:

In this section we consider only the quantum convolution product without non linearity, tomography nor pooling, hence writing its output directly as X +1 . Regarding derivatives, the quantum convolution product is equivalent to the classical one. Gradient relations ( 17) and ( 18) remain the same. Note that the -approximation from Section D.1.2 doesn't participate in gradient considerations.

The gradient relations being the same, we still have to specify the quantum algorithm that implements the backpropagation and outputs classical description of ∂L ∂X and ∂L ∂F . We have seen that the two main calculations ( 17) and ( 18) are in fact matrix-matrix multiplications both involving ∂L ∂Y +1 , the reshaped form of ∂L ∂X +1 . For each, the classical running time is O(H +1 W +1 D +1 HW D ). We know from Theorem F.7 and Theorem G.1 a quantum algorithm to perform efficiently a matrix-vector multiplication and return a classical state with ∞ norm guarantees. For a matrix V and a vector b, both accessible from the QRAM, the running time to perform this operation is O µ(V )κ(V ) log 1/δ δ 2

, where κ(V ) is the condition number of the matrix and µ(V ) is a matrix parameter defined in Equation ( 5). Precision parameter δ > 0 is the error committed in the approximation for both Theorems F.7 and G.1.

We can therefore apply theses theorems to perform matrix-matrix multiplications, by simply decomposing them in several matrix-vector multiplications. For instance, in Equation ( 17), the matrix could be (A ) T and the different vectors would be each column of ∂L ∂Y +1 . The global running time to perform quantumly Equation ( 17) is obtained by replacing µ(V ) by µ( ∂L ∂Y +1 ) + µ(A ) and κ(V ) by κ((A ) T • ∂L ∂Y +1 ). Likewise, for Equation ( 18), we have µ( ∂L ∂Y +1 )+µ(F ) and κ( ∂L ∂Y +1 •(F ) T ). Note that the dimension of the matrix doesn't appear in the running time since we tolerate a ∞ norm guarantee for the error, instead of a 2 guarantee (see Section G for details). The reason why ∞ tomography is the right approximation here is because the result of these linear algebra operations are rows of the gradient matrices, that are not vectors in an euclidean space, but a series of numbers for which we want to be δ-close to the exact values. See next section for more details.

It is a open question to see if one can apply the same sub-sampling technique as in the forward pass (Section D.1) and sample only the highest derivatives of ∂L ∂X , to reduce the computation cost while maintaining a good optimization. We then have to understand which elements of ∂L ∂X +1 must be set to zero to take into account the effects the non linearity, tomography and pooling.

E.2.2 QUANTUM NON LINEARITY AND TOMOGRAPHY:

To include the impact of the non linearity, one could apply the same rule as in ( 19), and simply replace ReLu by capReLu. After the non linearity, we obtain f (X +1 ), and the gradient relation would be given by

∂L ∂X +1 i +1 ,j +1 ,d +1 =    ∂L ∂f (X +1 ) i +1 ,j +1 ,d +1 if 0 ≤ X +1 i +1 ,j +1 ,d +1 ≤ C 0 otherwise (22)

If an element of X +1 was negative or bigger than the cap C, its derivative should be zero during the backpropagation. However, this operation was performed in quantum superposition. In the quantum algorithm, one cannot record at which positions (i +1 , j +1 , d +1 ) the activation function was selective or not. The gradient relation ( 22) cannot be implemented a posteriori. We provide a partial solution to this problem, using the fact that quantum tomography must also be taken into account for some derivatives. Indeed, only the points (i +1 , j +1 , d +1 ) that have been sampled should have an impact on the gradient of the loss. Therefore we replace the previous relation by

∂L ∂X +1 i +1 ,j +1 ,d +1 = ∂L ∂X +1 i +1 ,j +1 ,d +1 if (i +1 , j +1 , d +1 ) was sampled 0 otherwise (23)

Nonetheless, we can argue that this approximation will be tolerable. In the first case where X +1 i +1 ,j +1 ,d +1 < 0, the derivatives can not be set to zero as they should. But in practice, their values will be zero after the activation function and such points would not have a chance to be sampled. In conclusion their derivatives would be zero as required. In the other case, where X +1 i +1 ,j +1 ,d +1 > C, the derivatives can not be set to zero as well but the points have a high probability of being sampled. Therefore their derivative will remain unchanged, as if we were using a ReLu instead of a capReLu. However in cases where the cap C is high enough, this shouldn't be a source of disadvantage in practice.

E.2.3 QUANTUM POOLING:

From relation (23), we can take into account the impact of quantum pooling (see Section D.2.2) on the derivatives. This case is easier since one can record the selected positions during the QRAM update. Therefore, applying the backpropagation is similar to the classical setting with Equation (20).

∂L ∂X:

+1 i +1 ,j +1 ,d +1 = ∂L ∂ X +1 ĩ +1 , j +1 , d +1 if (i +1 , j +1 , d +1 ) was selected during pooling 0 otherwise (24)

Note that we know ∂L ∂ X +1 as it is equal to ∂L ∂X +1 , the gradient with respect to the input of layer +1, known by assumption and stored in the QRAM.

E.3 CONCLUSION AND RUNNING TIME:

In conclusion, given ∂L ∂Y +1 in the QRAM, the quantum backpropagation first consists in applying the relations (24) followed by ( 23). The effective gradient now take into account non linearity, tomography and pooling that occurred during layer . We can know use apply the quantum algorithm for matrix-matrix multiplication that implements relations (18) and (17).

Note that the steps in Algorithm 2 could also be reversed: during backpropagation of layer + 1, when storing values for each elements of ∂L ∂Y +1 in the QRAM, one can already take into account ( 24) and ( 23) of layer . In this case we directly store ∂L ∂X +1 , at no supplementary cost. Therefore, the running time of the quantum backpropagation for one layer , given as Algorithm 2, corresponds to the sum of the running times of the circuits for implementing relations ( 17) and ( 18).

We finally obtain

O µ(A ) + µ( ∂L ∂Y +1 ) κ((A ) T • ∂L ∂Y +1 ) + µ( ∂L ∂Y +1 ) + µ(F ) κ( ∂L ∂Y +1 • (F ) T ) log 1/δ δ 2

, which can be rewritten as

EQUATION

Besides storing ∂L ∂X , the main output is a classical description of ∂L ∂F , necessary to perform gradient descent of the parameters of F . In the Appendix (Section E.4), which details the impact of the quantum backpropagation compared to the classical case, which can be reduced to a simple noise addition during the gradient descent.

E.4 QUANTUM GRADIENT DESCENT AND CLASSICAL EQUIVALENCE:

In this part we will see the impact of the quantum backpropagation compared to the classical case, which can be reduced to a simple noise addition during the gradient descent. Recall that gradient descent, in our case, would consist in applying the following update rule F ← F -λ ∂L ∂F with the learning rate λ.

Let's note x = ∂L ∂F and its elements x s,q = ∂L ∂F s,q

. From the first result of Theorem F.7 with error δ < 0, and the tomography procedure from Theorem G.1, with same error δ, we can obtain a classical description of x x 2 with ∞ norm guarantee, such that:

x x 2 - x x 2 ∞ ≤ δ in time O( κ(V )µ(V ) log(δ) δ 2

), where we note V is the matrix stored in the QRAM that allows to obtain x, as explained in Section E.2. The ∞ norm tomography is used so that the error δ is at most the same for each component

∀(s, q), x s,q x 2 - x s,q x 2 ≤ δ

From the second result of the Theorem F.7 we can also obtain an estimate x 2 of the norm, for the same error δ, such that

| x 2 -x 2 | ≤ δ x 2 in time O( κ(V )µ(V )

δ log(δ)) (which does not affect the overall asymptotic running time). Using both results we can obtain an unnormalized state close to x such that, by the triangular inequality

x -x ∞ = x x 2 x 2 - x x 2 x 2 ∞ ≤ x x 2 x 2 - x x 2 x 2 ∞ + x x 2 x 2 - x x 2 x 2 ∞ ≤ 1 • | x 2 -x 2 | + x 2 • x x 2 - x x 2 ∞ ≤ δ x 2 + x 2 δ ≤ 2δ x 2 in time O( κ(V )µ(V ) log(δ) δ 2

). In conclusion, with ∞ norm guarantee, having also access to the norm of the result is costless.

Finally, the noisy gradient descent update rule, expressed as F s,q ← F s,q -λ ∂L ∂F s,q can written in the worst case with

EQUATION

To summarize, using the quantum linear algebra from Theroem F.7 with ∞ norm tomography from Theroem G.1, both with error δ, along with norm estimation with relative error δ too, we can obtain classically the unnormalized values ∂L ∂F such that ∂L

EQUATION

Therefore the gradient descent update rule in the quantum case becomes F s,q ← F s,q -λ ∂L ∂F s,q , which in the worst case becomes

EQUATION

This proves the Theorem E.1. This update rule can be simulated by the addition of a random relative noise given as a gaussian centered on 0, with standard deviation equal to δ. This is how we will simulate quantum backpropagation in the Numerical Simulations.

Compared to the classical update rule, this corresponds to the addition of noise during the optimization step. This noise decreases as ∂L ∂F 2 , which is expected to happen while converging. Recall that the gradient descent is already a stochastic process. Therefore, we expect that such noise, with acceptable values of δ, will not disturb the convergence of the gradient, as the following numerical simulations tend to confirm.

APPENDIX F PRELIMINARIES IN QUANTUM INFORMATION:

We introduce a basic and broad-audience quantum information background necessary for this work. For a more detailed introduction we recommend Nielsen & Chuang (2002a).

F.1 QUANTUM INFORMATION:

Quantum Bits and Quantum Registers: The bit is the most basic unit of classical information. It can be either in state 0 or 1. Similarly a quantum bit or qubit, is a quantum system that can be is state |0 , |1 (the braket notation |• is a reminder that the bit considered is a quantum system) or in superposition of both states α |0 + β |1 with coefficients α, β ∈ C such that |α| 2 + |β| 2 = 1. The amplitudes α and β are linked to the probabilities of observing either 0 or 1 when measuring the qubit, since P (0) = |α| 2 and P (1) = |β| 2 . Before the measurement, any superposition is possible, which gives quantum information special abilities in terms of computation. With n qubits, the 2 n possible binary combinations can exist simultaneously, each with a specific amplitude. For instance we can consider an uniform distribution In its most general formulation, a quantum state with n qubits can be seen as vector in a complex Hilbert space of dimension 2 n . This vector must be normalized under 2 -norm, to guarantee that the squared amplitudes sum to 1.

Quantum Computation: To process qubits and therefore quantum registers, we use quantum gates. These gates are unitary operators in the Hilbert space as they should map unit-norm vectors to unit-norm vectors. Formally, we can see a quantum gate acting on n qubits as a matrix U ∈ C 2 n such that U U † = U † U = I, where U † is the conjugate transpose of U . Some basic single qubit gates includes the NOT gate 0 1 1 0 that inverts |0 and |1 , or the Hadamard gate 1

√ 2 1 1 1 -1 that maps |0 → 1 √ 2 (|0 + |1 ) and |1 → 1 √ 2 (|0 -|1

), creating the quantum superposition. Finally, multiple qubits gates exist, such as the Controlled-NOT that applies a NOT gate on a target qubit conditioned on the state of a control qubit.

The main advantage of quantum gates is their ability to be applied to a superposition of inputs. Indeed, given a gate U such that U |x → |f (x) , we can apply it to all possible combinations of x at once U ( 1

C x |x ) → 1 C x |f (x) .

We now state some primitive quantum circuits, which we will use in our algorithm: 2018) ) Let M ∈ R d×d and x ∈ R d . Let δ 1 , δ 2 > 0. If M is stored in appropriate QRAM data structures and the time to prepare |x is T x , then there exist quantum algorithms that with probability at least 1 -1/poly(d) return

1. A state |z such that |z -|M x 2 ≤ δ 1 in time O((κ(M )µ(M ) + T x κ(M )) log(1/δ 1 )).

Note that this also implies

|z -|M x ∞ ≤ δ 1 2. Norm estimate z ∈ (1 ± δ 2 ) M x 2 , with relative error δ 2 , in time O(T x κ(M )µ(M ) δ2

log(1/δ 1 )).

The linear algebra procedures above can also be applied to any rectangular matrix V ∈ R n×d by considering instead the symmetric matrix V = 0 V V T 0 .

APPENDIX G ALGORITHM AND PROOF FOR ∞ NORM TOMOGRAPHY Finally, we present a logarithmic time algorithm for vector state tomography that will be used to recover classical information from the quantum states with ∞ norm guarantee. Given a unitary U that produces a quantum state |x = 1

x 2 d-1 j=0 x j |j , by calling O(log d/δ 2 ) times U , the tomography algorithm is able to reconstruct a vector X that approximates |x with ∞ norm guarantee, such that | X -|x ∞ ≤ δ, or equivalently that ∀i ∈ [d], |x i -X i | ≤ δ. Such a tomography is of interest when the components x i of a quantum state are not the coordinates of an meaningful vector in some linear space, but just a series of values, such that we don't want an overall guarantee on the vector (which is the case with usual 2 tomography) but a similar error guarantee for each component in the estimation. The proof of this theorem is similar to the proof of the 2 -norm tomography by Kerenidis & Prakash (2018) . However the ∞ norm tomography introduced in this paper depends only logarithmically and not linearly in the dimension d. Note that in our case, T (U ) will be logarithmic in the dimension.

Theorem G.2 [ 2 Vector state tomography Kerenidis & Prakash (2018) ] Given access to unitary U such that U |0 = |x and its controlled version in time T (U ), there is an algorithm that allows to output a classical vector X ∈ R d with 2 -norm guarantee X -x 2 ≤ δ for any δ > 0, in time

O(T (U ) × d log(d) δ 2

).

In the following we consider a quantum state |x = i∈ [d] x i |i , with x ∈ R d and x 2 = 1.

The following version of the Chernoff Bound will be used for analysis of algorithm 3.

Theorem G.3 (Chernoff Bound) Let X j , for j ∈ [N ], be independent random variables such that X j ∈ [0, 1] and let X = j∈[N ] X j . We have the three following inqualities: Case 1 : i ∈ S. We will show that if i ∈ S, we correctly have σ(i) = sgn(x i ) with high probability. Therefore we will need to bound |x i -σ(i)

1. For 0 < β < 1, P[X < (1 -β)E[X]] ≤ e -β 2 E[X]/2

√ p i | = ||x i | - √ p i |.

We suppose that x i > 0. The value of σ(i) correctly determines sgn(x i ) if the number of times we have measured (0, i) at Step 4. is more than half of the outcomes, i.e. n(0, i) > 1 2 E[n(0, i)]. If x i < 0, the same arguments holds for n(1, i). We consider the random variable that represents the outcome of a measurement on state |φ . The Chernoff Bound, part 1 with β = 1/2 gives P[n(0, i) ≤ 1 2 E[n(0, i)]] ≤ e -E[n(0,i)]/8 (29)

From the definition of |φ we have E[n(0, i)] = N 4 (x i + √ p i ) 2 . We will lower bound this value with the following argument.

For the k th measurement of |x , with k ∈ [N ], let X k be a random variable such that X k = 1 if the outcome is i, and 0 otherwise. We define X = k∈[N ] X k . Note that X = n i = N p i and E[X] = N x 2 i . We can apply the Chernoff Bound, part 3 on X for β = 1/2 to obtain, We conclude that for i ∈ S, if n(0, i) > 0.41N p i , the sign of x i is determined correctly by σ(i) with high probability 1 - 

P[|X -E[X]| ≥ E[X]/2] ≤ e -E[

P |x i | - √ p i ≥ β x 2 i |x i | + √ p i ≤ e β 2 N x 2 i /3

Since √ p i > 0, we have β We conclude that, if i ∈ S, we have |x i -Xi | ≤ δ with high probability.

Since |S| ≤ d, the probability for this result to be true for all i ∈ S is 1 -1 d 0.83 . This can be proved by using the Union Bound on the correctness of σ(i).

Case 2 : i / ∈ S.

If i /:

∈ S, we need to separate again in two cases. When the estimated sign is wrong, i.e. σ(i) = -sgn(x i ), we have to bound |x i -σ(i)

√ p i | = ||x i | + √ p i |.

On the contrary, if it is correct, i.e.

σ(i) = sgn(x i ), we have to bound |x i -σ(i)

√ p i | = ||x i | - √ p i | ≤ ||x i | + √ p i |.

Therefore only one bound is necessary.

We use Chernoff Bound (2.) on the random variable X with β > 0 to obtain

P[p i > (1 + β)x 2 i ] ≤ e β 2 2+β N x 2 i

We chose β = δ 2 /x 2 i and obtain P[p i > x 2 i + δ 2 ] ≤ e δ 4 3δ 2 N = 1/d 12 . Therefore, if i / ∈ S, with very high probability 1 -1 d 12 we have p i ≤ x 2 i + δ 2 ≤ 2δ 2 . We can conclude and bound the error:

|x i -Xi | ≤ ||x i | + √ p i | ≤ δ + √ 2δ = (1 + √ 2)δ

Since |S| ≤ d, the probability for this result to be true for all i / ∈ S is 1 - 

:



Figures and Tables:

Step 2: Quantum Sampling Use Conditional Rotation and Amplitude Amplification to encode the values α pq := α pq |p |q |f (Y +1 pq ) |g pq . Perform ∞ tomography from Theorem G.1 with precision η, and obtain classically all positions and values (p, q, f (Y +1 pq )) such that, with high probability, values above η are known exactly, while others are set to 0. 3: Step 3: QRAM Update and Pooling

Theorem B.5 [Quantum Distance Estimation Wiebe et al. (2014b); Kerenidis et al. (2019)]

Figure 4: Convolution of a 3-tensor input (Left) by one 3-tensor kernel (Center). The ouput (Right) is a matrix for which each entry is a inner product between the kernel and the corresponding overlapping region of the input.

Figure 5: Convolutions of the 3-tensor input X (Left) by one 4-tensor kernel K (Center). Each channel of the output X +1 (Right) corresponds to the output matrix of the convolution with one of the 3-tensor kernel.

to denote a pixel's value, with n ∈ [H +1 W +1 D +1 ]. Recall that Y +1 and X +1 are reshaped version of the same object.

Figure 7: Activation functions: ReLu (Left) and capReLu (Right) with a cap C at 5.

is stored in such a tree T +1 p . Each leafA +1p r correspond to a value sampled from the previous quantum state |f (Y +1 ) , output of the layer . The question is to know where to store a sample from |f (Y +1 ) in the tree T +1 p .

|i where |i represents the i th binary combination (e.g. |01 • • • 1001 ). Multiple qubits together are often called a quantum register.

two integers i and j, we can check their equality with the mapping |i |j |0 → |i |j |[i = j] . For two real value numbers a > 0 and δ > 0, we can compare them using |a |δ |0 → |a |δ |[a ≤ δ] . Finally, for a real value numbers a > 0, we can obtain its square |a |0 → |a |a 2 . Theorem F.7 (Quantum Matrix Operations, Chakraborty et al. (

∞ Vector state tomography) Given access to unitary U such that U |0 = |x and its controlled version in time T (U ), there is a tomography algorithm with time complexity O(T (U ) log d δ 2 ) that produces unit vector X ∈ R d such that X -x ∞ ≤ δ with probability at least (1 -1/poly(d)).

For β > 0, P[X > (1 + β)E[X]] ≤ e -β 2 2+β E[X] 3. For 0 < β < 1, P[|X -E[X]| ≥ βE[X]] ≤ e -β 2 E[X]/3 , by composing 1. and 2. Algorithm 3 ∞ norm tomography Require: Error δ > 0, access to unitary U : |0 → |x = i∈[d] x i |i , the controlled version of U , QRAM access. Ensure: Classical vector X ∈ R d , such that X = 1 and X -x ∞ < δ. 1: Measure N = 36 ln d δ 2copies of |x in the standard basis and count n i , the number of times the outcome i is observed. Store√ p i = n i /N in QRAM data structure. i + √ p i ) |0, i + (x i -√ p i ) |1, i )4: Measure both registers of each copy in the standard basis, and count n(0, i) the number of time the outcome (0, i) is observed.5: Set σ(i) = +1 if n(0, i) > 0.4N p i and σ(i) = -1 otherwise. 6: Output the unit vector X such that ∀i ∈ [N ], X i = σ i √ p i Theorem G.4 Algorithm 3 produces an estimate X ∈ R d such that X -x δ) is equivalent to show that for all i ∈ [d], we have |x i -X i | = |x i -σ(i) √ p i | ≤ O(δ).Let S be the set of indices defined by S = {i ∈ [d]; |x i | > δ}. We will separate the proof for the two cases where i ∈ S and i / ∈ S.

i |, therefore P |x i | -√ p i ≥ β|x i | ≤ e β 2 N x 2 i /3 . Finally, by chosing β = δ/|x i | < 1 we have P |x i | -√ p i ≥ δ ≤ e 36 ln d/3 = 1/d 12

(|p |q |0 |A p + |p |q |1 |F q ). By applying a Hadamard gate on the third register we obtain the following state, 1 2 |p |q |0 (|A p +|F q )+ |1 (|A p -|F q ) . The probability of measuring 0 on the third register is given by P pq = Thus we can rewrite the previous state as |p |q P pq |0, y pq 1 -P pq |1, y pq , where |y pq and |y pq are some garbage states. We can perform the previous circuit in superposition. Since A has H +1 W +1 rows, and F has D +1 columns, we obtain the state: |u =

Table content: <table><tr><td/><td colspan="4">1 √ 2 (|p |q |0 |0 + |p |q |1 |0 ) → 1 √ 2 1+ Ap|Fq 2</td><td>.</td></tr><tr><td>√</td><td>1 H +1 W +1 D +1</td><td>p</td><td>q |p |q</td><td>P</td></tr></table>

where |g pq is another garbage state. The overall probability of measuring |0 on the last register is P (0) =

Table content: <table><tr><td>1 HW D</td><td>pq |α pq | 2 . The</td></tr><tr><td>number of queries required to amplify the state |0 is O( 1 √</td><td/></tr></table>

This proves that the event |x 2 i -p i | ≤ x 2 i /2 occurs with probability at least 1 -1d 3 if i ∈ S. This previous inequality is equivalent to 2p i /3 ≤ |x i | ≤ √ 2p i .Thus, with high probability we haveE[n(0, i)] = N 4 (x i + √ p i ) 2 ≥ 0.82N p i , since 2p i /3 ≤ |x i |. Moreover, since |p i | ≤ x 2 i /2, E[n(0, i)] ≥ 0.82N x 2 i /2 ≥ 14.7 ln d. Therefore, equation equation 29 becomes P[n(0, i) ≤ 0.41N p i ] ≤ e -1.83 ln d = 1/d 1.83

Table content: <table><tr><td>X]/12</td><td>(30)</td></tr><tr><td>P[|x 2</td><td/></tr></table>

1 d 1.83 , as indicated in Step 5. We finally show |xi -σ(i) √ p i | = ||x i | -√ p i | is bounded.Again by the Chernoff Bound (3.) we have, for 0 < β < 1:P[|x 2 i -p i | ≥ βx 2 i ] ≤ e β 2 N x 2 By the identity |x 2 i -p i | = (|x i | -√ p i )(|x i | + √ p i ) we have

Table content: <table><tr><td>i /3</td></tr></table>

1 d 11 . This follows from applying the Union Bound on the eventp i > x 2 i + δ 2 .the accuracy is just over 96%, lower than the best possible accuracy with larger CNN. QCNN trained with quantum backpropagation on MNIST dataset. With C = 10 fixed.

Table content: <table><tr><td/><td/><td colspan="4">QCNN Test -Classification</td><td/><td/></tr><tr><td>σ</td><td>δ</td><td>0.01</td><td>0.01</td><td>0.1</td><td>0.01</td><td>0.1</td><td>0.1</td></tr><tr><td>0.1</td><td colspan="7">Loss Accuracy 82.8% 74.8% 11.5% 11.7% 0.519 0.773 2.30 2.30</td></tr><tr><td>0.2</td><td colspan="7">Loss Accuracy 89.5% 89.0% 86.2% 54.1% 0.334 0.348 0.439 1.367</td></tr><tr><td>0.3</td><td colspan="7">Loss Accuracy 93.4% 90.3% 87.9% 76.8% 0.213 0.314 0.381 0.762</td></tr><tr><td>0.4</td><td colspan="7">Loss Accuracy 94.7% 93.3% 91.8% 34.9% 0.177 0.215 0.263 1.798</td></tr><tr><td>0.5</td><td colspan="7">Loss Accuracy 95.4% 93.5% 89.2% 52.8% 0.142 0.211 0.337 1.457</td></tr><tr><td/><td/><td colspan="4">QCNN Test -Classification</td><td/><td/></tr><tr><td/><td>σ</td><td/><td/><td>0.01</td><td>0.1</td><td/><td/></tr><tr><td/><td>0.1</td><td colspan="5">Loss Accuracy 86.1% 78.6% 1.07 1.33</td><td/></tr><tr><td/><td>0.2</td><td colspan="5">Loss Accuracy 92.8% 86.5% 0.552 0.840</td><td/></tr><tr><td/><td>0.3</td><td colspan="5">Loss Accuracy 94,3% 85.8% 0.391 0.706</td><td/></tr><tr><td/><td>0.4</td><td colspan="5">Loss Accuracy 94.4% 84.0% 0.327 0.670</td><td/></tr><tr><td/><td>0.5</td><td colspan="5">Loss Accuracy 95.9% 93.5% 0.163 0.292</td><td/></tr></table>

QCNN created from a classical CNN trained on MNIST dataset. With δ = 0.01 and C = 10 fixed.

Table content: <table/>

Peer review (iclr2020):
The authors present a quantum algorithm for approximating the forward pass and gradient computation of a classical convolutional neural network layer with pooling and a bounded rectifier activation. This algorithm has complexity bounds that would open up (for instance) the possibility of exponentially large filter banks, and the authors show through a simple, classical simulation approach that the resulting network is also likely to be trainable. Feedback: A few typos/formatting issues: - The title accidentally includes *Conference Submissions* - The in-text citation format frequently has the parentheses in the wrong place; this is surprisingly distracting! Preliminaries: - Maybe explain what the ith vector in the standard basis is in terms of |0> and |1>? I assume the answer is along the lines of |000>, |001>, |010>, etc.? Main results: - The sentence *a speedup compared to the classical CNN for both the forward pass and for training using backpropagation in certain cases* is ambiguous; does *in certain cases* qualify only training speed or also forward pass speed? - There*s a clear separation of background (which is concise and well explained) and contributions, but maybe it would be worth connecting the introduced algorithm more closely to existing work in non-convolutional quantum neural networks? - Can you briefly justify (or cite) the claim that *most of the non linear functions in the machine learning literature can be implemented using small sized boolean circuits*? - I*m a little confused about the discussion of quantum importance sampling on page 4. Could you give some intuition for the relationship between eta and the fraction of output values that are on average flushed to zero (is this 1 minus sigma?), and perhaps connect this to the literature about activation pruning and sparse NNs? - Maybe define what you mean by *tomography* for ML folks without the quantum background? - I*m convinced by the simulations, even though I shouldn*t really be convinced by anything on MNIST... It just seems like the perturbations you*re applying are all things that modern neural networks take in stride. - The discussion of using a sigma-based classical sampling rather than the eta-based quantum importance sampling mentions a *Section C.1.15* which does not exist (I think you mean the end of Section C.1.5). - Re: *We will use this analogy in the numerical simulations (Section 6) to estimate, for a particular QCNN architecture and a particular dataset of images, which values of ? are enough to allow the neural network to learn.* My understanding is that you*re getting empirical estimates of which values of sigma are enough; it would be valuable to convert those to estimates of which values of eta would be enough (given quantum networks of the size used in the classical simulation experiment, or given larger networks). - The sampling procedure based on sigma might be inefficient in your PyTorch implementation, but it*s certainly something that GPUs are fairly well suited to computing. There might be other PyTorch operators that would help here (perhaps Bernoulli sampling?) or if nothing else you could write a small custom CUDA kernel.