Zakari Denis

Ivan Favero

Cristiano Ciuti
email: cristiano.ciuti@u-paris.fr

Photonic kernel machine learning for ultrafast spectral analysis

We introduce photonic kernel machines, a scheme for ultrafast spectral analysis of noisy radiofrequency signals from single-shot optical intensity measurements. The approach combines the versatility of machine learning and the speed of photonic hardware to reach unprecedented throughput rates. We theoretically describe some of the key underlying principles, and then numerically illustrate the reached performances on a photonic lattice-based implementation. We apply the technique both to picosecond pulsed radio-frequency signals, on energy-spectral-density estimation and a shape classification task, and to continuous signals, on a frequency tracking task. The presented optical computing scheme is resilient to noise while requiring minimal control on the photonic-lattice parameters, making it readily implementable in realistic state-of-the-art photonic platforms.

I. INTRODUCTION

As a result of its intrinsically faster timescales, photonics was very soon envisioned as a promising tool to outperform integrated electronics in terms of data processing rates [START_REF] Ambs | Optical Computing: A 60-Year Adventure[END_REF]. In this perspective, optical setups ranging from matrix-vector multipliers [START_REF] Mengert | Control apparatus[END_REF], function convolvers [START_REF] Weaver | A Technique for Optically Convolving Two Functions[END_REF] and discrete Fourier-transforming processors [START_REF] Goodman | Fully parallel, high-speed incoherent optical method for performing discrete Fourier transforms[END_REF] to non-von Neumann parallel digital processors [START_REF] Sawchuk | Digital optical computing[END_REF][START_REF] Huang | Architectural considerations involved in the design of an optical digital computer[END_REF] were proposed.

Although these ideas quickly became outdated as a consequence of the fast rise of silicon-based electronic processors, the latter have started to exhibit some of their limitations. In particular, the emergence of machine learning applications operating on ever increasing amounts of data, involving deeper and deeper neuralnetwork architectures with an increasing degree of complexity [START_REF] Khan | A survey of the recent architectures of deep convolutional neural networks[END_REF] has led to a situation where the progress of available digital processor technology no longer keeps pace with the demand in computing capabilities [START_REF] Waldrop | The chips are down for Moore's law[END_REF]. This trend has moreover gone hand in hand with an increase in the consumption of computational and energy resources, casting doubt on its sustainability [START_REF] Strubell | Energy and Policy Considerations for Deep Learning in NLP[END_REF]. This context has stimulated very interesting proposals aiming at surrogating the realization of specific tasks that are computationally and energetically very demanding to very specialized (electro-)optical devices. This process, reminiscent of current trends in hardware acceleration technologies, such as graphical (GPUs) and tensor processing units (TPUs), has already led to commercially available optical co-processors [START_REF] Dong | Scaling Up Echo-State Networks With Multiple Light Scattering[END_REF][START_REF] Ohana | Kernel Computations from Large-Scale Random Features Obtained by Optical Processing Units[END_REF][START_REF] Schneider | Lightmatter's Mars Chip Performs Neural-Network Calculations at the Speed of Light: MIT spinoff harnesses optical computing to make neural-networks run faster and more efficiently[END_REF][START_REF] Schneider | Deep Learning at the Speed of Light: Lightmatter bets that optical computing can solve AI's efficiency problem[END_REF]. In the recent years, this approach was scaled to deep architectures [START_REF] Wetzstein | Inference in artificial intelligence with deep optics and photonics[END_REF] and has proven spectacularly powerful in the field of computer vision, by exploiting diffraction [START_REF] Lin | All-optical machine learning using diffractive deep neural networks[END_REF] or by optical implementations of convolutional neural networks, both in free-space [START_REF] Miscuglio | Massively parallel amplitude-only Fourier neural network[END_REF] and on chip [START_REF] Feldmann | Parallel convolutional processing using an integrated photonic tensor core[END_REF]. Such architectures are able to extract increasingly abstract representations [START_REF] Zeiler | Visualizing and Understanding Convolutional Networks[END_REF] of the input images fed into the network by subsequent applications of pooled optically-operated linear convolutions.

The most standard neural-network-based machine learning schemes present roughly the following architecture: the digital data to be processed are transformed by a series of consecutive layers that consist in a parametrized affine transformation followed by the pointwise application of an elementary nonlinear activation function. Such a composition of parametrized functions is then expected to approximate some target function of the input upon a proper training process involving the optimization of the parameters in order to minimize the error of the model over some set of training examples. While this sequential architecture is ideally suited for standard processors, which offer arbitrary levels of programmability by design, the amount of parameters to be addressed during the training process is in practice a hurdle to flexible optical implementations, having progressed from roughly 63 thousands in the paradigmatic LeNet-5 convolutional neural network [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF] to several tens of millions on its nowadays counterparts [START_REF] Khan | A survey of the recent architectures of deep convolutional neural networks[END_REF].

Therefore, new machine learning paradigms that relax the above constraints, such as extreme-learning machines [START_REF] Huang | Extreme learning machine: Theory and applications[END_REF][START_REF] Huang | Extreme Learning Machine for Regression and Multiclass Classification[END_REF], echo-state networks or reservoir computing [START_REF] Jaeger | The "echo state" approach to analysing and training recurrent neural networks[END_REF][START_REF] Jaeger | Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication[END_REF][START_REF] Maass | Real-Time Computing Without Stable States: A New Framework for Neural Computation Based on Perturbations[END_REF][START_REF] Verstraeten | Isolated word recognition with the Liquid State Machine: a case study[END_REF][START_REF] Verstraeten | Reservoir-based techniques for speech recognition[END_REF][START_REF] Lukoševičius | Reservoir computing approaches to recurrent neural network training[END_REF][START_REF] Rodan | Minimum Complexity Echo State Network[END_REF][START_REF] Tanaka | Recent advances in physical reservoir computing: A review[END_REF][START_REF] Nakajima | Physical reservoir computing-an introductory perspective[END_REF][START_REF] Marcucci | Theory of Neuromorphic Computing by Waves: Machine Learning by Rogue Waves, Dispersive Shocks, and Solitons[END_REF], have inspired theoretical proposals and experimental realizations in a variety of settings, in optics [START_REF] Der Sande | Advances in photonic reservoir computing[END_REF][START_REF] Sunada | Using multidimensional speckle dynamics for high-speed, large-scale, parallel photonic computing[END_REF][START_REF] Pierangeli | Photonic extreme learning machine by free-space optical propagation[END_REF], integrated photonics [START_REF] Vandoorne | Experimental demonstration of reservoir computing on a silicon photonics chip[END_REF][START_REF] Denis-Le Coarer | All-Optical Reservoir Computing on a Photonic Chip Using Silicon-Based Ring Resonators[END_REF], memristors [START_REF] Kulkarni | Memristor-based reservoir computing[END_REF][START_REF] Du | Reservoir computing using dynamic memristors for temporal information processing[END_REF] and beyond [START_REF] Boyn | Learning through ferroelectric domain dynamics in solid-state synapses[END_REF][START_REF] Nakane | Reservoir Computing With Spin Waves Excited in a Garnet Film[END_REF][START_REF] Marković | Reservoir computing with the frequency, phase, and amplitude of spintorque nano-oscillators[END_REF]. In the context of optics, many fruitful configurations have been investigated, such as delay-line-based setups [START_REF] Vandoorne | Toward optical signal processing using Photonic Reservoir Computing[END_REF][START_REF] Paquot | Reservoir computing: A photonic neural network for information processing[END_REF][START_REF] Appeltant | Information processing using a single dynamical node as complex system[END_REF][START_REF] Vandoorne | Parallel Reservoir Computing Using Optical Amplifiers[END_REF][START_REF] Duport | All-optical reservoir computing[END_REF][47][START_REF] Larger | Photonic information processing beyond Turing: An optoelectronic implementation of reservoir computing[END_REF][START_REF] Soriano | Optoelectronic reservoir computing: Tackling noise-induced performance degradation[END_REF][START_REF] Nguimdo | Simultaneous Computation of Two Independent Tasks Using Reservoir Computing Based on a Single Photonic Nonlinear Node With Optical Feedback[END_REF][START_REF] Ortín | A Unified Framework for Reservoir Computing and Extreme Learning Machines based on a Single Time-delayed Neuron[END_REF][START_REF] Larger | High-Speed Photonic Reservoir Computing Using a Time-Delay-Based Architecture: Million Words per Second Classification[END_REF][START_REF] Chembo | Machine learning based on reservoir computing with time-delayed optoelectronic and photonic systems[END_REF], nonlinear polariton lattices [START_REF] Opala | Neuromorphic Computing in Ginzburg-Landau Polariton-Lattice Systems[END_REF][START_REF] Ballarini | Polaritonic Neuromorphic Computing Outperforms Linear Classifiers[END_REF][START_REF] Mirek | Neuromorphic Binarized Polariton Networks[END_REF] and systems combining linear light scattering and the measurement nonlinearity [START_REF] Pierangeli | Photonic extreme learning machine by free-space optical propagation[END_REF][START_REF] Dong | Optical Reservoir Computing Using Multiple Light Scattering for Chaotic Systems Prediction[END_REF][START_REF] Rafayelyan | Large-Scale Optical Reservoir Computing for Spatiotemporal Chaotic Systems Prediction[END_REF]. Beyond liquid-state machines/echo-state networks, several proposals have been put forward to reconcile generalpurpose machine learning and distributed memory architectures. This, in particular, is the realm of neuromorphic computing, that has led to significant work in the field of spintronics [START_REF] Sengupta | Proposal for an All-Spin Artificial Neural Network: Emulating Neural and Synaptic Functionalities Through Domain Wall Motion in Ferromagnets[END_REF][START_REF] Grollier | Spintronic Nanodevices for Bioinspired Computing[END_REF][START_REF] Torrejon | Neuromorphic computing with nanoscale spintronic oscillators[END_REF][START_REF] Romera | Vowel recognition with four coupled spintorque nano-oscillators[END_REF][START_REF] Grollier | Neuromorphic spintronics[END_REF] and memristors [START_REF] Jo | Nanoscale Memristor Device as Synapse in Neuromorphic Systems[END_REF][START_REF] Yu | An Electronic Synapse Device Based on Metal Oxide Resistive Switching Memory for Neuromorphic Computation[END_REF][START_REF] Yang | Memristive Physically Evolving Networks Enabling the Emulation of Heterosynaptic Plasticity[END_REF][START_REF] Kim | Experimental Demonstration of a Second-Order Memristor and Its Ability to Biorealistically Implement Synaptic Plasticity[END_REF][START_REF] Chu | Neuromorphic Hardware System for Visual Pattern Recognition With Memristor Array and CMOS Neuron[END_REF][START_REF] Prezioso | Training and operation of an integrated neuromorphic network based on metal-oxide memristors[END_REF][START_REF] Jeong | Memristors for Energy-Efficient New Computing Paradigms[END_REF][START_REF] Hsieh | A sub-1-volt analog metal oxide memristive-based synaptic device with large conductance change for energy-efficient spike-based computing systems[END_REF][START_REF] Wang | Memristors with diffusive dynamics as synaptic emulators for neuromorphic computing[END_REF][START_REF] Yoon | An artificial nociceptor based on a diffusive memristor[END_REF][START_REF] Wang | Fully memristive neural networks for pattern classification with unsupervised learning[END_REF][START_REF] Li | Efficient and self-adaptive in-situ learning in multilayer memristor neural networks[END_REF][START_REF] Kim | Nociceptive Memristor[END_REF][START_REF] Jeong | Nonvolatile Memory Materials for Neuromorphic Intelligent Machines[END_REF], notably that articulated around spike-based machinelearning schemes [START_REF] Afifi | Implementation of biologically plausible spiking neural network models on the memristor crossbar-based CMOS/nano circuits[END_REF][START_REF] Querlioz | Simulation of a memristor-based spiking neural network immune to device variations[END_REF][START_REF] Querlioz | Immunity to Device Variations in a Spiking Neural Network With Memristive Nanodevices[END_REF][START_REF] Srinivasan | Magnetic Tunnel Junction Based Long-Term Short-Term Stochastic Synapse for a Spiking Neural Network with On-Chip STDP Learning[END_REF][START_REF] Roy | Towards spike-based machine intelligence with neuromorphic computing[END_REF][START_REF] Sengupta | Going Deeper in Spiking Neural Networks: VGG and Residual Architectures[END_REF]. Such architectures come with their own technical difficulties. Because they rely on co-located memory and processing resources, standard optimization algorithms can prove impractical, although novel optimization procedures have recently been put forward to overcome this obstacle by meeting their peculiar hardware design [START_REF] Ernoult | Updates of Equilibrium Prop Match Gradients of Backprop Through Time in an RNN with Static Input[END_REF][START_REF] Ernoult | Equilibrium Propagation with Continual Weight Updates[END_REF][START_REF] Laborieux | Scaling Equilibrium Propagation to Deep ConvNets by Drastically Reducing Its Gradient Estimator Bias[END_REF][START_REF] Laydevant | Training Dynamical Binary Neural Networks with Equilibrium Propagation[END_REF][START_REF] Martin | EqSpike: Spike-driven equilibrium propagation for neuromorphic implementations[END_REF].

A second shortcoming of standard software machine learning relates to applications involving analog data that cannot be suitably interfaced with digital processors. Instances of this problem are, for example, situations where the input data to be analyzed are supplied at a throughput too high to be properly sampled in real time. This also occurs when the data are intrinsically analog, or when direct measurement processes add noise or perturb the system being measured. The utmost example of this is provided by genuinely quantum tasks with no classical counterparts, that involve quantum inputs. This has stimulated many original works, ranging from quantum metrology [START_REF] Ghosh | Quantum reservoir processing[END_REF][START_REF] Ghosh | Reconstructing Quantum States With Quantum Reservoir Networks[END_REF] and quantum state control [START_REF] Ghosh | Quantum Neuromorphic Platform for Quantum State Preparation[END_REF][START_REF] Ghosh | Realising and compressing quantum circuits with quantum reservoir computing[END_REF][START_REF] Krisnanda | Creating and concentrating quantum resource states in noisy environments using a quantum neural network[END_REF] to classical image recognition tasks with quantum hardware [START_REF] Xu | Superpolynomial quantum enhancement in polaritonic neuromorphic computing[END_REF], under the name of quantum neuromorphic computing [START_REF] Marković | Quantum neuromorphic computing[END_REF],

Here, we take a different path and propose to implement a distinct machine-learning paradigm, namely kernel machines [START_REF] Boser | A training algorithm for optimal margin classifiers[END_REF][START_REF] Schölkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF], in the world of photonics. The approach appears to circumvent the two difficulties listed above. We describe photonic kernel machines as well as the associated theoretical framework under very general assumptions. We explore the links between key concepts in support-vector-machine (SVM) theory and those of the machine proposed here. We show that, in contrast with general reservoir computing schemes, knowledge about the underlying internal representations may be revealed from measurable data, providing direct understanding of the learning process of actual hardware. By introducing a realistic physical model for a photonic kernel machine based on a two-dimensional lattice of coupled linear optical cavities, we numerically examine the performance of photonic kernel machines on regression and classification tasks involving ultrafast spectral analysis of analog and noisy radio-frequency (RF) signals, which are imprinted on an optical carrier wave.

This paper is organized as follows. After presenting some general concepts of kernel-machine theory in Section II, photonic kernel machines are theoretically introduced in Section III; their learning mechanism is discussed therein. A model for a physical implementation based on a photonic lattice is then described in Section IV. This model is numerically simulated in Section V and applied to the ultrafast spectral analysis of noisy radio-frequency signals from single-shot optical intensity measurements. Finally, conclusions are drawn in Sec. VI.

II. GENERAL FRAMEWORK

Supervised machine-learning problems can be formulated as an optimization problem in which one tries to best approximate a target quantity y = f (x) of some given input x = [x 1 , x 2 , . . .] T with a parametrized function f . The input data are distributed according to some unknown distribution p(x) from which only a restricted set of examples {(x (i) , y (i) )} N i=1 is known.

More precisely, in regression problems one seeks for the optimal parameters that make f the best fit for the known examples by approximating their unknown true functional dependence f . At variance, in classification problems one aims at determining a predictor f that best associates to any given input x a set of labels y that characterizes its belonging to one or more classes. These approximation problems are expressed in practice as the minimization of some cost function J with respect to the parameters of the model to be trained, for a given set of examples.

A specific choice of parametrization defines the architecture of the model. Many such architectures exist, ranging from shallow models [START_REF] Hastie | The Elements of Statistical Learning[END_REF], such as support vector machines (SVM), tree-based models and reservoir computing, to deep-learning models [START_REF] Goodfellow | Deep Learning[END_REF], such as feedforward, convolutional or recurrent neural networks. Here, we will exclusively consider shallow models, with a strong emphasis towards kernel machines, which fundamental concepts and results will be introduced below. These notions will be particularly useful in the understanding of the photonic kernel machine discussed in Sec. III and the associated physical-implementation model presented in Sec. IV.

A. Shallow models

The principle underlying the learning process of shallow models is sketched in Fig. 1. Inputs are first embedded from their original input space into a typically higherdimensional feature space where the optimization problem becomes linear. The optimization process is then realized therein, for instance by means of standard convex optimization algorithms, and often reduces to the geometrical problem of identifying an optimal hyperplane.

More formally, the trial function f can be expressed in terms of a set of transformed (feature-space) coordinates as

f (x) = M m=1 w m xm + b = w T x + b, (1)
where a higher-dimensional feature space spanned by nonlinear transformations of the inputs. The feature-space embedding maps the original datapoints onto a hyperplane therein. The model of Eq. ( 1) can thus perform regression by approximating the target function as a plane parametric curve x = ψ(x) lying on this possibly infinitedimensional hyperplane, whose parameters (w, b) are to be determined.

x := ψ(x), (2)
with orthogonal components ψ m , ψ n p = E p [ψ m (x)ψ n (x)] ∝ δ m,n , and E p [ψ m (x)ψ n (x)] := dxp(x)ψ m (x)ψ n (x).
Binary classification can be operated in a similar way. Inputs x are associated to some class a if f (x) > 0 and to the complementary class b otherwise. While this function might be involuted in input space, it may also become linear in some suitable feature space, the frontier between the two classes, f (x) = w T x + b = 0, thus becoming the equation of a plane. In this situation, the objective of the optimization becomes finding the proper plane in feature space that separates data belonging to the two distinct classes into two separate clusters, as illustrated in Fig. 1 (b). This linear decision boundary translates back into a potentially non-trivial one in input space [dashed circle in Fig. 1 (b)]. The same mechanism can be exploited to perform k-class classification via the one-vsone and one-vs-rest techniques, which split the problem into respectively k(k -1)/2 and k binary-classification problems.

Such a shallow architecture presents a number of advantages. First, the optimization process is most often convex, in many cases analytic. Second, because it relies on a fixed feature-space embedding process, independent from any parameter to be updated during the training process, this architecture can be implemented on fast physical hardware, thus lifting any further analog-todigital interfacing overhead.

Many possible choices exist to realize the above feature-space embedding. One popular choice is that of reservoir computing [START_REF] Tanaka | Recent advances in physical reservoir computing: A review[END_REF][START_REF] Nakajima | Physical reservoir computing-an introductory perspective[END_REF]. This scheme exploits the nonlinear response of a dynamical system with many de-grees of freedom-the reservoir -to a driving input that encodes the data to be analyzed. As the set of produced independent nonlinear outputs for a single input is increased, the spanned feature space is expected to "include" the optimal one where the problem becomes linear [START_REF] Huang | Extreme learning machine: Theory and applications[END_REF][START_REF] Huang | Extreme Learning Machine for Regression and Multiclass Classification[END_REF]. While no control on the encoding of the inputs and the parameters of the reservoir is in principle required, in practice, this may impact the convergence of the performance of the model in the number of extracted nonlinear outputs.

In the following, a similar yet alternative approach will be considered, rooted in the theory of support-vector machines (SVM): the so-called kernel machines. These rely on an educated guess : the feature-space embedding is induced based upon a choice of similarity metric, encoded in a kernel K, which associates to any two given inputs x and x a measure of their similarity K(x, x ). Beyond providing heuristics on the way of building a featurespace embedding best suited for a specific task, such models come with theoretical bounds on the probability of generalization error [START_REF] Schölkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF]101,[START_REF] Mohri | Adaptive Computation and Machine Learning Series[END_REF].

In practice, kernel-method practitioners choose a kernel suitable for their applications; among popular choices are linear kernels, K(x, x ) = x T x ; polynomial kernels, K(x, x ) = (1 + x T x /c) d ; radial basis functions (RBF) [START_REF] Powell | Radial basis functions for multivariable interpolation : A review[END_REF][START_REF] Broomhead | Multivariable Functional Interpolation and Adaptive Networks[END_REF], K(x, x ) = expxx 2 2 /σ 2 ; or sigmoid kernels, K(x, x ) = tanh(κx T x + θ); to name a few [START_REF] Van Gestel | Benchmarking Least Squares Support Vector Machine Classifiers[END_REF]. Very recently, "quantum" kernels of the form K(x, x ) = Tr ρ(x) † ρ(x ) or K(x, x ) = Â(x) † Â(x ) , where ρ(x) and Â(x) denote quantum operators evolved through some input-dependent unitary transformation Û [x] were proposed [START_REF] Schuld | Quantum Machine Learning in Feature Hilbert Spaces[END_REF][START_REF] Bartkiewicz | Experimental kernel-based quantum machine learning in finite feature space[END_REF] and even experimentally implemented [START_REF] Bartkiewicz | Experimental kernel-based quantum machine learning in finite feature space[END_REF][START_REF] Havlíček | Supervised learning with quantum-enhanced feature spaces[END_REF][START_REF] Kusumoto | Experimental quantum kernel trick with nuclear spins in a solid[END_REF]. Quantum advantage on a classification task was demonstrated on such kernel machines [START_REF] Liu | A rigorous and robust quantum speed-up in supervised machine learning[END_REF]; more recently, it was shown that circuitbased "quantum neural networks" could be described as such quantum kernel machines [START_REF] Schuld | Supervised quantum machine learning models are kernel methods[END_REF].

Although in principle large-scale data processing with kernels can be computationally very expensive, several approximate techniques have been developed to overcome this limitation [START_REF] Williams | Using the Nyström Method to Speed Up Kernel Machines[END_REF][START_REF] Drineas | On the Nystrom Method for Approximating a Gram Matrix for Improved Kernel-Based Learning[END_REF][START_REF] Rahimi | Random features for largescale kernel machines[END_REF][START_REF] Rahimi | Weighted sums of random kitchen sinks: Replacing minimization with randomization in learning[END_REF], some of which have recently inspired physical implementations capable of efficiently performing approximate kernel evaluations of high-dimensional vectorial data (e.g. images) on optical hardware [START_REF] Ohana | Kernel Computations from Large-Scale Random Features Obtained by Optical Processing Units[END_REF][START_REF] Saade | Random projections through multiple optical scattering: Approximating Kernels at the speed of light[END_REF]; see also Ref. [START_REF] Pierangeli | Photonic extreme learning machine by free-space optical propagation[END_REF] for interesting related strategies using free-space photonics. In the present work, we both define a kernel suitable for processing (infinite-dimensional) continuous time signals and provide a way to efficiently approximate it with photonic hardware.

B. Kernel machines

We here recall some fundamental identities from the theory of kernel machines that will prove useful in the following.

Upon making a specific choice of symmetric positive semi-definite kernel K, the kernel machine's predictions f can be evaluated in two different forms. First, as an average over a set of trainable model parameters α = [α 1 , α 2 , . . . , α N ] T , weighted over the similarity between the considered input x and all those composing the training set, of size N :

f (x) = N i=1 α i K(x, x (i) ) + b, (3)
where b is a bias parameter. Alternatively, predictions can be obtained as an expansion over some set of

M ≤ +∞ possibly-nonlinear functions {h m } M m=1 such that K(x, x ) = m h m (x)h m (x ), ∀(x, x ). Explicitly: f (x) = M m=1 β m h m (x) = β T h(x), (4)
where

β = [β 1 , . . . , β M ] T , of size M , is the set of param- eters of the model [117].
The optimal parameters of the model ( α and β, respectively) are determined by minimizing a cost function. Each approach leads to a different yet equivalent optimization problem, respectively known as the dual and the primal problem, of the form:

J(ξ) = N i=1 V y (i) , f (x (i) ) + λR(ξ), (5)
with either ξ = α (dual problem) or ξ = β (primal problem). Here, V (y (i) , f (x (i) )) is some pointwise error function on the predictions made by the model that one wants to minimize and R a regularization term whose strength is controlled by means of a "bias" hyperparameter λ.

Regularization prevents the trained model from overfitting the dataset, thus increasing its generalization power. In the following, ridge regularization will be considered, defined by

R(α) = 1 2 α T Kα, R(β) = 1 2 β 2 ; (6)
where

K ij = K(x (i) , x (j)
) denotes the kernel matrix.

Both the dual and the primal approaches can be shown to lead to a same feature-space embedding, of the form of Eq. ( 1), with feature maps given by x → ψ m (x) = √ γ m φ m (x), where γ m and φ m simply derive from the kernel's eigendecomposition:

K(x, x ) = M m=1 γ m φ m (x)φ m (x ), (7)
with

γ m+1 ≤ γ m and φ m , φ n p = E p [φ m (x)φ n (x)] ≡ dxp(x)ψ m (x)ψ n (x) = δ m,n
. This eigendecomposition can be given empirical estimates by approximating the actual probability density distribution of inputs p(x) by the empirical one p(x) [START_REF] Williams | Using the Nyström Method to Speed Up Kernel Machines[END_REF], such that dxp(x)f (x) = (1/N ) i f (x (i) ). As shown in the appendix, given a primal problem parametrized by {h m } M m=1 , the empirical eigenvalues {γ m } rank(K)≤M m=1 are simply given by those of the matrix k/N = H T H/N , where H im = h m (x (i) ), with the following associated empirical eigenfunctions:

φm (x) = 1 √ N γm u T m h(x), (8)
where u m denotes k's mth eigenvector. These altogether lead to the empirical feature map ψm (x) = √ γm φm (x). The ability of kernel machines to filter out noise in input space can now be understood in very simple terms. To this aim, let us start from the primal approach and consider a kernel machine whose predictions are obtained as in Eq. ( 4), with outputs centered by the empirical mean over the training set, i.e. replacing h m (x) by hm

(x) = h m (x) -E p[h m (x)], where E p[X (x)] = (1/N ) i X(x (i) )
, for any X. It then follows that the corresponding kernel's empirical eigenfunctions are centered as well, (1/N ) i φm (x (i) ) = 0, and the eigendecomposition ( 7) is in all respects analogous to principal component analysis on the feature-space-embedded data inputs [START_REF] Schölkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF][START_REF] Schölkopf | Kernel principal component analysis[END_REF][START_REF] Mika | Kernel PCA and de-noising in feature spaces[END_REF]. In this picture, schematically illustrated in Fig. 2, the empirical eigenvalues exactly correspond to the variance of the data along every featurespace direction:

γm = E p ψm (x) -E p[ ψm (x)] 2 . ( 9
)
Large variances of this type correspond to very correlated features that maximally deviate between inputs considered as "dissimilar" according to one's initial choice of similarity metric. Conversely, low variances correspond either to redundant features or to uncorrelated noise contained in the inputs.

By now decomposing the trial function in the kernel basis identified above:

f (x) = m φm , f p φm (x) ≡ m fm φm (x), (10)
where the components fm may be interpreted as degrees of freedom of the model, it follows from Eq. ( 6) that the regularization penalty can be rewritten as

λR(α) = m λ m 2 | fm | 2 , λ m = λ γm . ( 11
)
This means that the penalty on the magnitude of the mth component of f during the optimization is inversely proportional to the variance of the associated feature-space coordinate ψ m (x). Therefore, the effect of the regularization on the optimization is to impose a soft cut-off on the dimensionality of the feature-space hyperplane schematically illustrated in Fig. 1, while preserving most of the original variance contained in the inputs by only keeping the most statistically relevant directions, expected to encode the correlated information of interest in the inputs, while filtering out the noise [START_REF] Mika | Kernel PCA and de-noising in feature spaces[END_REF]. This can be understood as a filtering process as will clearly appear in the following.

Feature space

√ γ1 √ γm 1 ψ1 ψm 1 FIG.
2. Schematic representation of the projection of the feature-space-embedded data onto the plane spanned by the leading principal component of some kernel K and some subleading one. The variance stemming from the correlations of the embedded training examples is expectected to be mostly contained within the dimensions spanned by the kernel's principal components; that of the noise, along orthogonal directions.

III. PHOTONIC KERNEL MACHINES

The spectral analysis of radio-frequency signals is a very broad field with countless applications. We notably focus on the analysis of ultrashort pulsed signals. This is relevant for instance in the field of pulse-Doppler radars. Such devices emit a pulsed radio-frequency signal around a carrier frequency with a given repetition rate, and receive it back via an antenna after being reflected by a moving target. This reflected signal is measured and numerically processed to extract information about the position and the velocity of the tracked target. The repetition rate is limited by the processing time, which is in turn bounded by the acquisition time. Nowadays, the state of the art performances are obtained by direct sampling of the signal sensed by the antenna. Yet, the sampling rate of state-of-the-art RF analogue-to-digital converters lies typically below 10 gigasamples per second [START_REF]See for instance Texas Instruments' ADC12DL3200 (6.4 GSPS) or Analog Devices[END_REF]. In order to measure 1000 samples of a pulse of interest with such a converter, one needs the pulse length to exceed 25 µs, limiting our ability to analyze shorter pulses. This issue may be circumvented by operating several such samplers in parallel. However, the resulting devices are very expensive, heavy and bulky. In the optical frequency domain, spectral analysis can be performed by means of monochromators, which perform a frequency sweep with a power measurement at the selected frequency. The control of the mechanical element allowing the frequency sweep greatly limits the spectrum acquisition rate, which can be at best of the order of 10 kHz.

This context motivates the quest for technologies that go past these limitations. In what follows, we propose to resort to a learning all-optical processing device based upon the kernel-machine concept. We will first introduce a similarity kernel ideally suited for such spectral analysis tasks. We will show that this kernel can be evaluated rather naturally by means of an optical single-shot intensity measurement at the output of a photonic lattice. Such a photonic kernel machine makes it possible to process complex radio-frequency signals in times of the order of tens or hundreds of picoseconds.

A. Theoretical description

As we have seen above, learning with kernels starts from the introduction of a kernel that serves as a measure of similarity between inputs. In particular, we here focus the discussion on the spectral analysis of ultrashort pulsed RF signals. The spectral information of such signals s(t) is encoded in their energy spectral density

S[ω] = |s[ω]| 2 , where s[ω] = (2π) -1/2 dte -iωt s(t).
Rather naturally, the similarity between two signals s(t) and s (t) can thus be given a general expression of the form:

K(S, S ) = dωdω S[ω]K(ω -ω )S [ω ], (12)
where K(ωω ) is a function peaked around ω = ω and with a typical width δω. Such a kernel compares two given input signals by contrasting their energy spectral densities at each frequency with a certain tolerance on their fine structure at frequency scales below δω. While this provides a good and flexible similarity metric, the numerical evaluation of such a kernel is rather unsuitable in practice using current electronics. Indeed, this would require several costly steps: (i) each pulse would have

s (i) (t) ∆t dω|s[ω]| 2 × . . . dω|s (i) [ω]| 2 × { ψ } Feature-space embedding . . . • . . . { ŵ } Linear regression ŷ(i) [ω]
FIG. 3. Illustration of the photonic-kernel-machine processing mechanism. The measurement of the pulse-induced optical populations embeds the pulse's energy spectral density into some reciprocal space by projection onto a set of orthogonal base functions { ψ } . The spectrum of the ith pulse is reconstructed by a linear combination ŷ(i

) [ω] = w [ω] ψ |S (i)
of its reciprocal-space components (features) over some set of learned functions { ŵ } . On this example, the learned functions { ŵ } are filtered analogues of { ψ } and the photonic kernel machine is able to extract the spectrum of the incoming pulse from its noisy background.

to be sampled in time at similar sampling rates over a time interval larger than 2π/δω, (ii) the digitized signals would then have to be numerically Fourier-transformed and stored, (iii) finally, each kernel evaluation would involve calculating a double integral. In contrast, we will see that the feature-space embedding of this kernel can be performed rather naturally on photonic hardware. To that purpose, the RF signal to be analyzed can first be imprinted on an optical carrier injected into the hardware.

Let us consider an optical system consisting of a set of M generic normal modes as described by the set of bosonic operators {α } M =1 . Their response to any ith pulsed input signal s (i) (t) of interest is well described by a quantum Langevin equation [START_REF] Gardiner | Input and output in damped quantum systems: Quantum stochastic differential equations and the master equation[END_REF][START_REF] Gardiner | Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics[END_REF]. In frequency space, this takes the form:

α [ω] = χ [ω] is (i) [ω] + κ/2 αin [ω] , (13)
where α [ω] = (2π) -1/2 dte -iωt α (t) is the Fouriertransformed annihilation operator associated to the th normal mode, χ [ω] its (model-dependent) susceptibility at frequency ω, and κ the loss rate, taken independent of . Here, the Langevin input fields αin (t) account for zero-temperature quantum noise and satisfy the usual expectation values:

αin (t) = αin † (t)α in (t ) = 0 and [α in (t), αin † (t )] = δ , δ(t -t ).
Optical populations in the modes, induced by the input pulse, are measured via the radiated optical power collected by a detector, with some integration time ∆t much larger than the length of the pulse, ∆t 2π/∆ω, where ∆ω is the bandwidth of the signal to be analyzed. In this limit, the measured normal-mode populations are expressed as

n(i) 1 ∆t dt α † (t)α (t) = 1 ∆t dω α † [ω]α [ω] = dωh (ω)S (i) [ω] ≡ h |S (i) , (14)
where the integration window was approximately extended to R in the first equality,

h (ω) = ∆t -1 |χ [ω]| 2
is the optical population susceptibility and S (i) [ω] = |s (i) [ω]| 2 the energy spectral density of the ith pulse.

From such a population measurement, vector-valued predictions on any ith input energy spectral density S take the form of those of a kernel machine as in Eq. ( 4):

f (S) = B T n(S), (15)
with components:

fn (S) = B n h |S , (16)
where B is a matrix of parameters, S plays the role of the input, and the feature index m in Eq. ( 4) is replaced by that of the optical modes . For optical relaxation times (2πκ -1 ) of the order of tens of picoseconds, such predictions can be realized for ultrashort RF pulses at a throughput above the tens of gigahertz.

Let us check that the corresponding dual-picture kernel is indeed of the form of Eq. ( 12) under some general assumptions on the optical-mode density spectrum of the photonic system. Our considered kernel takes the form:

K(S, S ) = nT (S) n(S ) ≡ S|h h |S = dωdω S[ω]K(ω, ω )S [ω ], (17)
with K(ω, ω ) := h (ω)h (ω ) an integral kernel [START_REF] |ω | may be thought of as components of a same operator K = |h h |. Indeed, one has K(S[END_REF]. Provided all normal modes share the same susceptibility, h (ω) = h(ωω ), and for a continuous optical-mode density spectrum ρ(ω), this amounts to:

K(ω, ω ) = dΩρ(Ω)h(ω -Ω)h(ω -Ω). (18)
By further assuming a smooth spectral density of width larger than the normal-mode linewidth, one has:

K(ω, ω ) ρ ω+ω 2 dth(t)h(-t)e -i|ω-ω |t . ( 19
)
The integral kernel is thus completely determined by the observable susceptibility, here the mode population via

h(ω) = |χ[ω]| 2 /∆t.
For instance, the susceptibility of a linear optical mode χ [ω] = 1/(-i(ωω ) + κ/2) translates into a Lorentzian:

K(ω, ω ) ∝ 1 (ω -ω ) 2 + κ 2 , (20)
that was shown to perform better on some tasks than the more common radial basis function [124]. This integral kernel indeed bears the desired form introduced in Eq. ( 12), with a spectral tolerance δω ∼ κ given by the mode's linewidth.

B. Inspecting the feature space

One can access to the feature space associated to such a photonic kernel machine as well. Indeed, the above kernel admits an eigendecomposition analogous to Eq. ( 7):

K(S, S ) = S|h h |S = γ S|φ φ |S , (21)
with empirical eigenfunctionals and feature maps as given by Eq. ( 8):

φ (ω) = 1 √ N γ u T h(ω), ψ (ω) = γ φ (ω), (22)
where u and the empirical eigenvalues γ are obtained following the previous section as the th eigenvector and eigenvalue of the matrix k/N = H T H/N , with now

H i = h |S (i) = n(i) .
Strikingly, this kernel matrix, and thus the above feature maps, can be directly constructed in experiments from the optical-population measurements over the training set as

1 N k = 1 N i n(i) n(i) . (23)
Now that the feature map is identified, let us examine how the geometric picture introduced in Sec. II A translates to the photonic kernel case. In particular, we consider a regression problem where the targeted quantity is a function of the angular frequency y[ω], as is the case, for instance, when trying to estimate the spectrum of a pulse (y[ω] = S[ω]). This scenario can be dealt with by frequency binning the target as y n = y[ω n ] and trying to best fit it using an estimator of the form ŷ(i) [ω n ] := fn (S (i) ), for any pulse i. Then, it follows from Eqs. ( 16) and [START_REF] Jaeger | The "echo state" approach to analysing and training recurrent neural networks[END_REF], that predictions can be rewritten in the form of Eq. ( 1), as an expansion over feature-space coordinates:

ŷ(i) [ω n ] = ŵ [ω n ] ψ |S (i) (24)
where

ŵ [ω n ] = √ N [u T B] n , (25)
with u = [u 1 , . . . , u M ].
Here, the quantities S ≡ ψ |S play the role of feature-space coordinates [x m in Eq. ( 1)] and can be interpreted as reciprocal-space components of the input S[ω] with respect to the basis { ψ } . The quantities ŵ [ω n ] correspond to the parameters of the featurespace hyperplane [w m in Eq. ( 1)] and can be interpreted as a set of "functions" {ω → ŵ [ω]} learned during the optimization procedure. These may be evaluated explicitly at the chosen frequency bins {ω n } n from the measurement data and the trained parameters B by making use of the above expression. This picture is schematically illustrated in Fig. 3, in analogy with Fig. 1(a), on the spectral-analysis task.

Let us note that by setting the trainable parameters to ŵ [ω n ] = ψ (ω n )/γ n in the absence of any optimization step, one directly has that ŷ(i

) [ω n ] = S (i) [ω n ].
This means that a photonic kernel machine is, at least, able to reproduce the energy spectral density of a pulse s(t) from a single-shot intensity measurement, provided the energy spectral density belongs to the linear span of { ψ } , as will be indeed verified in Sec. V. The action of the ridge regularization of Eq. ( 6) on such a spectral-analysis task can be analytically uncovered. Indeed, upon choosing mean square error as error function,

V (y, ŷ) = n |y[ω n ] -ŷ[ω n ]| 2 ,
the functions learned by the model can be shown to be given by

ŵ [ω] λ = 1 1 + λ/γ × ŵ [ω] λ=0 , (26)
and are thus filtered by the photonic kernel machine in accordance with the inverse of the variance of their associated features over the training set:

γ = E p S| ψ -E p[ S| ψ ] 2 . ( 27
)

IV. PHYSICAL IMPLEMENTATION

We apply the above theory to the ultrafast processing of radio-frequency pulsed signals with a photonic lattice. The task consists in analyzing a set of baseband radiofrequency signals {s (i) (t)} over a bandwidth ∆ω around some reference angular frequency ω 0 by extracting some associated quantity of interest y, such as the energy spectrum or the peak frequency. To this aim, let us introduce an adapted physical implementation of a photoniclattice-based kernel machine. The setup is schematically illustrated in Fig. 4. It involves four successive elements.

Input signals enter the first unit of the system via an electro-optic modulator through frequency mixing with an optical carrier, c(t) = c 0 exp(-iω p t), resulting in a modulated signal of the form F (i) (t) = s (i) (t)e -iωpt whose central angular frequency ω p + ω 0 may be shifted to accommodate the processing of signals in very different bands. The angular frequency of the local oscillator is set to ω p = ωω 0 , where ω denotes the central angular frequency of the lattice's normal modes, to maximize

LO ω p Noise s(t) |s[ω]| 2 ω 0 ∆ω v L×L v 2 v 1 n1 n2 nL×L 1 B = n = ŷ FIG. 4
. Schematic representation of a trained photonic-lattice-based kernel machine estimating the energy spectral density of a noisy pulse from single-shot local intensity measurements. The processing mechanism is the following: (i) A pulsed radiofrequency signal s (i) (t) is sent to the device in the presence of white noise; (ii) the noisy signal is then transferred to an optical carrier of angular frequency ωp thanks to an electro-optic modulator ( ) and (iii) coherently injected into a lattice of L × L linear cavities with random parameters. (iv) The resulting optical populations are then measured through the light radiated by the cavities and, finally, (v) the spectrum of the original noiseless signal is reconstructed by linearly combining the measured quantities thanks to the fixed matrix B, obtained from the previous training protocol.

the response of the system. The now optical modulated signal is then routed to the photonic lattice, entering the cavities as a coherent drive.

The second unit consists of a L × L quadratic photonic lattice whose cavities are mutually coupled via near-field nearest-neighbor interactions. These cavities are coupled to the external modulated drive with some arbitrarily spatially-dependent weight v . In a frame rotating at the drive frequency, the dynamics of such a lattice is described by the following set of quantum Langevin equations [START_REF] Gardiner | Input and output in damped quantum systems: Quantum stochastic differential equations and the master equation[END_REF][START_REF] Gardiner | Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics[END_REF]:

∂ t â (t) = i(∆ + ω 0 ) -κ /2 â (t) + i ∈v( ) J , âm (t) + iv s(t) + κ /2 âin (t), ( 28
)
where â is the photon annihilation operator for site , ∆ = ω pω is the detuning of the local oscillator, κ the optical relaxation rate, J , the linear coupling rate between cavities and , s(t) the broadband driving signal, and v the local weight of the coupling of the th cavity to the external drive. Operators âin (t) account for quantum noise and satisfy âin (t) = âin † (t)â in (t ) = 0 and [â in (t), âin † (t )] = δ , δ(tt ). A high degree of control of the parameters is not required. As a matter of fact, in the following numerical simulations, while angular frequencies ω will be uniformly set for simplicity, with ∆ = -ω 0 , the remaining parameters will be random variables. In particular, J m, will be uniformly drawn in the interval [0, J max ], κ normally distributed around κ = zJ max /40 (z = 4) with a standard deviation of 10%, ∆ω = zJ max and v will be set to a normalized random real vector.

The third unit consists of a sensor that measures the intensity of the light radiated by the optical cavities as they are externally driven by an input signal. These populations are measured by time-averaging over some detector integration time ∆t and collected into a vector

n(i) = [n (i) 1 , . . . , n(i) L×L , 1]
T , where a unit entry is added to get a supplementary trainable parameter acting as a bias. In what follows, two alternative measurement settings will be considered:

(i) Measure of the local populations:

n(i) = 1 ∆t
dt â † (t)â (t) . This corresponds, for instance, to the experimental situation where the intensity that is vertically emitted by the cavities is measured by a camera facing the lattice.

(ii) Measure of the normal-mode populations: n(i) = 1 ∆t dt α † (t)α (t) . This corresponds, for example, to the experimental situation where the field leaking from the cavities is collected by an evanescently coupled waveguide and frequency-demultiplexed into L × L frequency channels coupled to photodetectors.

Finally, a fourth unit performs a linear combination of the measured populations by acting with a n × (L 2 + 1) matrix of parameters B, to be optimized by training. The output of this last unit is thus of the form ŷ(i) = B T n(i) , as in Eq. [START_REF] Lin | All-optical machine learning using diffractive deep neural networks[END_REF].

A. Training the lattice-based photonic kernel

The training starts from a training and a testing sets composed of N train and N test signals, respectively, each of the form {(s (i) , y (i) )} i , where s (i) denotes the ith signal and y (i) a set of known associated features we want our system to learn how to estimate. In regression tasks, these labels y (i) may take arbitrary values, whereas in binary classification tasks y (i) = ±1 depending on whether the associated input s (i) belongs to a targeted class or not.

The trial function of the untrained model is initially given by f (s

(i) ) = B T n(i) [125].
The training is carried out as follows:

(i) Construct the matrix Y. From the known features {y (i) } i associated to the samples of the training set, Y is first computed as Y im = y (i) m . (ii) Obtain the matrix H. Each input signal s (i) of the training set is fed into the device at the electro-optic modulator and the resulting cavity populations n(i) are measured by the sensor. All these measured populations are stored in a matrix

H im = n(i) m .
Note that this evaluation is only to be performed once.

(iii) Determine the optimal parameters B. For a chosen value of the regularization hyperparameter λ, one obtains the optimal parameters by minimizing the cost function as B = arg min B J(B), with

J(B) = Ntrain i=1 V (y (i) , f (s (i) )) + λ 2 B 2 2
, where f (s (i) ) = B T n(i) , with n(i) as obtained in the previous step. The choice of loss function V depends on the task. For regression tasks, a popular choice is simply V (y, ŷ) = yŷ 2 , which corresponds to a least-square problem; this directly yields B = (H T H + λ1) -1 H T Y analytically, with Y and H as computed at steps (i) and (ii), respectively. In binary classification (y = ±1), it is customary to choose a margin-maximizing loss functions [START_REF] Rosset | Margin maximizing loss functions[END_REF]. Popular choices [START_REF] Hastie | The Elements of Statistical Learning[END_REF] are the hinge loss, V (y, f (s)) = max(0, 1y f (s)), which makes f directly approximate the class label f (s) = y, or the binomial deviance ln(1 + e -y f (s) ), which makes f instead approximate f (s) = ln(P(y = +1|s)/P(y = -1|s)). Upon choosing any of these loss functions, the problem is convex and can be solved either analytically, for the square error, or by means of a convex optimization solver, using standard iterative methods.

(iv) Evaluate the accuracy of the model. Once the model is trained, for any new input signal s fed into the system, features are estimated according to f (s) = BT n, from the measurement of the resulting populations n. Its accuracy can then be benchmarked on the testing set by comparing the predictions ŷ(j) = f (s (j) ) against the known features y (j) , for every input signal s (j) in the testing set, to which the model was yet never exposed. This is done through a metric that may differ from the cost function, in particular if regularization was employed.

B. Benchmarking the photonic-lattice kernel machine

In the following, two different approaches will be employed to benchmark the above-defined physical implementation of a photonic kernel.

In order to evaluate the ability of the model to estimate the spectrum of pulsed signals, we will first use the modulus square of the fast Fourier transform (FFT) of the input signals as a reference of energy spectral density. We will assume an ideal sampling of the pulse over a centered window of time length ∆T = 5 × 2π/κ at a sampling rate of f s = 200 × κ/2π. For cavities with 2π/κ ∼ 10 ps, this corresponds to a sampling rate of f s = 10 THz. Note that this is more than three orders of magnitude beyond the state of the art [START_REF]See for instance Texas Instruments' ADC12DL3200 (6.4 GSPS) or Analog Devices[END_REF]. Yet in the following the photonic kernel will be shown to outperform this rather fictional ideal device.

We will also use a reservoir computing approach based on a nonlinear polaritonic lattice as first introduced in Ref. [START_REF] Opala | Neuromorphic Computing in Ginzburg-Landau Polariton-Lattice Systems[END_REF]. This reservoir was numerically proven successful in image and speech recognition tasks and chaotic time-series forecasting, and was recently experimentally tested on an optical character recognition task [START_REF] Ballarini | Polaritonic Neuromorphic Computing Outperforms Linear Classifiers[END_REF]. This reservoir is modelled by the following discrete complex Ginzburg-Landau equation:

∂ t α (t) = + i(∆ + ω 0 ) + γ -(Γ + ig)|α (t)| 2 α (t) + i ∈v( ) J , α m (t) + iv s(t), (29)
where γ = Pκ/2 is a gain coefficient that accounts for the single-body decay rate κ/2 and the magnitude of the external pumping P , that brings the system close to instability. Here, Γ and g respectively account for twobody dissipation and interaction processes. The input signals are coherently injected in the same fashion as for the linear lattice. In the following numerical simulations, we will set ∆ = -ω 0 , J , uniformly drawn in the interval [0, J max ], Γ = zJ max /40 (z = 4), γ/Γ = 8 × 10 -4 /2π, g/Γ = 1.6/2π, ∆ω = zJ max and v to a normalized random real vector. The output of this model depends on the amplitude of the input. In what follows, the input energy will be set to 50Γ.

The training of this model is realized exactly as for the photonic kernel machine, from local population measurements n(i) FIG. 6. Relative absolute error ∆E/E0 on the testing set for three values of the signal-to-noise ratio (SNR) for the fast Fourier transform (horizontal dashed line), the complex Ginzburg-Landau reservoir computing model (dash-dotted) and the photonic kernel machine where either the local-mode (plain) or the normal-mode (dashed) populations are measured. Data is averaged over 5 realizations of the reservoir, error bars correspond to the intervals between the lowest and highest errors over these realizations. For each realization: Ntrain = 7000 and Ntest = 3000.

V. APPLICATIONS

A. Pulse spectral analysis

As a first illustration of the learning capabilities of the above-described photonic kernel machine, let us consider the extraction of the noiseless energy spectral density of an ultrashort pulsed radio-frequency signal embedded into a noisy background. Given an input noisy pulsed sig-nal s(t) + ξ(t), this task consists in giving the best possible estimation ŷn of its frequency-binned energy spectral density y n = |s[ω n ]| 2 regardless of the noisy background ξ(t).

In order to train and study the performance of the model, we generate training and testing sets of N train and N test pulses, respectively. To do so, we first generate known noiseless random spectra of same total energy E 0 = dωS ss [ω] by cubic B-spline interpolation of a set of N b signal bins sn := s[ω n ], with frequency bins {ω n } n equally spaced within the band [ω 0 -∆ω/2, ω 0 + ∆ω/2] and randomly sampled from a Boltzmann distribution p[s n ] = (1/Z)e -βV [sn] parametrized by the following potential:

V [s n ] = a N b N b -1 n=1 |s n+1 -sn | 2 + b max n |s n | 2 -S peak 2 + c std n |s n | 2 -S 2 , (30)
with si = 0, i = 1, N b , S peak = 8 S, S = E 0 /∆ω. Here, a acts as stiffness parameter whereas the b-term favors peaked spectra. Finally, the c-term prevents the sampled spectra from sharing too similar shapes by favoring the presence of secondary peaks. The parameters used throughout this section are N b = 20, βa = βc = 100 and βb = 50. The energy spectral densities of the obtained random noiseless spectra are shown in Fig. 5 (a), where the average energy distribution is compared to S and a few typical examples are plotted, exhibiting various peaks with different heights. The choice of bandwidth (∆ω = zJ max ) ensures that all the power of the signal to be analyzed can be sensed by the photonic lattice, as shown in Fig. 5 (b), where the average optical spectrum of a large random photonic lattice is plotted. The random spectra are then Fourier-transformed to the time domain and white noise ξ(t) is added to match some signal-tonoise ratio SNR, here defined as the relative contribution of the noise to the total energy in the analyzed band, i.e.

SNR = E -1 0 ω0+∆ω/2 ω0-∆ω/2 dω|ξ[ω]| 2 .
We then simulate the response of the photonic lattice to each of the driving input signals s (i) (t) of the training set by numerically integrating the coupled dynamical equations [START_REF] Rodan | Minimum Complexity Echo State Network[END_REF]. For each signal, either local or normalmode populations are then measured yielding a set of time-averaged populations n(i) . The weights B are then optimized over the training data so as to minimize the square error between the predictions of the photonic kernel machine ŷ and the known spectra of the noiseless pulses y. Upon fine-tuning of the hyperparameter λ by 10-fold cross validation, the optimal weights B are obtained analytically as explained above. The error of the trained model is then evaluated on the test set. We here quantify this error by the relative absolute error on the energy ∆E/E 0 , where ∆E represents the energy area between the estimated and the actual energy spectral density curves of the original noiseless spectra. In Fig. 6, we use this metric to benchmark our model against the ideal FFT of the noisy input signal and a the nonlinear polariton-based reservoir. For any of the chosen signalto-noise ratio values and lattice sizes, the photonic kernel outperforms the other two approaches. Interestingly, the maximum performance is already reached with lattices as little as 10 × 10, with an error in the reproduced spectrum about three times smaller than that obtained by the ideal FFT procedure. In addition, Fig. 7 shows the error frequency over the testing set for the three values of the signal-to-noise ratio and a 20 × 20 photonic lattice. Note that, for any amount of noise, the worst-case predictions remain more accurate than the ideal FFT. In order to understand how the device learns from training examples, let us make use of the theoretical concepts introduced above.

As expressed in Eq. ( 24), regression by the implemented photonic kernel machine is performed by a linear expansion ŷ(i) [ω n ] = ŵ [ω n ] ψ |S (i) of the feature-space components ψ |S (i) over a set of learned functions ŵ . The eigenvalues γ and eigenfunctions ψ / √ γ of the photonic kernel do not depend on the optimization procedure and can be given empirical estimations from the measured populations n(i) during the training process, as described above. In Fig. 8 (a), we show the convergence of the empirical eigenvalues as the amount of training examples is increased, that is already reached for roughly N train = 1000 with our training protocol. Figs. 8 (b) and (c) show the leading empirical feature maps as well as the corresponding learned eigenfunctions. On this task, one observes that the optimization procedure leads to a set of learned functions that correspond to filtered analogues of the empirical eigenfunctions of the kernel. One sees from Fig. 8 (b) that the model builds a Fourier-sine expansion with a spectral resolution cut-off at ∼ ∆ω/2πL 2 . While in principle the optimization procedure may be sensitive to any feature-space component within this bound, the ridge regularization introduced during the optimization induces a soft cut-off for those whose associated eigenvalues have magnitudes lower than λ. This reduction to only the most statistically relevant components of the functional basis prevents the model from overfitting the training set, which would undermine its generalization capacity. In Fig. 8 (c), the effect of the regularization is clearly visible on the highest-order represented learned function, that is completely filtered out. Hence, regularization here manifests itself as a low-pass filter on the learned decomposition. From this figure, the interpretation of the the photonic kernel machine regression mechanism becomes very clear: (i) the training set determines some optimal Fourier-like decomposition of the spectra; (ii) the regularization truncates the basis of such a decomposition to its most statistically relevant components, preventing the device from overfitting the noise; and, finally, (iii) the optimization finds a set of smooth filtered base functions on which to expand back those components to best reproduce the features. This results in the reconstruction of the noiseless spectrum, from which the noisy background uncorrelated with the training set is regularized out.

B. Pulse shape recognition

We now show the performance of the same photonic kernel machine trained on a noisy-pulse shape classification problem.

The task consists in determining whether the envelope of an input pulse is Gaussian (y = 1) or Lorentzian (y = -1) from the measurement of the photonic populations of the lattice in the presence of noise at the input port of the system. We prepare a set of pulses with central angular frequencies uniformly drawn at random in the band of interest [ω 0 -∆ω/2, ω 0 + ∆ω/2] and full widths at half maximum (FWHM) normally distributed around FWHM = ∆ω/20 with a relative standard deviation of 10%. The induced populations are then time-integrated over some detection time ∆t yielding a vector of intensities n that are finally linearly combined in such a way that the output of the photonic kernel machine is now a scalar of the form f (S) = β T n. The optimization process is realized by minimizing the hinge loss V (y (i) , f (S (i) )) = max(0, 1y (i) f (S (i) )). This is here achieved by means of a convex optimization solver [START_REF] Udell | Convex optimization in Julia[END_REF] for λ → 0 + .

The output of the ith pulse may be equivalently rewritten in terms of components of the feature map as f (S (i) ) = ŵT S(i) + b, with S(i) = ψ |S (i) . It is worth noting that f (S) = ŵ × d( S), that is, the prediction for any incoming pulse s is proportional to the signed distance d( S) of its energy spectral density's feature-space coordinates to the plane ( ŵ, b). As discussed above, this plane defines the feature-space decision boundary of the model. Input pulses are then classified into either of the two classes depending on whether their feature-space coordinates S(i) fall on either sides of this plane. Hence, predictions are of the form ŷ(i) = sign( f (S (i) )).

The classification process is illustrated in Fig. 9 for a 20 × 20 photonic kernel machine and intermediate noise strength (SNR = 20). In panel (a), one observes that features corresponding to either classes indeed cluster at either of the sides of the discriminating hyperplane independently from the value of the FWHM. As it becomes clear in panel (b), predictions become less accurate as the distance from the separating hyperplane becomes smaller, thereby giving an estimation of the likelihood of the prediction. This can be made more quantitative by calibrating the probability of the classifier. This probability is shown in panel (c) as given by P(s

(i) = "Gaussian"| f (S (i) )) = [1+exp(-A f (S (i) ) + B)] -1
, where the calibration parameters A and B were determined by Platt scaling [START_REF] Platt | Probabilistic Outputs for Support Vector Machines and Comparisons to Regularized Likelihood Methods[END_REF].

The performance of a 20×20 trained photonic kernel machine is shown via its receiver operating characteristic (ROC) curve in Fig. 10 for increasing noise strengths and the two population measurement scenarios. This displays the sensitivity (true positive rate) as one allows the specificity of the model to drop (higher false positive rates) by playing on some external bias added to the trained model, the best trade-off being found at the top left corner for no external bias. The ROC curve of an unbiased classifier that affects pulses randomly to either shape class is plotted as well for comparison. The high sensitivity and specificity values of the trained classifier are given in the right panel of Fig. 10 for both population-measurement protocols.

C. Frequency tracking

Above, only the case of pulsed input signals was investigated. In the following, we illustrate the performance of the photonic kernel machine presented above on the analysis of continuous radio-frequency signals by considering frequency estimation of noisy sinusoidal radio-frequency signals. This is of great relevance for instance in the context of short-timescale force sensing with optomechanical devices [START_REF] Allain | Optomechanical resonating probe for very high frequency sensing of atomic forces[END_REF][START_REF] Guha | Force Sensing with an Optomechanical Self-Oscillator[END_REF]. In such applications, forces are read out by frequency tracking of the harmonic radiofrequency modulation imprinted on an optical signal by a resonating mechanical probe.

To do so, we first generate a first set of training baseband signals consisting of complex exponentials with random initial phases and angular frequencies ω uniformly drawn between ω 0 -∆ω/2 and ω 0 +∆ω/2, to which white noise is added so as to match some signal-to-noise ratio SNR, here defined as the ratio between the average power of the sinusoidal baseband signal and that of the noise. We then measure the steady-state populations of the cavity resulting from the driving of the coupled cavities by the modulated signals, here after a time τ d = 10/κ, and use this vector of populations to make frequency predictions of the form f (s) = β T n. The vector β is then optimized so as to minimize the mean squared error between the estimated ω(i) = f (s (i) ) and the actual ω a ridge-regularization hyperparameter determined by 10fold cross validation. The performance of the trained photonic kernel machine is finally evaluated on a testing set composed of new random complex exponentials. The achieved average resolution is shown in Fig. 11 as a function of the lattice size for increasing values of the SNR, revealing the performance of the above-described photonic kernel machine on the spectral analysis of continuous radio-frequency signals. For cavities with 2π/κ ∼ 10 ps, the waiting time would be as little as τ d ∼ 10 ps.

VI. CONCLUSION

In this article, we presented photonic kernel machines, a new framework for optical ultrafast spectral analysis of noisy radio-frequency signals that translates kernel methods from support-vector machines to photonic hardware. Such devices realize regression or classification tasks on high-dimensional data with throughputs above the gigahertz by utilizing the optical response of a set of optical modes to input analog signals as a measure of similarity.

We first gave a theoretical description of photonic kernel machines under very general assumptions. We analytically investigated the similarity kernel built-in in such devices and were able to express it explicitly from the susceptibility of the measured observables. Furthermore, we explored the feature maps associated to photonic kernel machines and found that their expressions could be experimentally determined from population measurements and the knowledge of the single-mode susceptibility.

We then studied a model describing a physical implementation consisting of a lattice of coupled linear optical cavities. We numerically demonstrated its capabilities on various regression and classification tasks, comprising the analysis of both pulsed and continuous radio-frequency signals. In particular, the proposed setup proved efficient in predicting the spectrum of picosecond pulses with nontrivial spectral structure from single-shot intensity measurements, being able to predict spectra with higher fidelity than the FFT of the noisy input signal. Moreover, it was shown to be able to discriminate pulses with distinct shapes as well as to estimate the angular frequency of continuous harmonic input signals. We showed that, by adding noise at the input of the device during the training protocol, the spurious effect of background noise on the predictive performance of the device could be successfully mitigated. On the spectrum estimation task, we could extract the actual feature maps associated to the simulated kernel machines as well as the basis of learned functions the photonic kernel machine composes its predictions from. This allowed us to interpret the photonic kernel machine regression mechanism and revealed the ability of the system to filter out the uncorrelated background noise.

We believe that such devices, capable of analyzing above one million radio-frequency signals per second, may found applications in a broad variety of domains beyond spectroscopy. In the field of radio-frequency sensing, it could be used in pulse-Doppler radar systems as a way to optically analyze the reflected signals. In this way, the rate of emission could be increased by orders of magnitude, by relaxing the limiting dependence on the sampling rate of analog-to-digital converters. In telecommunications, these devices could be used, for instance, as a decoding means in FSK protocols involving high modulation rates in noisy environments. The fast frequencytracking ability of the described device could be exploited in the context of short-timescale force sensing with optomechanical resonators. Finally, photonic kernel machines could be integrated into more conventional machine learning pipelines as a means of extracting nontrivial digital features from analog signals, acting as both an analog-to-digital converter and a preprocessing stage.

ψ

  FIG. 1. Schematic representation of the embedding from input to feature space by a shallow model. (a) Example of regression: the quadratic function f (x) = 3x -x 2 + 4 is linearly fitted after the embedding of Eq. (1), with ψ1(x) = x 2 and ψ2(x) = x. Note that, in feature space, all data points lie on a same plane. (b) Example of binary classification: triangles and dots become linearly separable in feature space after the embedding of Eq. (1), with ψ1(x) = x1, ψ2(x) = x2 and ψ3(x) = x 2 1 + x 2 2 . The resulting input-space decision boundary f (x) = 0 is represented by a dashed line. Origins are shifted to improve the legibility.

∝ 1 ∆tFIG. 5 .

 15 FIG. 5. (a)Energy spectral density averaged over 10 5 realizations of the noiseless pulses, as given by Eq. (30), as well as that of a few individual realizations. (b) Mean and interval between 5% and 95% quantiles for the optical spectra of 100 realizations of a 30 × 30 random photonic lattice.

  ∝ â † â ) PKM (n ∝ α † α )

10 FIG. 7 .FIG. 8 .

 1078 FIG. 7.Frequency of the error score over the testing set for three values of the signal-to-noise ratio averaged over 5 realizations of a 20×20 photonic kernel machine. Ntest = 3000 for each realization. The error score of the FFT is shown as dashed lines for comparison and falls systematically above the upper bound of that of the photonic kernel machine, in spite of being slower and requiring an ideal sampling rate significantly beyond the state of the art.

FIG. 9 .

 9 FIG. 9. (a) Signed distance from the testing pulses features to the learned hyperplane (decision boundary) as a function of their FWHM for a 20 × 20 photonic kernel machine under normal-mode measurement of the photonic populations. Dots correspond to Gaussian pulses, squares to Lorentzian ones. Pulses falling above the decision boundary (dashed line) are categorized as Gaussian by the classifier. Parameters: SNR = 20, Ntrain = 14000 and Ntest = 6000. (b) Frequency of each class as a function of the signed distance to the discriminating hyperplane. (c) Probability that a pulse be Gaussian or Lorentzian at any given value of its associated signed distance to the discriminating hyperplane.

  (i) angular frequencies of all signals in the training set, with

FIG. 10 .

 10 FIG. 10. Receiver operating characteristic (ROC) curve for a 20 × 20 photonic kernel machine on the pulse shape classification task for three noise strengths and two measurement scenarios. The closer to the top left corner, the better. The numerical values of their associated sensitivities and specificities are compiled on the right tables. Parameters: Ntrain = 14000 and Ntest = 6000.

FIG. 11 .

 11 FIG. 11. Resolution of the photonic kernel machine on the testing set for either local-mode (plain) or normal-mode (dashed) population measurements for increasing lattice sizes L and three values of the signal-to-noise ratio. The performance of the nonlinear polariton-based reservoir-computing scheme (dash-dotted) is shown for comparison. Data is averaged over 5 realizations of the reservoir, error bars correspond to the intervals between the lowest and highest resolutions over these realizations. Parameters: Ntrain = 7000 and Ntest = 3000.

ACKNOWLEDGMENTS

This work was supported by the ERC Consolidator grant NOMLI No. 770933, by ANR via the project UNIQ, and the FET flagship project PhoQuS (grant agreement ID No. 820392).

Appendix A: Feature-space inspection Kernel machines' feature space embedding may be straightforwardly examined. Indeed, for a symmetric positive semidefinite kernel K, Mercer's theorem ensures that it admits an eigendecomposition of the form:

, where the measure is set to the probability density function of inputs p. The feature map of Eqs. ( 1) and ( 2) can then be easily shown to be given by ψ m (x) = √ γ m φ m (x). We shall now identify the elements of this eigendecompositon.

To identify its eigenvalues and eigenvectors, one has to solve in principle for

Yet these may be estimated from the N known training samples by replacing the input distribution p(x) with the empirical one p(x), such that dxp

. This yields a tractable discrete eigenvalue problem:

in terms of empirical eigenvalues γm and empirical eigenfunctions φm . It follows that the empirical eigenvalues correspond to the non-zero eigenvalues of either of the kernel matrices K = HH T and k = H T H (let us recall

. Indeed, because of their Gram matrix structure, their eigendecompositions,

with U im = φm (x (i) )/ √ N , share the same non-zero eigenvalues and can be related through the following simple algebraic identity:

following K's singular-value decomposition. In the physical model of the main text, one has typically M N and it becomes more suitable to work with the M × M matrix k.

Similarly, the empirical eigenfunctions can be determined from Eq. (A3). By making use of the property of Eq. (A5), these finally read:

From Eq. (A1) and Eq. (A6), the kernel can be expanded into a series of empirical eigen feature maps K(x, x ) = M m=1 ψm (x) ψm (x ), with orthogonal empirical feature maps simply given by ψm (x) = γm φm (x).

(A7) Thus, the internal functional representations of the model are independent from the optimization process, completely determined by the statistics of the training samples and can be estimated by diagonalizing the matrix k.

The predictions of the model in the abstract feature-space picture, f (x) = w T x ≡ w T ψ(x), take the form of a linear combination of features whose weights, determined through training, can be geometrically interpreted as the parameters of a hyperplane. As a function of the primal parameters, this learned hyperplane is characterized by

Therefore, provided a set of generating functions {h m } m , one is now able to access to a complete understanding of the model and its feature-space representations.