Quantum optical classifier with superexponential speedup
Simone Roncallo Angela Rosy Morgillo Chiara Macchiavello Lorenzo Maccone Seth Lloyd




	[Simone Roncallo: ]simone.roncallo01@ateneopv.it
	Dipartimento di Fisica, Università degli Studi di Pavia, Via Agostino Bassi 6, I-27100, Pavia, Italy
	INFN Sezione di Pavia, Via Agostino Bassi 6, I-27100, Pavia, Italy
	


	[Angela Rosy Morgillo: ]angelarosy.morgillo01@ateneopv.it
	Dipartimento di Fisica, Università degli Studi di Pavia, Via Agostino Bassi 6, I-27100, Pavia, Italy
	INFN Sezione di Pavia, Via Agostino Bassi 6, I-27100, Pavia, Italy
	


	[Chiara Macchiavello: ]chiara.macchiavello@unipv.it
	Dipartimento di Fisica, Università degli Studi di Pavia, Via Agostino Bassi 6, I-27100, Pavia, Italy
	INFN Sezione di Pavia, Via Agostino Bassi 6, I-27100, Pavia, Italy
	


	[Lorenzo Maccone: ]lorenzo.maccone@unipv.it
	Dipartimento di Fisica, Università degli Studi di Pavia, Via Agostino Bassi 6, I-27100, Pavia, Italy
	INFN Sezione di Pavia, Via Agostino Bassi 6, I-27100, Pavia, Italy
	


	[Seth Lloyd: ]slloyd@mit.edu
	Massachusetts Institute of Technology, Cambridge, MA 02139, USA
	

	We present a quantum optical pattern recognition method for binary classification tasks. Without direct image reconstruction, it classifies an object in terms of the rate of two-photon coincidences at the output of a Hong-Ou-Mandel interferometer, where both the input and the classifier parameters are encoded into single-photon states. Our method exhibits the same behaviour of a classical neuron of unit depth. Once trained, it shows a constant 𝒪(1) complexity in the number of computational operations and photons required by a single classification. This is a superexponential advantage over a classical neuron (that is at least linear in the image resolution). We provide simulations and analytical comparisons with analogous neural network architectures.


Quantum optical classifier with superexponential speedup
    Seth Lloyd 0000-0003-0353-4529
    April 29, 2024
========================================================





§ INTRODUCTION

Image classification has been significantly affected by the introduction of deep learning algorithms, providing several architectures that can learn and extract image features <cit.>. The large number of parameters involved is motivating a consistent effort in reducing the cost of these methods, e.g. by leveraging all-optical implementations that bypass hardware usage <cit.>, or quantum mechanical effects that can provide a significant speedup in these computations <cit.>. Quantum optical neural networks harness the best of both worlds, i.e. deep learning capabilities from quantum optics <cit.>.

In this paper, we introduce a quantum optical setup to classify objects without reconstructing their images. Our approach relies on the Hong-Ou-Mandel effect, for which the probability that two photons exit a beam splitter in different modes, depends on their distinguishability <cit.>. In our implementation, an input object is targeted by a single-photon source, and eventually followed by an arbitrary lens system. The single-photon state interferes with another one, which encodes a set of trainable parameters, e.g. through a spatial light modulator. Classification occurs by measuring the rate of two-photon coincidences at the Hong-Ou-Mandel output (see <ref>). The Hong-Ou-Mandel effect has been successfully applied to quantum kernel evaluation <cit.>, which can compute distances between pairs of data points in the feature space. In this case, each point is sent to one branch of the interferometer, encoded in the temporal modes of a single-photon state. In our method, the interferometer has only one independent branch, which takes the spatial modes of a single-photon state reflected off the target object. The other branch remains fixed after training, and contains the layer of parameters. After the measurement, the response function of our apparatus mathematically resembles that of a classical neuron. For this reason, we refer to our setup as quantum optical neuron. By analytically comparing the resource cost of the classical and quantum neurons, we show that our method requires constant 𝒪(1) computational operations and injected photons, whereas the classical methods are at least linear in the image resolution: a superexponential advantage.




§ METHOD

In this section, we discuss the apparatus of <ref>, without explicitly modelling the probe. Two single-photon states are fed into the left and top branches of a 50:50 beam splitter, acting as input and processing layers, respectively. In the left branch, the single-photon source reflects off the object, and reaches the beam splitter after a linear optical system. In the top branch, we consider a generic single-photon state, which depends on a set of trainable real parameters. We count the two-photon coincidences at the beam splitter output. We show how to interpret the Hong-Ou-Mandel response as the one produced by a single-layer neural network-like operation on the object image.

We call input and probe modes, i.e. a and b, those fed into the left and top branches of the interferometer. In the input branch, a single photon with spectrum ϕ is generated at the longitudinal origin z=0, followed by an object with two-dimensional shape Ø. An imaging system with transfer function ℒ_d, e.g a pinhole or a linear optical apparatus, is placed after the object. Here, z_o and z_i are the longitudinal positions of the object and the image plane, respectively, and d = z_i - z_o their displacement.

The output of the imaging optics reads (see <ref>)

    |Ψ_⟩ = ∫^2k  _ω(k | Ø) a^†_ω(k) |0⟩ ,

with _ω(·|Ø) = [(ϕ̂_ωℌ̂_z_o)*Ø̂] ℒ̂_d the total transfer function from the single-photon source to the image plane, and a^†_ω(k) the creation operator of a photon in the input mode, acting on the vacuum state |0⟩. The hat operator denotes the two-dimensional Fourier transform on the transverse coordinates plane, * the convolution operation, ℌ_z_o the transfer function from the source to the object plane, k = (k_x, k_y) the transverse momentum, and ω the frequency conjugated to the temporal degree of freedom of the electromagnetic potential.

In the probe branch, a generic quantum state is prepared, eventually followed by a linear optical system. At the beam splitter plane, the probe state reads

    |Ψ_⟩ = ∫^2k  _ω(k|λ) b_ω^†(k) |0⟩ ,

with λ = {λ_i_1 … i_n} a collection of (trainable) parameters,  the spatial spectrum of the probe, and b^†_ω(k) the creation operator of a photon in the probe mode.

A photodetector is placed at the output of each branch. After feeding both states into a 50:50 beam splitter, the rate of two-photon coincidences reads

    p(1_a ∩ 1_b|λ, Ø) = 1/2[α_λ(Ø) - _λ(Ø) ]  ,
        
       with  α_λ(Ø)    = || _ω(·|Ø) ||^2 || _ω(·|λ) ||^2   , 
    _λ(Ø)    =  | ⟨_ω(·|Ø),_ω (·|λ) ⟩|^2  ,

where || · || and ⟨· , ·⟩ denote the L^2-norm and inner product, respectively. Here, α_λ(Ø) depends on the normalization of the input and probe states, which can be α_λ < 1 in the presence of optical losses. Whenever the two spectra are indistinguishable, i.e. when  perfectly matches , coincidences are not observed. On the other hand, the more distinguishable the input and the probe states are, the smaller ⟨(·|Ø),(·|λ) ⟩ becomes and the rate of coincidences increases. See <ref> for a derivation, and <ref>, for a similar result in the Fourier domain.

At the image plane I, with transverse coordinates r = (x,y), we have

    _λ(Ø) = | ∫_I  r  _ω(r|Ø)_ω^*(r|λ) |^2  .

This integral measures the point-wise overlap between the input image and the probe. We interpret it as the prediction of our classification model, where _λ∈ [0,1] represents the probability that  belongs to the class of . In particular, _λ→ 0 (_λ→ 1) when the class of  is orthogonal to (is the same of) . In the next section, we show how to encode a generic class in , by means of the optimization of the set of parameters λ.

The output measurement introduces a non-linear operation after the beam splitter, represented by the squared absolute value in the left-hand side of <ref>. We increase the predictability of our model, by enhancing this non-linearity through the following post-processing operations. Consider the sigmoid (logistic) function

    σ(x) := 1/1+e^-β x + γ ,

where β, γ are hyperparameters, i.e. constants with respect to the training process. We introduce an additional trainable parameter b ∈ℝ, called bias, which, combined with f_λ and σ, yields

    _bλ(Ø) = σ(_λ(Ø) + b)  ,

which determines the label predicted by the Hong-Ou-Mandel apparatus. These modifications can improve the performance of the neuron. The sigmoid increases the non-linearity introduced by the squared absolute value, and so the predictability of the model. In addition, the bias is introduced on heuristic motivations: it compensates the constraint given by the normalization in <ref>, while enhancing the robustness of our protocol against optical losses (which may affect the above-mentioned normalizability, yielding α_λ < 1).

We now discuss the training stage. Consider a training set, i.e. an ensemble of objects {Ø_j} with target labels {y_j ∈{0,1}}. We separately feed each object into the input branch of the interferometer. Predicted and target classes are compared in terms of their binary cross-entropy, which is used as loss function of a gradient descent optimizer. The optimizer updates λ through the derivative of the loss function, whose only model-dependent contribution is

    ∂_λ = 2 [ ⟨_ω, _ω⟩⟨_ω, ∂_λ_ω⟩^* ]  .
	
Ideally, the training is complete after finding a set of parameters  that minimizes the loss. Notice that our model is resilient against the issue of gradient explosion <cit.>, since it depends on physical data and functions only (see <ref> for a discussion).

There is a formal relationship between the post-processed output of the Hong-Ou-Mandel interferometer of <ref> and that of a classical neuron. Consider f_λ(Ø) discretized and vectorized in a mesh of N cells, either in the spatial or in the Fourier domain. Then, <ref> corresponds to the composition of a real-valued neuron, with N trainable weights, square absolute value activation function and no bias, and a second neuron, with a scalar unit weight, sigmoid activation function and a trainable bias. Namely

    G_bw(x) = σ(|w · x|^2 + b )  ,

where x ∈ℂ^N is the input, while w ∈ℂ^N and b ∈ℝ are the weights and bias, respectively. We can formally identify G_bw(x) with F_bλ(Ø) under the substitution

    ( x, w ) ( _ω (r|Ø),  _ω (r|λ) )   ,

where  is the discretization and vectorization to ℂ^N. This analogy is represented in <ref>.


A classical neuron requires at least N photons and N computational operations to classify an image composed of N pixels. Our setup bypasses both costs, by leveraging two essential features. On the one hand, it is completely optical, avoiding the computational need of processing the image. On the other hand, it classifies patterns through the Hong-Ou-Mandel effect, reducing the photon cost of imaging. In both ways, it provides a superexponential speedup, from 𝒪(N) to 𝒪(1). Photon losses due to absorption introduce a constant overhead in both the classical and quantum strategies, which depends on the total reflectivity of the object. We summarize this discussion in <ref>. See <ref> for a detailed derivation.




§ AMPLITUDE MODULATED PROBE

We specialize our discussion by replacing the generic probe state  with a toy model of an amplitude spatial light modulator (SLM), placed in the top branch of the Hong-Ou-Mandel interferometer, e.g. a liquid crystal (LC) grid with negligible losses <cit.>. Different approaches can be investigated, such as phase-only SLM <cit.>, which may exhibit superior resiliency against losses.


Consider a pattern on a greyscale LC grid with N real amplitudes {λ_μν}. Each pixel, labelled by (μ,ν), is represented by an L × L square with center r_μν = (μ + 1/2, ν + 1/2)L. Upon an overall parameter-independent normalization, the probe can be approximated as a combination of top-hat functions

    _ω(r|λ) = ∑_μ,ν u(r-r_μν)λ_μν/||λ|| ,

where ||λ||^2 = ∑_μ,νλ_μν^2 and u(r) := θ(r + L/2)-θ(r - L/2), with θ the two-dimensional Heaviside step function. Under this choice, <ref> simplifies to

    _λ(Ø) = | ∑_μ,ν (u ⋆_ω)(r_μν) λ_μν/||λ|||^2  ,

where ⋆ is the cross-correlation operation. We introduce a bias and a sigmoid activation function, so that the post-processed output reads _bλ(Ø) = σ(_λ(Ø) + b). Assuming that  is real, <ref> simplifies to

    ∂_μν≃ 2 √(f)/||λ||[(u ⋆_ω)(r_μν) - √(f)λ_μν/||λ||]  ,

with ∂_μν  := ∂ f / ∂λ_μν. This expression can be evaluated in an all-optical way, by taking the amplitude measurement of  directly in the left branch of the interferometer, before the beam splitter. This operation can be done off-line, and once per training object. In the next section, we present a simulation of these results, for different choices of the dataset.



 §.§ Simulations

We present a simulation of the model introduced above, comparing its performance against those of classical neural network-based techniques, for different datasets. All the simulations are run in Python and TensorFlow <cit.>, and summarized in <ref>.

We tested our model using two widely recognized datasets: the MNIST, which contains 28 × 28 images of handwritten digits from 0 to 9, and the CIFAR-10, comprised of 32 × 32 color images, distributed across 10 different classes. We guaranteed a fair comparison by increasing the MNIST resolution to 32 × 32 pixels (separately padding each image of the dataset), while converting the CIFAR-10 to greyscale. We represent each element of the dataset as (x_j,  y_j), where y_j ∈{0, 1} is the true class label, and x_j is the input vector, obtained by discretizing and vectorizing either the amplitudes  or their Fourier spectrum , thus bypassing the simulation of the imaging optics. We adopt the binary cross-entropy as loss function, combined with the standard (non-stochastic) gradient descent optimizer. We use the accuracy, i.e. the proportion of correct predictions over the total ones, as figure of merit of our results.

Our model demonstrates significative performances in both datasets (see <ref>). In the MNIST, it achieves accuracy rates exceeding 99%, when discerning between zeros and ones. In the CIFAR-10, it reaches accuracy above 58%, when distinguishing between cats and dogs. This difference reflects the complexity of the two classification tasks.

We compared our model against conventional neural network designs with a similar number of parameters. Specifically, we considered a single neuron and a convolutional neural network, commonly employed in pattern recognition tasks <cit.>. Adopting the TensorFlow notation, the convolutional structure is: Conv2D (10, 3 × 3) → Conv2D (4, 2 × 2) → MaxPooling2D (2 × 2). Roughly, all the architectures have ∼ 10^3 trainable parameters. All the models equally perform in the MNIST dataset, both in terms of trainability and final accuracy. When applied to the CIFAR-10 dataset, our classifier outperforms the conventional ones, showing superior efficiency under a strongly-constrained parameters count. All the findings emphasize the competitive accuracy of our method, and also its comparative advantage in pattern recognition tasks with a limited number of parameters.



§ CONCLUSIONS

In summary, we introduced an interferometric setup of a quantum optical classifier, with the Hong-Ou-Mandel effect as cornerstone of our classification method. We demonstrated the mathematical relation between our model and a classical neuron, constrained to unit depth, showing their similarity in terms of structure and response function. Our design is completely optical and single-photon based: it provides a superexponential speedup with respect to its classical counterpart, in terms of number of photons and computational resources. After modelling the classifier in terms of a spatial light modulator, we numerically compared our performances against those of standard neural network architectures, showing compatible to superior capabilities in terms of accuracy and training convergence, under the same number of parameters and depending on the pattern complexity.



§ ACKNOWLEDGEMENTS

S.R. acknowledges support from the PRIN MUR Project 2022SW3RPY. A.R.M. acknowledges support from the PNRR MUR Project PE0000023-NQSTI. C.M. acknowledges support from the National Research Centre for HPC, Big Data and Quantum Computing, PNRR MUR Project CN0000013-ICSC. L.M. acknowledges support from the PRIN MUR Project 2022RATBS4 and from the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Superconducting Quantum Materials and Systems Center (SQMS) under Contract No. DE-AC02-07CH11359. S.L. acknowledges support from ARO, DOE, and DARPA.



§ DATA AVAILABILITY

The underlying code that generated the data for this study is openly available in GitHub <cit.>.




§ SINGLE-PHOTON ENCODING

In this section, we consider the single-photon state obtained at the output of the left branch of the Hong-Ou-Mandel apparatus, providing a detailed discussion of <ref>. We adopt units in which c = 1.

Consider a generic single-photon state, generated by a monochromatic source with longitudinal position z

    |Ψ⟩ = ∫^3k  Φ̂(𝐤) a^†(𝐤)|0⟩ ,

with momentum spectrum Φ and 𝐤 = (k_x,k_y,k_z). We neglect the polarization of the photon and consider the single-frequency-mode assumption <cit.>, i.e. we assume that the wavefront propagates along definite-sign z-directions only. Then, k = (k_x,k_y) represents the only independent degrees of freedom of the single-photon state, which reads

    |Ψ⟩   = ∫^2k  ϕ̂_ω(k) a^†_ω(k)|0⟩ , 
       = ∫_S ^2r  ϕ_ω(r) a_ω^†(r)|0⟩ ,

where ϕ̂_ω(k) = Φ̂(k_x,k_y,√(ω^2-k_x^2-k_y^2)) and r = (r_x,r_y) labels the transverse coordinates on the source plane S.

For simplicity, we assume that the source is placed at the longitudinal origin z = 0. Consider an object with two dimensional shape Ø, placed at longitudinal position z_o. After free-space propagation occurs, the single-photon spectrum undergoes spatial amplitude modulation <cit.>, that is Ψ_Ø(r) = Ø(r) Ψ(r)_→, with Ψ(r)_→ the spatial input wavefront on the object plane O. Namely

    |Ψ_Ø⟩ = ∫_O ^2r  [ϕ_ω*ℌ_z_o](r)Ø(r) a^†_ω(r)|0⟩ ,

where ℌ_z_o denotes the free-space transfer function between the S and O planes. Using twice the convolution theorem, it follows that 

    |Ψ_Ø⟩ = ∫^2k  [(ϕ̂_ωℌ_z_o)*Ø̂](k) a^†_ω(k) |0⟩ .

Consider a linear optical system with transfer function ℒ, with image plane at longitudinal position z_i. By applying again the convolution theorem to _ω(·|Ø) = ((ϕ_ω*ℌ_z_o) Ø)*ℒ_z_o - z, we obtain

    |Ψ_⟩ = ∫^2k  _ω(k | Ø) a^†_ω(k) |0⟩ ,

with _ω(k|Ø) = [(ϕ̂_ωℌ_z_o)*Ø̂] ℒ̂_d. Notice that _ω(r|Ø) describes the image formed on a screen placed at distance d from the object.



§ HONG-OU-MANDEL COINCIDENCES

In this section, we compute the rate of coincidences at the output of the Hong-Ou-Mandel interferometer of <ref>, with left and top branch states given by <ref>, respectively. We write the input-probe bipartite state as

    |Ψ_⟩⊗|Ψ_⟩ = ∫^2k_1 ^2k_2  Ψ̂(k_1,k_2) a^†(k_1) b^†(k_2) |0⟩ ,

with Ψ̂(k_1,k_2) = (k_1 | Ø) (k_2|λ), where we dropped the ω subscript for simplicity. The 50:50 beam splitter acts as the unitary operation <cit.>

    a^†→1/√(2)( a^† + b^†) 
    
    		b^†→1/√(2)( a^† - b^†)
    	 ,
 
yielding

    |Ψ_⟩⊗|Ψ_⟩→|Φ⟩ = 1/2∫^2k_1 ^2k_2  Ψ̂(k_1,k_2) [ a^†(k_1) + b^†(k_1) ][ a^†(k_2) - b^†(k_2) ] |0⟩ .

Detection of mode m ∈{a,b} is described by the projector Π_m = ∫^2k  m^†(k)|0⟩⟨0|m(k). The rate of coincidences, i.e. the probability that one and only one photon is detected in each mode, reads

    p(1_a ∩ 1_b) = [|Φ⟩⟨Φ|Π_a ⊗Π_b]  , 
    with Π_a ⊗Π_b = ∫^2k_3 ^2k_4  a^†(k_3) b^†(k_4) |0⟩⟨0| a(k_3) b(k_4)  .

By substitution of <ref>, we get

    p(1_a ∩ 1_b) = 1/4∫∏_i=1^6 ^2k_i  Ψ̂(k_1,k_2) Ψ̂^*(k_5,k_6) W_1(k_1,k_2,k_3,k_4)W_2(k_3,k_4,k_5,k_6)  ,

where

    W_1(k_1,k_2,k_3,k_4)    = ⟨0| a(k_3) b(k_4)[a^†(k_1)a^†(k_2) - a^†(k_1)b^†(k_2) + b^†(k_1)a^†(k_2) - b^†(k_1)b^†(k_2) ] |0⟩
       = δ(k_2 - k_3)δ(k_1 - k_4) - δ(k_1 - k_3)δ(k_2 - k_4)  , 
    	
    
    		W_2(k_3,k_4,k_5,k_6)    = ⟨0|[a(k_6)a(k_5) - b(k_6)a(k_5) + a(k_6)b(k_5) - b(k_6)b(k_5) ] a^†(k_3) b^†(k_4) |0⟩
       = δ(k_3 - k_6)δ(k_4 - k_5) - δ(k_3 - k_5)δ(k_4 - k_6)  .

By integrating out the Dirac deltas in <ref>, we obtain

    p(1_a ∩ 1_b) = 1/2∫^2k_1 ^2k_2 ^2k_5 ^2k_6  Ψ̂(k_1,k_2) Ψ̂^*(k_5,k_6) [δ(k_1 - k_5)δ(k_2-k_6) - δ(k_1 - k_6)δ(k_2 - k_5) ]  .

Finally, the rate of coincidences reads

    p(1_a ∩ 1_b|λ, Ø) = 1/2∫^2k_1  |(k_1 | Ø)|^2 ∫^2k_2  |(k_2|λ)|^2 - 1/2| ∫^2k  (k | Ø) ^*(k|λ) |^2  .

More compactly, 

    p(1_a ∩ 1_b|λ, Ø) = 1/2[|| _ω(·|Ø) ||^2 || _ω(·|λ) ||^2 -  | ⟨_ω(·|Ø),_ω (·|λ) ⟩|^2 ]  ,

with || · || and ⟨· , ·⟩ denoting the L^2-norm and inner product, which is precisely the results of <ref>.



§ TRAINING

In this section, we discuss how to train the Hong-Ou-Mandel interferometer as a binary classifier. We separately feed each element of the training set (an ensemble of objects with known labels) into the input branch of <ref>, comparing the predicted classes with the target ones. We optimize the probe parameters λ by means of the gradient descent algorithm, and using the binary cross-entropy as loss function.

Consider a training set made of M objects {Ø_j}, each associated to a binary target label y_j ∈{0,1}, with 0 ≤ j ≤ M - 1. We denote ^(j)_λ = _λ(Ø_j) our model prediction. After feeding Ø_j into the input branch of the interferometer

    ^(j)_λ = C - 2p(1_a ∩ 1_b|λ, Ø_j)  , 
    _bλ^(j) = σ(^(j)_λ + b)  ,

where p ∈ [0, 1/2]. For simplicity, we assumed that the losses are independent on both the input and the probe, that is C := α_λ(Ø_j)  ∀λ, j.

Given a sample object, the binary cross-entropy between the target label and the predicted one reads

    H(y_j,_bλ^(j)) = - y_jlog (_bλ^(j)) - (1-y_j)log (1-_bλ^(j))  .
 
We optimize the probe parameters by means of the gradient descent algorithm, where the binary cross-entropy, averaged on the training set, is used as loss function. Namely

    λ→λ - η_λ/M∑_j=0^M-1∂_λ H(y_j,^(j)_bλ)  , 
    
    	b → b - η_b/M∑_j=0^M-1∂_b H(y_j,^(j)_bλ)  ,

with η_λ, η_b the learning rates of the probe and bias parameters, respectively. The derivatives with respect to the parameters and the bias yield

    ∂_λ H = (∂_ H )(∂_ξσ) ∂_λ , 
    ∂_b H = (∂_ H )∂_ξσ ,

with ξ_bλ = _λ + b. Then,

    ∂_ H = -y/(1-) , 
    ∂_ξσ = β (1 - )  ,

with β the hyperparameter of <ref>. For any complex function of real variable h: ℝ→ℂ, it follows that ∂_λ|h(λ)| = [h(λ)(∂_λ h(λ))^*]/|h(λ)|. Hence, 

    ∂_λ = 2 [ ⟨_ω, _ω⟩⟨_ω, ∂_λ_ω⟩^* ]  .
	
Neglecting the phase of ⟨, ⟩,

    ∂_λ≃ 2 √(f)[ ⟨_ω, ∂_λ_ω⟩]  .

This assumption, which we verified in our simulations under a self-consistency test, simplifies the computation of the first factor of <ref>, which is directly determined at the output of the Hong-Ou-Mandel interferometer.



§ CLASSIFICATION IN THE FOURIER DOMAIN

In this section, we discuss the effect of adding a single lens in the probe branch of the Hong-Ou-Mandel interferometer, as shown in <ref>. We summarize the main calculations, which closely follow that of <ref>.

A thin lens is placed at one focal length ℓ from both the probe image plane and the beam splitter. In the near-field limit, the lens performs a Fourier transform of the probe state <cit.>, yielding |Ψ_⟩→|Ψ_'⟩, where

    '_ω(r|λ) = -i ω/ℓ e^2i ω f_ω(ω/ℓr|λ)  .

After the beam splitter, the rate of coincidences is

    p(1_a ∩ 1_b|λ, Ø) = 1/2[α_ω(Ø) - _λ(Ø) ] 
         , 
    _λ(Ø) =  | ⟨_ω(·|Ø),_ω (·|λ) ⟩|^2  ,

yielding

    _bλ(Ø) = σ(_λ(Ø) + b)  , 
    _λ(Ø) = | ∫_I  r  r'  _ω(r|Ø)_ω^*(r'|λ)e^i r· r'|^2  ,

with σ and b the sigmoid activation function and bias, already introduced in <ref>. In contrast to <ref>, _λ(Ø) is not a point-wise evaluation: it combines the image spatial modes with the momentum spectrum of the probe state. Using the duality of the Fourier transform, it follows that

    _λ(Ø) =  | ⟨_ω(·|Ø),_ω (·|λ) ⟩|^2  ,

which corresponds to the output of the same scheme of <ref>, but with the thin lens placed in the left branch, before the beam splitter. Equivalently, this takes the Fourier transform of the image, instead of that of the probe. In the next section, we leverage this symmetry to simplify both the training process and the numerical simulations.

The training of the model follows the same procedure of <ref>. By placing the lens on the top branch of the interferometer, while using the duality of the Fourier transform, we get

    ∂_λ≃ 2 √()[ ⟨_ω, ∂_λ_ω⟩]  .

Under the same conditions of <ref>, the last two equations become

    _λ(Ø) = | ∑_μ,ν (u ⋆^*_ω)(r_μν) λ_μν/||λ|||^2 
    	 , 
    ∂_μν≃ 2 √()/||λ||[ (u ⋆^*_ω)(r_μν) - √()λ_μν/||λ||]   ,

where in the last step we neglected the phase of ⟨,⟩. Similarly to <ref>, <ref> can be evaluated in an all-optical way through the characterization of the real part of , namely, by performing an amplitude and phase measurement at the output of a thin lens, placed in the left branch, before the beam splitter. In <ref>, we compare the predictability of the neuron in the spatial and Fourier domains.



§ OPTICAL AND COMPUTATIONAL ADVANTAGE

In this section, we discuss the optical and computational advantage as the number of photons and operations required by a single image classification. Assuming that all the parameters have been previously trained with optimal accuracy, we show that our protocol requires a constant number of resources, i.e. 𝒪(1) complexity, independently of the input image resolution: it provides a superexponential speedup over its classical counterpart. 

We first discuss the computational advantage when substituting a classical neuron with a quantum optical one. From now on, we denote Ω, Θ and 𝒪, respectively the lower, tight and upper bounds on the number of resources needed by a certain (optical or computational) operation. Consider a digital image x of N pixels, fed into a neuron

    G_bw(x) = σ(w· x + b)  ,

where x, w ∈ℝ^N, b ∈ℝ and σ is the sigmoid  activation of <ref>, with hyperparameters β = 1 and γ = 0. <ref> costs N operations to compute w · x. The Hong-Ou-Mandel interferometer performs the same operation in an all-optical way, leaving the computational cost of the activation function and bias only, which is 𝒪(1).

We now discuss the optical advantage when using coincidences to classify single-photon states instead of a classical neuron on fully reconstructed images. After targeting an object with light, a digital image x is an ensemble of grey levels obtained by counting the number of photons collected by different pixels on a sensor grid, e.g. a charge-coupled device <cit.>. Let n_p be the average number of photons in the input state, and μ_i the average number of photons collected by the i-th pixel of the grid, with i ∈{0,…,N}. Assuming perfect quantum efficiency and sufficiently low exposure times to neglect the saturation of the sensor, the grey values at each pixel read

    x_i =  μ_iL/μ_w ,

with L the number of grey levels, i.e. the depth of the image, and μ_w = max_iμ_i the maximum number of photons collected in a single pixel. Indeed, x_i ∈{0, 1, …, L-1} with 0 and L labelling the black and white colors, respectively. Each pixel has variance ς_i^2 = Δμ_i^2 L^2/μ_w^2, with Δμ_i^2 the variance on the number of collected photons. For coherent light, the photo-detection process undergoes the standard quantum limit (SQL) <cit.>, with Poissonian fluctuations that satisfy Δμ_i^2 ≃μ_i. The average uncertainty reads

    ς^2 := 1/N∑_i=0^Nς_i^2 SQL≃⟨ x ⟩^2  N n_p^-1 ,

with ⟨ x ⟩ = N^-1∑_i x_i ∈ [0,L] the average brightness of the image, which we assume to be independent of its resolution. Hence, the number of photons n_p required by a full image reconstruction with average variance ς^2 is Θ(ς^-2⟨ x ⟩ N). This is the cost of image reconstruction only. We now take into account the information propagation through the neuron of <ref>.

Consider a neuron with sigmoid activation function. Suppose that there exists a sequence of parameters { (w_N , b_N) ∈ℝ^N+1}_N ≫ 1 that optimally solve the N-pixel image classification task, with b_N and the ℓ^1-norm ||w_N||_1 asymptotically bounded for N →∞. Then, the number of photons n_p required to classify an image x with uncertainty ε, is Ω(ε^-2⟨ x ⟩ N).


Consider the output of the neuron G_bw(x) = σ(w_N· x + b_N), and its derivative ∂ G(x) = G_bw(x)(1-G_bw(x)). By neglecting the spatial neighbourhood correlations, which may introduce at most a constant overhead in our estimation, we propagate the uncertainty of x as

    ε^2 = ⟨ x ⟩ (∂ G)^2(x) ∑_i=0^N-1(w_N)_i^2 x_i Nñ^-1_p  ,

where ñ_p = n_r n_p, with n_r is the number of independent image acquisition and classification. Since black pixels do not contribute to this summation, we get 

    ∑_i=0^N-1(w_N)_i^2 x_i ≥∑_i∉ℬ(w_N)_i^2 = ||w_N||^2 - ∑_i∈ℬ(w_N)_i^2  ,

with ℬ = {i∈ℕ |  x_i = 0  for 0 ≤ i ≤ N - 1 } the set of black pixels labels. However, ||w_N||^2 ≫∑_i∈ℬ(w_N)_i^2. Otherwise, ||w_N||^2 ≃∑_i∈ℬ(w_N)_i^2 would imply either that the image is mostly black, independently of its resolution, or that (w_N)_i ≃ 0 for all non-black pixels, which are both conditions that prevent the learnability of the neuron. By substitution into <ref> we get

    ñ_p ≥ε^-2⟨ x ⟩ (∂ G)^2(x) ||w_N||^2 N  .

Since w_N is a sequence of non-trivial solutions of the classification problem, the ℓ^2-norm ||w_N||^2 cannot go to zero for N →∞. Finally, we show that (∂ G(x))^2 does not converge to 0 for N →∞. Consider

    (∂ G)^2(x) = e^-2(w_N· x + b_N)/[1+e^-(w_N· x + b_N)]^4 .

If b_N is asymptotically limited, (∂ G)^2 converges to zero if and only if w_N· x →±∞. By splitting this scalar product into positive and negative contributions w_N· x = ∑_(w_N)_i > 0(w_N)_i x_i - ∑_(w_N)_i < 0|(w_N)_i| x_i, it follows that

    w_N· x ≤∑_(w_N)_i > 0(w_N)_i x_i ≤ L||w_N||_1  , 
       w_N· x ≥ - ∑_(w_N)_i < 0|(w_N)_i| x_i ≥ - L||w_N||_1  ,

namely that |w_N· x| ≤ L ||w_N||_1. Since the ℓ^1-norm is limited, (∂ G)^2 admits strictly positive lower bound for N →∞. Finally, this imply that ñ_p = Ω(ε^-2⟨ x ⟩ N).

In the previous discussion, two conditions lead to the above lower bound. On the one hand, that ||w_N||^2 ↛0 for N →∞, which is essential to guarantee that the neuron is trainable at any resolution. On the other hand, that ||w_N||_1 is bounded for N →∞, which is compatible with LASSO and Tikhonov's regularization techniques <cit.>.

We show that our protocol exponentially reduces this cost, requiring only the estimation of the rate of coincidences of the Hong-Ou-Mandel interferometer of <ref>. Let ñ_p = 2n_p be the number input photons, and p̃∈ [0,1/2] the empirical rate of coincidences. Under the normal approximation, with the 95% confidence level <cit.>, the estimation uncertainty reads

    ε = 2 √(p̃(1-p̃)/ñ_p) .

Since 4 p̃(1-p̃) ≤ 1, the total number of photons is Ø(ε^-2), which is constant with respect to the resolution of the image.

In conclusion, the quantum optical neuron provides a superexponential advantage over its classical counterpart, both in the number of operations and photons saved to classify a single image. We summarize these results in <ref>.




