Photonic reservoir computer using speckle in multimode waveguide ring resonators

Matthew N. Ashner; Uttam Paudel; Marta Luengo-Kovac; Jacob Pilawa; George C. Valley

doi:10.1364/OE.425062

1. Introduction

Photonic machine learning promises an energy-efficient platform for real-time processing of wide-bandwidth signals and images [1–6], such as classification, prediction, and denoising. Among various classes of photonic machine learning, photonic reservoir computing exploits physical phenomena to generate large random matrices and perform matrix multiplications passively in an analog domain. Recently, several researchers [7–10] showed how to use multimode propagation at the speed of light for reservoir computing without the need to pass information back and forth between memory and processing units [11,12].

Broadly speaking, there are two photonic reservoir computing architectures [2]: one uses a single-mode ring with virtual neurons that are read out serially in time [13–19], and the other uses speckle in a multimode waveguide [8,10,20], scattering medium in free-space [6,7,21,22], or an optical cavity [23]. Photonic implementations of the single-mode ring resonator reservoir architecture with virtual neurons have been demonstrated to perform a few orders of magnitude better than electronic systems [17]. However, the readout rate scales as the input signal bandwidth times the number of neurons. This makes chip-scale integration of such a system impractical for wideband signals since this would require an analog-to-digital converter (ADC) with an input bandwidth greater than tens of GHz. The fastest ADCs in the current literature have bandwidths just above 10 GHz but have signal-to-noise-plus-distortion ratios of 30-35 dB [24–26]. The 5-6 effective bits of these ADCs are insufficient for machine learning tasks, so current demonstrations of these systems rely on digital oscilloscopes [14]. In addition, because the length of the delay line scales as the number of neurons, building a chip-scale system with more than about 100 neurons is also impractical.

The second photonic architecture processes an array of neurons and addresses this problem by using parallelized ADCs. In this architecture, the number of neurons would equal the number of ADCs operating at the speed of the signal of interest. One of the most attractive types of array processing uses laser speckle from rough surfaces or multimode waveguides. However, one of the key bottlenecks in implementations of a speckle-based RC is the requirement for feedback. In our prior work [8,10] and work by Dong et al. and Rafayelyan [7,21], a spatial light modulator (SLM) was used with electrical/digital feedback. Such implementations are inherently inefficient as they require one costly optical-to-electrical-to-optical (OEO) conversion per pixel per time step. They also require as many SLM elements as the number of neurons, making chip-scale integration challenging due to the difficulty of fabricating and controlling hundreds or thousands of modulators on one chip and the power required to operate them at high speed.

To circumvent these limitations, we propose and demonstrate through simulations a photonic reservoir computer (RC) architecture that uses a multimode ring resonator with a built-in nonlinear element. The neuron vector or matrix is fed back optically avoiding the OEO conversion at each time step. This hybrid architecture combines the all-optical feedback of the single-mode ring resonator architectures with the parallelized neuron readout of the speckle architectures into a design exhibiting the advantages of both. This architecture reduces the number of required electro-optic modulators (EOMs) from the number of neurons to the dimensions of the input signal or image. This represents a large reduction in EOMs since the number of neurons must be much larger than the dimensions of the input signal to achieve the RCs dimensionality expansion.

2. Ring resonator reservoir computer simulation

A generic reservoir computer is comprised of 3 layers (Fig. 1(a)): an input layer containing vectors of input data, written u_n for the n^th timestep, a reservoir layer consisting of neurons with random, fixed connections, and an output layer yielding one or more output predictions [27,28]. The reservoir connections are defined by a random, square matrix, W, and the input weights are defined by another fixed, random matrix, W_in. The states of the reservoir at the n+1^th time step, x_n₊₁, can then be defined recursively as

(1)$${{\mathbf x}_{n + 1}} = {f_{NL}}({{\mathbf W}{{\mathbf x}_n} + {{\mathbf W}_{in}}{{\mathbf u}_n}} )$$

where f_NL is a nonlinear activation function such as a hyperbolic tangent or a reLU. The predictions yielded by the output layer are then defined as

(2)$${{\mathbf y}_n} = {{\mathbf W}_{out}}{f_{NL,out}}({{{\mathbf x}_n}} )$$

where W_out is a matrix containing the output weights trained by some optimization routine, and f_NL,out is an optional output nonlinear activation function. The one change we make to the typical presentation of these equations is to reverse the nesting of the output nonlinearity and multiplication by W_out to reflect the physics of the photonic system. This change does not affect the ESN’s predictive capabilities as demonstrated in the following sections. In this work, we use ridge regression for training W_out in all the examples described.

Fig. 1. a) Topology of an echo state network. The connections between the input channels and the reservoir neurons and the connections within the reservoir are fixed and random. The connections to the output layer are trained. b) Physical architecture of an optical ring resonator RC with the mathematical operations corresponding to each physical component.

Download Full Size | PDF

The choice of W and W_in along with the functional forms of the reservoir and output layer nonlinear activation functions fully define a generic RC. We simulate an optical RC by calculating these matrices and functions in terms of properties of a waveguide ring resonator and performing a generic RC calculation. Beyond this, there are two slight differences between the optical and conventional RCs. First, W, W_in, and the reservoir state vectors x_n contain complex numbers representing the amplitude and phase of the optical field rather than the strictly real numbers typically used. Second, because the input data are injected using the intensity of an input laser, the optical RC cannot accept negative numbers for input data, and the inputs must be recast into all positive numbers. Figure 1(b) shows the architecture of the optical ring resonator RC and the physical elements that correspond to the mathematical elements that define the RC. The training and evaluation of W_out is performed computationally on the image data read out from the camera.

To describe this process, we begin by calculating the W matrix. In the optical ring resonator, each “neuron” is the portion of the waveguide cross-sectional area sampled by one pixel of the camera or photodiode array. Thus, the state vector x_n in Eq. (1) contains the complex field values at each discrete spatial position immediately after application of the nonlinearity as shown in Fig. 1. The action of propagating the light one lap around the waveguide transforms one pseudorandom interference pattern into another with a fixed transfer function, thus replicating the action of the W matrix. In our simulations, the W matrix is a spatial mapping of electric field pattern at the beginning of the resonator loop to the field pattern at the end of the resonator loop that encompasses the physics of propagating the light through the waveguide. The optical field propagating through a waveguide is described in terms of discrete modes, each with its own propagation constant, β_k, and transverse spatial field pattern in the waveguide, ${E_k}(x,y)$,

(3)$$E(x,y,z) = \sum\limits_{k = 1}^N {{a_k}{E_k}(x,y)\exp ({\textrm{ - i}{\beta_k}z} )} .$$

The first step in performing a space-space mapping is to calculate the complex expansion coefficients, a_k, representing the amplitude and phase for each mode such that the sum in Eq. (3) reproduces the projection of the input field onto the propagating modes. The second step is to reconstruct the field after propagation around the ring resonator using Eq. (3), at z = L, the round-trip length of the waveguide.

To represent this process mathematically, we first calculate a matrix M containing the field pattern, ${E_k}(x,y)$, for each mode evaluated on a fixed grid of spatial points. This matrix will have dimensions P×N, where P is the number of spatial points in the grid and N is the number of transverse modes supported by the waveguide, and it transforms a set of expansion coefficients, {a_k}, into an equivalent spatial field pattern. Thus, each element of the M matrix, M_p,k, is the electric field of the k^th mode evaluated at the p^th spatial position. Sections 3 and 4 describe the field pattern calculation for each type of ring resonator. We construct the matrix, W, that represents the composite process of diverting a fraction of the intensity from the ring resonator, calculating the expansion coefficients for each mode from the spatial field pattern, propagating the modes around the ring resonator, and converting back to the spatial domain as

(4)$${\mathbf W} = \sqrt {1 - \eta } {\mathbf{MB}}{{\mathbf M}^ + }.$$

M⁺ represents the Moore-Penrose pseudoinverse of M, B is an N×N diagonal matrix containing the phase factors applied to each mode as it propagates around the resonator, $\exp ({ - \textrm{i}{\beta_n}L} )$, and η is the outcoupling fraction, the fraction of optical power split to the camera at the output. In this construction, W is a P×P matrix that maps each spatial point at the beginning of the resonator to a field profile at the end of the resonator. The construction of W_in is similar.

(5)$${{\mathbf W}_{in}} = {\mathbf{MB}}{{\mathbf M}^ + }{\mathbf C}$$

where C is a matrix that maps the input and bias channels to the spatial points and is determined by the spatial profile of the input signals injected into the waveguide. Thus, C has dimensions P×(U+1) where U is the number of input data channels and the additional 1 corresponds to the optional bias channel.

Because the neuron feedback bypasses the optical detection step, we can no longer rely on the intensity measurement to perform the nonlinear activation. Instead, we must apply the nonlinear activation in the optical domain. In the model, we can choose any activation function that can be reasonably replicated in an optical process. For our simulations, we use the hyperbolic tangent function frequently utilized in RC calculations, which can be approximated with an optical limiter or two-photon absorption medium, for example. Because the neuron states are complex, the exact definition of f_NL needs to be reformulated slightly. Saturable gain media and two-photon absorption stages will apply an amplitude gain or loss that depends on the incident intensity, but the phase accumulated by the electric field will be roughly constant with intensity [29]. Because a uniform phase factor will not affect the speckle pattern formed by interference between the modes, we ignore this unknown phase factor and preserve the phase of the incoming field at each point. We apply a hyperbolic tangent function to the field intensity at each position on the spatial grid to approximate the action of the optical nonlinear activation, so the nonlinear activation function is

(6)$${f_{NL}}(x) = {\left[ {\gamma \tanh \left( {\frac{{{{\left\Vert x \right\Vert}^2}}}{\gamma }} \right)} \right]^{\frac{1}{2}}}\exp [{\textrm{iArg}(x )} ],$$

where γ is the saturation intensity. That said, optical devices can replicate many different non-linear activation functions, and not all of them apply a uniform phase factor [30]. However, the precise functional form of the nonlinear activation often has little effect on RC performance [30], so the hyperbolic tangent is a suitable choice for this model.

Finally, because a camera measures the optical intensity and not the electric field, there is an output nonlinearity that includes the modulus squared operation at a minimum

(7)$${f_{NL,out}}(x) = {{\big \Vert} x {\big \Vert}^2}.$$

Other output nonlinearities carried out in digital computations are possible, but these must be applied to the measured field intensity in the geometries for the ring resonator RC that we have investigated. For example, the simulations of the multimode optical fiber RC shown in Section 4 apply W_out to the measured intensity squared in the trained model.

This completes the general description of an optical ring resonator RC. For each of the waveguide geometries, the matrices M and B are calculated from the transverse spatial modes and propagation constants derived from wave optics models as described in the following sections, and the optical RC is simulated with standard RC calculations using the formulation given in Eqs. (4)–(7).

3. Planar waveguide reservoir computer simulations

For a straight planar multimode waveguide, there is a well-known procedure for calculating the field profiles and propagation constants for each transverse mode [31]. Here we assume a resonator with a sufficiently large radius of curvature as in Fig. 1 that we can ignore the bends. Thus, the optical field at the end of the resonator is derived by propagating the optical field from the beginning of the resonator as if the guide were straight. This approximation enables the use of analytical solutions for the mode field profiles and avoids the need for computationally expensive finite-element simulations. We argue that the bends will increase the mode mixing needed for an RC. If there is an effect, omitting them will underestimate the performance of an RC in an optical resonator as it may give a less dense effective mixing matrices for the neuron vector and the input signal.

The transverse magnetic (TM) modes of a planar waveguide are given by

(8)$$\begin{array}{l} {E_j}({x,z} )= \left[ {\sin ({{h_j}x} )- {{\left( {\frac{n}{{{n_c}}}} \right)}^2}\frac{{{h_j}}}{{\kappa ({{h_j}} )}}\cos ({{h_j}x} )} \right]\exp [{\textrm{i}\beta ({{h_j}} )z} ]\\ \kappa (h )= {\left[ {\beta {{(h )}^2} - {{\left( {\frac{{2\pi n}}{\lambda }} \right)}^2}} \right]^{\frac{1}{2}}}\quad ,\quad \beta (h )= {\left[ {{{\left( {\frac{{2\pi {n_c}}}{\lambda }} \right)}^2} - {h^2}} \right]^{\frac{1}{2}}} \end{array}$$

where x is the transverse coordinate across the guide, z is the coordinate in the direction of propagation, j is the mode index, n and n_c are the refractive indices of the core and cladding respectively. The transverse wavenumbers h_j for the jth mode are determined by solving the standard transcendental equation obtained from the boundary conditions at the edge of the guide,

(9)$$\begin{array}{l} \tan ({hd} )= \frac{{2hd{{({{V^2} - h{d^2}} )}^{\frac{1}{2}}}{{\left( {\frac{{{n_c}}}{n}} \right)}^2}}}{{{{({hd} )}^2}\left( {1 + {{\left( {\frac{{{n_c}}}{n}} \right)}^4}} \right) - {V^2}{{\left( {\frac{{{n_c}}}{n}} \right)}^4}}}\\ V = \left( {\frac{{2\pi }}{\lambda }} \right)d{({n_c^2 - {n^2}} )^{\frac{1}{2}}} \end{array}$$

where V is the V-number of the waveguide and d is the transverse width of the waveguide. To construct the matrix M in Eqs. (4)–(5) we evaluate Eq. (8) at each spatial grid position for each mode found by solving Eq. (9). For a 50-µm guide, vacuum wavelength of 1.537 µm, a core index of 3.48 (e.g., Si at this wavelength), and a cladding index of 1.44, there are N=206 transverse modes. We simulate this pattern at 1000 equally spaced points across the waveguide giving an M matrix with dimensions of 1000×206 and a W matrix with dimensions of 1000×1000. Each mode has a different propagation coefficient β_j and these range from 14.2 µm^-1 for the first mode down to 5.9 µm^-1 for the last. Since the guide length L is taken to be 20 cm, the phase β_jL modulo 2π is effectively random for each mode. Because it would be unrealistic to split out 1000 output waveguides 50 nm wide from the 50-µm multimode waveguide, we simulate the readout using detector pixels that bin the output from 5 spatial points.

As a test of the ring waveguide resonator for reservoir computing we used solutions to the Mackey-Glass (MG) delay differential equation given by

(10)$$\frac{{\textrm{d}u}}{{\textrm{d}t}} = \frac{{u({t - \tau } )}}{{5[{u{{({t - \tau } )}^{10}} + 1} ]}} - \frac{{u(t )}}{{10}}$$

where u(t) is the amplitude of the MG waveform and $\tau $ is the delay parameter [32]. Typical reservoir computing calculations train on the MG waveform for multiple delays and then perform operations such as replication, prediction or labelling on test samples of the MG waveform. Here we illustrate the performance of the ring resonator RC for prediction as we have found prediction to be the hardest of these three tasks. Prediction is performed by using the trained RC to perform free-running predictions of the subsequent values to the waveform and recursively feeding those results back into the RC as the input data. Calculations of the predictive performance for a conventional RC as described by reviews such as Lukoševičius et. al. [27] with about 200 neurons and 2500 training points show that it is possible to predict about 500 time steps for this problem as shown in Fig. 2.

Fig. 2. Prediction of a Mackey-Glass waveform using simulations of a conventional RC with pseudo-random matrices W and W_in. a) The true (blue trace) and predicted (red trace) MG waveforms for a delay of τ = 17. b) Prediction error (difference of red and blue curves in (a)) for the MG waveform with delay parameter of 17.

Download Full Size | PDF

To compare performance of a conventional RC to the ring resonator RC, we model the 50-µm waveguide described above. The waveform u(t) is injected across the waveguide with a random pattern whose spatial scale is 0.25 µm as shown in Fig. 3(a), and no additional bias is used (we note that the numerical solutions to Eq. (1)0 are naturally biased). Injecting this random pattern means that nearly all 206 modes of the 50-µm guide are excited although not with equal amplitude (Fig. 3(b)). The exact spatial profile of the injected signal does not seem to matter much as long as a large fraction of the modes are excited. The nonlinear activation and non-idealities in an experimental system such as bending and strain will excite the small minority of modes that the injection pattern may not excite. For example, injecting the waveform with a single 5-µm spot excites a large fraction of the modes (Fig. 3(c)-(d)) and gives similar performance to the random pattern although the number of time steps for good prediction is less (about 300 steps) and the hyperparameters are different.

Fig. 3. Spatial and modal distributions for signal injection patterns. a) Example spatial pattern for a random pattern injection with each random value occupying a 0.25-µm section of the waveguide input. b) Mode amplitudes for the random injection condition. c) Spatial pattern for injecting the input signal uniformly across a 5-µm section of the waveguide. d) Mode amplitudes for injection into the 5-µm section.

Download Full Size | PDF

Figure 4 shows the predicted waveform, the true waveform and the error for MG delays of 15, 17, 19, and 21 using hyperparameters shown in Table 1 and 2500 training points. The predictive performance for an MG delay of 17 shown in Fig. 4(c)-(d) for the ring resonator RC is very similar to the performance shown in Fig. 2 for the conventional RC with a number of neurons on the order of the number of waveguide modes. The quality of the prediction degrades as the MG delay increases as expected as the MG waveform becomes more chaotic at larger delays. We note that the optimal hyperparameters in Table 1 are obtained by finding the best prediction of the test data, and that, as is typical for machine learning tasks, the optimal hyperparameters can vary substantially for different tasks. Since systems that are extremely sensitive to physical hyperparameters such as outcoupling and saturation intensity are not very useful, we evaluated the predictive performance over the first 100 steps as a function of outcoupling, saturation intensity, and the ridge regression regularization parameter as shown in Fig. 5. Note the broad operating ranges in all three parameters.

Fig. 4. Predictions of Mackey-Glass waveforms using the simulated planar waveguide ring resonator RC. a) The true (blue trace) and predicted (red trace) MG waveforms for a delay parameter of 15. b) Prediction error for the MG waveform with delay parameter of 15. The remaining rows of panels show the same for MG delay parameters of 17 (c-d), 19 (e-f), and 21 (g-h).

Download Full Size | PDF

Fig. 5. Root mean square error for the first 100 points to the MG prediction with a delay parameter of 17 as a function of the a) outcoupling fraction, b) saturation intensity, and c) regularization hyperparameters.

Download Full Size | PDF

Table 1. Hyperparameters used for the Mackey-Glass prediction in the planar waveguide RC.

View Table | View all tables in this article

4. Multimode fiber reservoir computer simulation

A key advantage of the ring resonator RC is that its power and processing time requirements scale much more slowly with number of neurons than a conventional RC. The number of effective neurons can be increased by simply using a waveguide with more modes and an optical sensor with more pixels. The power and time performance of the camera digitization and prediction evaluation will scale as O(N) at worst, while the time complexity of repeated multiplications by a dense N×N matrix in a typical RC scales as O(N²) [21]. Thus, a key application of the photonic ring resonator RC is to a problem requiring tens of thousands of neurons or more with dense connectivity. To demonstrate the ability of a large optical ring resonator RC to perform more complex predictions with large multivariate input signals u(t), we simulate an RC based on speckle in multimode optical fiber and use it to perform free-running predictions of the KS equation, a benchmark RC test problem [33,34].

Like the planar waveguide, the field profiles and propagation constants are calculated for the LP_lm modes of a cylindrical waveguide in the weakly guiding case, which can be cast in the same general form as Eq. (3) [35]. Briefly, the valid modes and their propagation constants are calculated from the characteristic equation:

(11)$$\begin{array}{l} u\left[ {\frac{{{J_{l - 1}}(u)}}{{{J_l}(u)}}} \right] ={-} w\left[ {\frac{{{K_{l - 1}}(w)}}{{{K_l}(w)}}} \right]\\ {u_{lm}} = a{({{k^2}n_c^2 - \beta_{lm}^2} )^{{\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}} \right.}\!\lower0.7ex\hbox{$2$}}}}\\ {w_{lm}} = a{({\beta_{lm}^2 - {k^2}n_{}^2} )^{{\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}} \right.}\!\lower0.7ex\hbox{$2$}}}} \end{array}$$

where a is the fiber core radius, n_c is the core refractive index, n is the cladding refractive index assumed to extend to infinity, β_lm is the propagation constant, and k is the laser wavenumber in vacuum. The characteristic equation has multiple solutions for each value of l indexed by the m subscript. J_ν(x) and K_ν(x) are the Bessel function and modified Bessel function of the second kind, respectively. Each (l,m) combination corresponds to a total of 4 modes, 2 orientations for each of two polarization directions. For this model, we assume the light is polarized and just model the two orientations for a single polarization. The transverse electric fields are given by:

(12)$${E_{lm}} = {E_0}\left\{ {\begin{array}{c} {\frac{{{J_l}({{{{u_{lm}}r} / a}} )}}{{{J_l}({u_{lm}})}}\cos (l\phi ),\quad r < a}\\ {\frac{{{K_l}({{{{w_{lm}}r} / a}} )}}{{{K_l}({w_{lm}})}}\cos (l\phi ),\quad r > a} \end{array}} \right.$$

and

(13)$${E_{lm}} = {E_0}\left\{ {\begin{array}{c} {\frac{{{J_l}({{{{u_{lm}}r} / a}} )}}{{{J_l}({u_{lm}})}}\sin (l\phi ),\quad r < a}\\ {\frac{{{K_l}({{{{w_{lm}}r} / a}} )}}{{{K_l}({w_{lm}})}}\sin (l\phi ),\quad r > a} \end{array}} \right..$$

The remainder of the simulation is identical to the linear ring resonator, with the electric fields of each mode and the propagation constants used to construct a complex W matrix that propagates the electric field in the fiber one lap around the ring resonator. To efficiently calculate the spatial profiles of each mode, we use a polar grid of points instead of the typical cartesian grid. This choice leverages the polar symmetry of the modes to substantially reduce the computation time. The polar grid contains 75 equally spaced radial points extending out to 1.05 times the fiber core radius and 236 equally spaced angular points for a total of P=17,700 spatial points. As with the planar waveguide in Section 3, we calculate the M matrix by evaluating the field equations at each position for each solution to the boundary conditions, yielding a matrix with dimensions 17700×N, where N is the number of modes supported by the fiber (see Table 2). The resulting W matrix has dimensions of 17,700×17,700. Unlike the planar waveguide from Section 3, we can magnify the image of the fiber output to an arbitrary size, so we do not bin the spatial points into larger pixels for this simulation.

Table 2. List of parameters for the simulated multimode optical fibers

View Table | View all tables in this article

The test problem for this system is the Kuramoto-Sivashinsky (KS) equation, which defines a scalar field y(x,t) as

(14)$$\frac{{\partial y}}{{\partial t}} ={-} y\frac{{\partial y}}{{\partial x}} - \frac{{{\partial ^2}y}}{{\partial {x^2}}} - \frac{{{\partial ^4}y}}{{\partial {x^4}}}.$$

The input data, u_n consists of the solution to the KS equation at time point n evaluated at 64 spatial points. Each of the 64 spatial input channels plus the bias channel for a total of 65 inputs, are assigned to illuminate one square of a 9 × 9 grid fully enclosed in the fiber core as shown in Fig. 6(a). The input signal is calculated by evenly dividing a total intensity proportional to the input value among all the spatial points lying inside the corresponding square for each channel.

Fig. 6. a) Illustration of the distribution of square pixels used as input channels in the circular fiber core. b) Excerpt from the K-S equation solution used as the dataset for testing the fiber ring resonator RC. The time axis is normalized to the Lyapunov time.

Download Full Size | PDF

A solution to the KS equation with periodic boundary conditions is generated using the method in [36], the system parameters used in previous reservoir computing work [33,34], and a random initial condition with the last 2 points on each side of the spatial domain set to enforce the periodic boundary conditions. The solution forms braided wave structures containing irregular intensity peaks and valleys that split and merge chaotically. As with the Mackey-Glass prediction problem in the previous section, the goal is to train the network on a set of data to predict the state of the waveform at the subsequent time point then make free-running predictions by recursively feeding the predicted state back into the network.

As shown in Fig. 6(b), the waveform dynamics vary significantly over time, with long periods of irregular but relatively stable oscillations interleaved with periods of more chaotic behavior. To enable comparisons with other RC predictions of the KS equation in the literature [21,33] time axes in Fig. 6(b) and the figures below are normalized to the Lyapunov time, a measure of the time it takes for two similar initial conditions to decorrelate, following the convention in those references. The Lyapunov time of this system is 23.2, and we use a time step of 0.25 in these simulations, so each time step is about 0.011 times the Lyapunov time. To get a reliable understanding of the predictive power of the ring resonator RC, we use the following scheme to perform 100 different test runs from a single training run and analyze the error statistics across the runs. A long data set of 250,000 time points is calculated by solving the KS equation starting from a single initial condition, and a constant is added to the resulting data so that all the numbers are positive as required by the optical RC.

To train the network, we initialize the reservoir with all neuron states set to 0 and run the simulated RC feeding the first 100,000 points of the dataset as u_in. The first 500 points are discarded, and we train W_out using ridge regression with the remaining reservoir states for a total of 99,500 training points. The state of the network (x_n) at the end of the training stage is saved for reuse. The remaining 150,000 time points in the generated KS solution are divided into 100 test sections of 1500 points each. For each test, we reinitialize the reservoir state using the saved state from the end of the training phase creating a scenario for each test analogous to the typical procedure for a single test phase, where the RC transitions to testing immediately after the training phase without restarting. We found that the RC performs better with this initialization than it does when all neurons are reinitialized at 0. The network is then run by feeding in truth data for the first 500 time points in the test section. This discards any “history” data in the network from the training phase and replaces it with the correct state for that test. Finally, the free-running prediction test is performed by running the network for 1000 time steps while feeding the prediction back into the input. The prediction results are then compared to the remaining 1000 points in each test section.

Although the number of neurons in the W matrix is nominally equal to the number of spatial pixels of the measured speckle pattern, the predictive power of the photonic RC depends strongly on the number of modes in the waveguide. To demonstrate this, we run the simulation for 2 fibers with parameters shown in Table 2 using the same spatial discretization with significantly more pixels than modes to ensure that all of the information in the speckle patterns is captured. Table 3 shows the RC hyperparameters used in these calculations. As with the planar waveguide, the prediction quality is robust to small changes in these hyperparameters. The spatially averaged root-mean-square (RMS) error is calculated for each time point in each test, and the tests are sorted by the first time point when the RMS error exceeds 1.5 to systematically aggregate the results. The network shows highly variable performance across the test sets, generally correlated with the complexity of the waveform. Figure 7 shows some typical test results for the 4238-mode fiber that illustrate the range of possibilities. Figure 7(a)-(c) show the 10^th test to cross the 1.5 RMS error threshold, effectively the 10^th worst performing test case out of 100. Figure 7(d)-(f) show the 50^th case to cross the error threshold, and Fig. 7(g)-(i) show one of the tests that had not crossed the threshold by the end of the 1000 point window. The plots show how the trajectories tend to have a discrete moment where the waveform pattern clearly diverges from the truth data, and setting an RMS error threshold of 1.5 for sorting the tests is intended to approximately track that divergence. The performance of the median case (Fig. 7(d)-(f)) is similar to other RC predictions of the KS equations shown in the literature with a numbers of neurons on the order of the number of modes in the fiber [21,33], demonstrating the effectiveness of the simulated ring resonator photonic RC.

Fig. 7. Example predictions of the solution to the KS equation using the simulated fiber ring resonator with 4238 modes. The panels in the top row a,d,g) show the truth data, the middle panels b,e,h) show the simulated fiber RC predictions, and the bottom panels c,f,i) show the prediction error. Panels a-c) show the 10^th test sample to cross the RMS error threshold of 1.5, panels d-f) show the 50^th sample to cross the error threshold, and panels g-i) show one of the samples to reach the end of the 1000 time step test window without crossing the threshold.

Download Full Size | PDF

Table 3. Hyperparameters used for both optical fiber RC simulations

View Table | View all tables in this article

The plots in Fig. 8 show the aggregated statistics from all 100 test cases for both fibers. Figure 8(a)-(b) show several percentiles of the RMS error distribution at each time point for the larger and smaller fiber respectively, and Fig. 8(c) shows the cumulative distribution of the first time point where the test cases exceed the 1.5 RMS error threshold. The fiber with 4238 modes shows excellent performance on this problem, with the median test case accurately tracking the waveform for around 550 time points, and the 10 best test cases tracking all the way to the end of 1000 time point window. The 246-mode fiber showed somewhat worse performance, with the median test case crossing the RMS error threshold after 350 time points, and only 3 cases reaching the end of the 1000-point window. This data shows that the number of modes in the waveguide has a substantial influence on the RC performance and affects the number of “effective neurons” in the reservoir.

Fig. 8. Aggregate statistics for all 100 test runs. a) Percentile values of RMS error for the 4238-mode fiber across the 100 test runs. b) Percentile values for the 246 mode fiber. c) Cumulative number of test runs that have crossed the RMS error threshold of 1.5 at each time step for both fiber sizes.

Download Full Size | PDF

5. Physical realizations of a multimode ring resonator reservoir computer

The multimode ring resonator RC is a hybrid of photonic and electronic parts, where the matrix multiplication, feedback, and nonlinear activation are performed entirely in the optical domain while the output speckle patterns are measured using an array of photodetectors or a camera with the appropriate analog-digital converters (ADCs). As is the case in most photonic RC architectures, the output layer evaluation and training are done digitally using a dedicated FPGA or digital computer. Most of the passive photonic components required to realize the proposed architecture can be designed and fabricated on an SOI platform for a chip-scale device. However, due to the indirect bandgap of silicon, light generation, amplification/saturation, and detection will require the heterogenous integration of III-V materials or other active media to realize a full implementation of the proposed architecture [37].

The input data, u(t), is injected into the reservoir using an array of electro-optical modulators (such as a ring modulator in SOI) that transfer the electrical signal of interest into the optical domain through amplitude or phase modulation. Unlike existing speckle-based RCs [6–8,10,21] that require the number of modulators equal the sum of the number of feedback neurons plus the dimensions of the input signal, this proposed architecture requires that the number of modulators equal the dimensions of the input, which has to be much less than the number of neurons for effective RC operation. The input information gets mapped into the larger reservoir as the injected light propagates around the MM waveguide. Our earlier studies have shown that the MM waveguide only needs to be a few cm long to function as a pseudo-random matrix [10,38]. Several technologies for implementing the nonlinear activation have been proposed using semiconductor optical amplifier (SOA) configurations [15,39–42], cavity-enhanced 2-photon absorption in GaAs [29], micro-ring resonators [43], or nanophotonic structures [30]. Nonlinear elements requiring a reflective geometry can be used with a fiber optic circulator to maintain the direction of propagation around the ring resonator [2]. For readout, the speckle pattern at the output of the multimode waveguide is mapped to an array of single mode waveguides [38] where they can then be split using an array of Y-splitters with one arm routed to on-chip germanium detectors for reading out neurons and the other fed back into the multimode waveguide to create a recurrent reservoir computer. The light intensity incident on the detectors can be optimized as necessary by adjusting the input pump laser power, or with an additional attenuation or amplification stage in the output waveguides. Alternatively, an integrated beamsplitter [44,45] in the MM waveguide would eliminate the need for the MM to single mode transition. As the majority of the existing literature on PIC-related applications are focused around single-mode waveguides, further work will be needed to demonstrate the large-scale multimode beamsplitter, the evanescently coupled detector array, and the various proposed nonlinear activation systems.

Most of the elements described above are well-developed in optical fiber systems and can be used to assemble a device as depicted in Fig. 1. Fiber beamsplitters are readily available to handle the splitting of optical power out of the ring resonator for detection and coupling the modulated input signals into the ring resonator, although custom solutions may be required to accommodate the large MM fiber. For the output, light in the ring resonator can simply be split out into a short length of MM fiber, and the light exiting the end of the fiber can be imaged onto a camera or photodiode array in free space. As with the PIC geometry, the output light can be attenuated with a neutral density filter or amplified with a fiber amplifier as necessary to match the camera’s dynamic range. For the input, the modulated inputs can be imaged onto the end of the MM fiber via a lens system as used in our previous work [10]. A nonlinear activation element similar to those described above for the chip-scale device would have to be designed and included in the system.

The final design consideration is digitization of the optical signals, and how that relates to the length of the ring resonator. Our calculations show that a minimal length of MM waveguide can mimic the random matrix multiplication. Rather, the length of the waveguide must be set to approximately match the digitization rate of the ADCs, and likewise the ADCs must be chosen to match a reasonable resonator length. For a PIC device, reasonable arrays of ADCs could operate at tens of MSa/s, so a resonator of 1-5 m would be required. For the optical fiber system where tens of thousands of pixels are desired, commercially available high-speed CMOS cameras can operate at speeds of 50-100 kHz when windowed to that number of pixels. This would require a 2-4 km fiber spool. Optical losses are low enough to support these lengths, but it will be important to properly pot the fiber and mount it in a thermally controlled housing to ensure stable operation of the RC. Challenges of fabricating a PIC waveguide or stabilizing an optical fiber of the required lengths could be addressed by alternate reformulations that allow the signal to make multiple trips around the ring resonator for each input data point.

6. Conclusions

In this work, we demonstrated that a multimode waveguide ring resonator is mathematically equivalent to a reservoir computer by showing how the equations for light propagation through the waveguide combined with a nonlinear activation process are mathematically equivalent to those of a standard reservoir computer. We then simulated reservoir computers derived from both planar and cylindrical waveguides to demonstrate that this system can accomplish tasks with similar performance to a conventional RC. We used the simulated planar waveguide RC to produce robust predictions of the Mackey-Glass equation and the cylindrical fiber waveguide to do the same for the 2-dimensional Kuramoto-Sivashinsky equation. We conclude by describing how these results inform the design of physical devices. Lightweight computational tasks requiring hundreds of neurons could be accomplished with a chip-scale PIC device, while larger tasks requiring tens of thousands of neurons could be performed using a relatively low-power fiber device instead of high-power computational hardware.

Acknowledgements

The authors would like to thank Mushegh Rafayelyan for providing MATLAB code for generating KS equation solutions and T. Justin Shaw for helpful discussions.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. G. Tanaka, T. Yamane, J. B. Héroux, R. Nakane, N. Kanazawa, S. Takeda, H. Numata, D. Nakano, and A. Hirose, “Recent advances in physical reservoir computing: A review,” Neural Networks 115, 100–123 (2019). [CrossRef]

2. G. Van der Sande, D. Brunner, and M. C. Soriano, “Advances in photonic reservoir computing,” Nanophotonics 6(3), 561–576 (2017). [CrossRef]

3. X. Lin, Y. Rivenson, N. T. Yardimci, M. Veli, Y. Luo, M. Jarrahi, and A. Ozcan, “All-optical machine learning using diffractive deep neural networks,” Science 361(6406), 1004–1008 (2018). [CrossRef]

4. Y. Shen, N. C. Harris, S. Skirlo, M. Prabhu, T. Baehr-Jones, M. Hochberg, X. Sun, S. Zhao, H. Larochelle, D. Englund, and M. Soljačić, “Deep learning with coherent nanophotonic circuits,” Nat. Photonics 11(7), 441–446 (2017). [CrossRef]

5. J. Bueno, S. Maktoobi, L. Froehly, I. Fischer, M. Jacquot, L. Larger, and D. Brunner, “Reinforcement learning in a large-scale photonic recurrent neural network,” Optica 5(6), 756–760 (2018). [CrossRef]

6. J. Dong, S. Gigan, F. Krzakala, and G. Wainrib, “Scaling Up Echo-State Networks With Multiple Light Scattering,” in 2018 IEEE Statistical Signal Processing Workshop (SSP) (2018), pp. 448–452.

7. J. Dong, M. Rafayelyan, F. Krzakala, and S. Gigan, “Optical Reservoir Computing Using Multiple Light Scattering for Chaotic Systems Prediction,” IEEE J. Sel. Top. Quantum Electron. 26(1), 1–12 (2020). [CrossRef]

8. A. C. Scofield, G. A. Sefler, T. J. Shaw, and G. C. Valley, “Recent results using laser speckle in multimode waveguides for random projections,” in Optical Data Science II (International Society for Optics and Photonics, 2019), 10937, p. 109370B.

9. K. Vandoorne, P. Mechet, T. Van Vaerenbergh, M. Fiers, G. Morthier, D. Verstraeten, B. Schrauwen, J. Dambre, and P. Bienstman, “Experimental demonstration of reservoir computing on a silicon photonics chip,” Nat. Commun. 5(1), 3541 (2014). [CrossRef]

10. U. Paudel, M. Luengo-Kovac, J. Pilawa, T. J. Shaw, and G. C. Valley, “Classification of time-domain waveforms using a speckle-based optical reservoir computer,” Opt. Express 28(2), 1225–1237 (2020). [CrossRef]

11. M. Miscuglio and V. J. Sorger, “Photonic tensor cores for machine learning,” Appl. Phys. Rev. 7(3), 031404 (2020). [CrossRef]

12. P. Prucnal and B. Shastri, Neuromorphic Photonics (CRC, 2017).

13. F. Duport, A. Smerieri, A. Akrout, M. Haelterman, and S. Massar, “Fully analogue photonic reservoir computer,” Sci. Rep. 6(1), 22381 (2016). [CrossRef]

14. A. Argyris, J. Bueno, and I. Fischer, “Photonic machine learning implementation for signal recovery in optical communications,” Sci. Rep. 8(1), 8487 (2018). [CrossRef]

15. F. Duport, B. Schneider, A. Smerieri, M. Haelterman, and S. Massar, “All-optical reservoir computing,” Opt. Express 20(20), 22783–22795 (2012). [CrossRef]

16. L. Larger, M. C. Soriano, D. Brunner, L. Appeltant, J. M. Gutierrez, L. Pesquera, C. R. Mirasso, and I. Fischer, “Photonic information processing beyond Turing: an optoelectronic implementation of reservoir computing,” Opt. Express 20(3), 3241–3249 (2012). [CrossRef]

17. L. Larger, A. Baylón-Fuentes, R. Martinenghi, V. S. Udaltsov, Y. K. Chembo, and M. Jacquot, “High-Speed Photonic Reservoir Computing Using a Time-Delay-Based Architecture: Million Words per Second Classification,” Phys. Rev. X 7(1), 011015 (2017). [CrossRef]

18. D. Brunner, M. C. Soriano, C. R. Mirasso, and I. Fischer, “Parallel photonic information processing at gigabyte per second data rates using transient states,” Nat. Commun. 4(1), 1364 (2013). [CrossRef]

19. S. Ortín, M. C. Soriano, L. Pesquera, D. Brunner, D. San-Martín, I. Fischer, C. R. Mirasso, and J. M. Gutiérrez, “A Unified Framework for Reservoir Computing and Extreme Learning Machines based on a Single Time-delayed Neuron,” Sci. Rep. 5(1), 14945 (2015). [CrossRef]

20. C. Mesaritakis and D. Syvridis, “Reservoir computing based on transverse modes in a single optical waveguide,” Opt. Lett. 44(5), 1218–1221 (2019). [CrossRef]

21. M. Rafayelyan, J. Dong, Y. Tan, F. Krzakala, and S. Gigan, “Large-Scale Optical Reservoir Computing for Spatiotemporal Chaotic Systems Prediction,” Phys. Rev. X 10(4), 041037 (2020). [CrossRef]

22. P. Antonik, N. Marsal, and D. Rontani, “Large-Scale Spatiotemporal Photonic Reservoir Computer for Image Classification,” IEEE J. Sel. Top. Quantum Electron. 26(1), 1–12 (2020). [CrossRef]

23. F. Laporte, A. Katumba, J. Dambre, and P. Bienstman, “Numerical demonstration of neuromorphic computing with photonic crystal cavities,” Opt. Express 26(7), 7955–7964 (2018). [CrossRef]

24. M. Pisati, F. D. Bernardinis, P. Pascale, C. Nani, M. Sosio, E. Pozzati, N. Ghittori, F. Magni, M. Garampazzi, G. Bollati, A. Milani, A. Minuti, F. Giunco, P. Uggetti, I. Fabiano, N. Codega, A. Bosi, N. Carta, D. Pellicone, G. Spelgatti, M. Cutrupi, A. Rossini, R. Massolini, G. Cesura, and I. Bietti, “6.3 A Sub-250mW 1-to-56Gb/s Continuous-Range PAM-4 42.5dB IL ADC/DAC-Based Transceiver in 7nm FinFET,” in 2019 IEEE International Solid- State Circuits Conference - (ISSCC) (2019), pp. 116–118.

25. L. Kull, D. Luu, C. Menolfi, M. Braendli, P. A. Francese, T. Morf, M. Kossel, A. Cevrero, I. Ozkaya, and T. Toifl, “A 24-to-72GS/s 8b time-interleaved SAR ADC with 2.0-to-3.3pJ/conversion and >30dB SNDR at nyquist in 14nm CMOS FinFET,” in 2018 IEEE International Solid - State Circuits Conference - (ISSCC) (2018), pp. 358–360.

26. Boris Murmann, “ADC Performance Survey 1997-2020,” http://web.stanford.edu/∼murmann/adcsurvey.html.

27. M. Lukoševičius and H. Jaeger, “Reservoir computing approaches to recurrent neural network training,” Comput. Sci. Rev. 3(3), 127–149 (2009). [CrossRef]

28. H. Jaeger, “The” echo state” approach to analysing and training recurrent neural networks-with an erratum note’,” Bonn, Germany: German National Research Center for Information Technology GMD Technical Report 148, (2001).

29. J. H. Vella, J. H. Goldsmith, A. T. Browning, N. I. Limberopoulos, I. Vitebskiy, E. Makri, and T. Kottos, “Experimental Realization of a Reflective Optical Limiter,” Phys. Rev. Appl. 5(6), 064010 (2016). [CrossRef]

30. M. Miscuglio, A. Mehrabian, Z. Hu, S. I. Azzam, J. George, A. V. Kildishev, M. Pelton, and V. J. Sorger, “All-optical nonlinear activation function for photonic neural networks [Invited],” Opt. Mater. Express 8(12), 3851–3863 (2018). [CrossRef]

31. A. Yariv, “Coupled-mode theory for guided-wave optics,” IEEE J. Quantum Electron. 9(9), 919–933 (1973). [CrossRef]

32. M. C. Mackey and L. Glass, “Oscillation and chaos in physiological control systems,” Science 197(4300), 287–289 (1977). [CrossRef]

33. J. Pathak, B. Hunt, M. Girvan, Z. Lu, and E. Ott, “Model-Free Prediction of Large Spatiotemporally Chaotic Systems from Data: A Reservoir Computing Approach,” Phys. Rev. Lett. 120(2), 024102 (2018). [CrossRef]

34. M. Rafayelyan, J. Dong, Y. Tan, F. Krzakala, and S. Gigan, “Optical reservoir computing for high-dimensional spatio-temporal chaotic systems prediction (Conference Presentation),” in AI and Optical Data Sciences (International Society for Optics and Photonics, 2020), 11299, p. 112990B.

35. D. Gloge, “Weakly Guiding Fibers,” Appl. Opt. 10(10), 2252–2258 (1971). [CrossRef]

36. A.-K. Kassam and L. N. Trefethen, “Fourth-Order Time-Stepping for Stiff PDEs,” SIAM J. Sci. Comput. 26(4), 1214–1233 (2005). [CrossRef]

37. J. E. Bowers, T. Komljenovic, M. Davenport, J. Hulme, A. Y. Liu, C. T. Santis, A. Spott, S. Srinivasan, E. J. Stanton, and C. Zhang, “Recent advances in silicon photonic integrated circuits,” in Next-Generation Optical Communication: Components, Sub-Systems, and Systems V (International Society for Optics and Photonics, 2016), 9774, p. 977402.

38. D. B. Borlaug, S. Estrella, C. Boone, G. A. Sefler, T. J. Shaw, A. Roy, L. Johansson, and G. C. Valley, “A Photonic Integrated Circuit Based Compressive Sensing Radio Frequency Receiver Using Waveguide Speckle,” (2020).

39. S. Keyvaninia, G. Roelkens, D. V. Thourhout, G. H. Duan, M. Lamponi, F. Lelarge, J.- Fedeli, S. Messaoudene, E. J. Geluk, and B. Smalbrugge, “A highly efficient electrically pumped optical amplifier integrated on a SOI waveguide circuit,” in The 9th International Conference on Group IV Photonics (GFP) (2012), pp. 222–224.

40. J. Rönn, W. Zhang, A. Autere, X. Leroux, L. Pakarinen, C. Alonso-Ramos, A. Säynätjoki, H. Lipsanen, L. Vivien, E. Cassan, and Z. Sun, “Ultra-high on-chip optical gain in erbium-based hybrid slot waveguides,” Nat. Commun. 10(1), 432 (2019). [CrossRef]

41. F. E. Doany, R. Budd, L. Schares, T. Huynh, M. Wood, D. Kuchta, N. Dupuis, C. Schow, B. Lee, M. Moehrle, A. Sigmund, W. Rehbein, T. Y. Liow, L. W. Luo, and G. Q. Lo, “A Four-Channel Silicon Photonic Carrier with Flip-Chip Integrated Semiconductor Optical Amplifier (SOA) Array Providing >10-dB Gain,” in 2016 IEEE 66th Electronic Components and Technology Conference (ECTC) (2016), pp. 1061–1068.

42. G. Mourgias-Alexandris, A. Tsakyridis, N. Passalis, A. Tefas, K. Vyrsokinos, and N. Pleros, “An all-optical neuron with sigmoid activation function,” Opt. Express 27(7), 9620–9630 (2019). [CrossRef]

43. S. Yan, J. Dong, A. Zheng, and X. Zhang, “Chip-integrated optical power limiter based on an all-passive micro-ring resonator,” Sci. Rep. 4(1), 6676 (2015). [CrossRef]

44. H. Xu, D. Dai, and Y. Shi, “Ultra-broadband on-chip multimode power splitter with an arbitrary splitting ratio,” OSA Continuum 3(5), 1212–1221 (2020). [CrossRef]

45. D. Ahn, L. C. Kimerling, and J. Michel, “Evanescent Coupling Device Design for Waveguide-Integrated Group IV Photodetectors,” J. Lightwave Technol. 28(23), 3387–3394 (2010).

Hyperparameter	Delay = 15	17	19	21
Outcoupling Fraction $(η)$	0.05	0.03	0.01	0.01
Saturation Intensity $(γ)$	0.81	0.83	0.83	0.8
Bias Amplitude	0	0	0	0
Regularization Parameter	10⁻⁷	10⁻⁷	5 × 10⁻⁶	5 × 10⁻⁶

Parameter	Fiber 1	Fiber 2
Core Refractive Index	1.45	1.45
Numerical Aperture	0.22	0.10
Core Diameter	200 μm	105 μm
Laser Wavelength	1064 nm	1064 nm
Number of Modes	4238	246
Fiber Length	1 m	1 m

Hyperparameter	Value
Outcoupling Fraction $(η)$	0.69
Saturation Intensity $(γ)$	0.28
Bias Amplitude	0.72
Regularization Parameter	10^−8.5

Hyperparameter	Delay = 15	17	19	21
Outcoupling Fraction $(η)$	0.05	0.03	0.01	0.01
Saturation Intensity $(γ)$	0.81	0.83	0.83	0.8
Bias Amplitude	0	0	0	0
Regularization Parameter	10⁻⁷	10⁻⁷	5 × 10⁻⁶	5 × 10⁻⁶

Parameter	Fiber 1	Fiber 2
Core Refractive Index	1.45	1.45
Numerical Aperture	0.22	0.10
Core Diameter	200 μm	105 μm
Laser Wavelength	1064 nm	1064 nm
Number of Modes	4238	246
Fiber Length	1 m	1 m

Photonic reservoir computer using speckle in multimode waveguide ring resonators

Abstract

1. Introduction

2. Ring resonator reservoir computer simulation

3. Planar waveguide reservoir computer simulations

4. Multimode fiber reservoir computer simulation

5. Physical realizations of a multimode ring resonator reservoir computer

6. Conclusions

Acknowledgements

Disclosures

Data availability

References

Data availability

Cited By

Figures (8)

Tables (3)

Equations (14)

Optics Express