Using multidimensional speckle dynamics for high-speed, large-scale, parallel photonic computing

The recent rapid increase in demand for data processing has resulted in the need for novel machine learning concepts and hardware. Physical reservoir computing and an extreme learning machine are novel computing paradigms based on physical systems themselves, where the high dimensionality and nonlinearity play a crucial role in the information processing. Herein, we propose the use of multidimensional speckle dynamics in multimode fibers for information processing, where input information is mapped into the space, frequency, and time domains by an optical phase modulation technique. The speckle-based mapping of the input information is high-dimensional and nonlinear and can be realized at the speed of light; thus, nonlinear time-dependent information processing can successfully be achieved at fast rates when applying a reservoir-computing-like-approach. As a proof-of-concept, we experimentally demonstrate chaotic time-series prediction at input rates of 12.5 Gigasamples per second. Moreover, we show that owing to the passivity of multimode fibers, multiple tasks can be simultaneously processed within a single system, i.e., multitasking. These results offer a novel approach toward realizing parallel, high-speed, and large-scale photonic computing.


Introduction
Reservoir computing (RC) [1,2,3] and an extreme learning machine (ELM) [4] are novel neuro-inspired computing methods that use high-dimensional and nonlinear systems, e.g., neural networks, to process input information. Unlike a conventional neural network, internal networks are fixed, and only the readout weights are trained, which greatly simplifies the training method and enables easy hardware implementation. To date, a variety of physical implementations of RC and ELM, including optoelectronic and photonic systems [5], memristors [6], spin waves [7], and soft-robots [8], have been proposed, and have demonstrated a performance comparable with that of digital computing based on other algorithms in a series of benchmark tasks, including time-series prediction, wireless channel equalization, and image and phoneme recognition [9].
Photonics-based RC is of particular interest. The use of photonic systems for RC implementation is expected to accelerate recurrent neural network processing with low energy costs [10,11]. Most photonic RC systems are based on delayed feedback structures, which can be simply realized using a nonlinear device with a time-delayed feedback loop [12,13,14,15,16,17,18]. Based on a time-multiplexing technique [19,20], the sampled signals of the nonlinear device are used as virtual neurons (network nodes) for information processing. With this approach, most complex tasks can be solved with few errors; however, these systems suffer from an inherent trade-off between the number of neurons and the required processing time.
Meanwhile, significant progress has been made with various multiplexing approaches that utilize other degrees of freedom of light, i.e., space and wavelength, to overcome the above trade-off [21,22,23,24,25,26,27]. Most multiplexing techniques are based on optical data transportation techniques used for optical fiber communication and have the potential to realize high-speed parallel operation. Spatial multiplexing schemes, based on optical node arrays, have been used for RC in silicon photonic chips and have demonstrated high-speed processing of the bit header recognition at up to 12.5 Gbit/s [21]. To generate a larger number of spatial nodes (i.e., a larger network scale), the use of speckles generated from light scattering [24] and multimode waveguides [25], and diffractively coupled systems [26], have been proposed, demonstrating good performance even for complex tasks, such as chaotic timeseries prediction. As other approaches, frequency multiplexing schemes have also been proposed, where the amplitude and phase of the frequency sidebands of a laser in a fiber loop are used as reservoir nodes, demonstrating nonlinear channel equalization and speech recognition [27]. These approaches are promising; however, there are certain limitations. For example, the constructions of the abovementioned spatial multiplexing networks essentially involve a direct trade-off between parallelism and footprint. In frequency multiplexing, the number of nodes may be practically limited by the amplitudes or bandwidth of the electronic modulators.
Herein, we propose a new approach to photonic parallel information processing based on highdimensional speckle dynamics generated from multimode fibers. This approach easily enables the combined use of space, wavelength, and time multiplexing to achieve high-speed, large-scale, and parallel processing. In the proposed system, high nonlinearity required to process nonlinear tasks can be introduced using fast optical phase modulation over a few Gigasamples/s (GS/s). The origin of the nonlinearity is the nonlinear optical mapping from phase-encoded input information into speckle patterns. This differs from other types of photonic RC using multimode fiber speckles [25], where nonlinearity is electronically introduced. As a proof-of-concept demonstration for high-speed parallel processing, we experimentally show that a chaotic time-series prediction can be performed at 12.5 GS/s. Moreover, the proposed approach enables a multitasking operation, i.e., simultaneous processing of a number of independent tasks with a single photonic processing unit, using both the space-and wavelength multiplexing techniques together. As a primitive demonstration, we also show the simultaneous processing of two independent tasks with respect to nonlinear channel equalization.
2 Speckle dynamics in multimode fibers

Multimode fiber speckles
Multimode fibers support hundreds of guided modes with different phase velocities. For monochromatic input light, a complex speckle pattern is generated in the intensity distribution at the end of the multimode fiber as a result of the interference of guided modes with different phase velocities [28,29]. A remarkable feature of a speckle pattern is its dependence on the wavelength of the input light. The longer the fiber length is, the more sensitive it becomes [29,30]. This suggests that when the wavelength (frequency) of the input light is modulated, the speckle pattern dynamically varies over time. In other words, information encoded by the wavelength modulation of the input light can be mapped as dynamically varying speckle patterns. See the Appendix for the theoretical model. Figure 1(a) shows the experimental setup used to generate and measure the modulation-induced speckle dynamics in a multimode fiber. We used a commercially available step-index multimode fiber with a core diameter of 50 µm, numerical aperture (NA) of 0.22, and length of 10 m. The multimode fiber supports approximately 250 guided modes at an input wavelength of λ 0 = 1550 nm. The multimode fiber was covered with an insulator to reduce the effects of thermal fluctuation and mechanical vibration. A narrow-linewidth tunable laser (Alnair Labs, TLG-200, 100 kHz linewidth, 10 mW) was used as a light source. The laser light was phase-modulated using a phase modulator (EO Space, AX-0MSS-20-LV, 16 GHz bandwidth) to encode an input signal, which was generated by an arbitrary waveform generator (Tektronix, AWG70002A, 25 GS/s), as an instantaneous frequency. The maximum phase modulation was approximately 0.9π. The modulated light was sent through a polarization-maintaining single-mode fiber to the multimode fiber. We used a lensed fiber (spot size of ∼ 2 µm with a working distance of ∼ 10 µm) as a probe to measure the time-variation of the output near the focal position at the end of the multimode fiber, and the output was sent to a photodetector (New Port, 1554-B, 12 GHz bandwidth) through an optical fiber amplifier (FiberLabs, AMP-FL80133-CB) and measured using a digital oscilloscope (Tektronix, DPO71604B, 50 GS/s, 16 GHz bandwidth).

Experimental setup
To obtain the time-variation of the speckle pattern, the measurement position was repeatedly changed by changing the position of the lensed fiber in parallel with the cross-section of the fiber end [ Fig. 1(b)] with a micro-positioning stage, whereas the input light was repeatedly modulated using the same signal. The origins of the time variation of the output signals at each measurement position were matched by adjusting each trigger signal to reconstruct the intensity patterns at each time. This method is effective because a multimode fiber is passive with respect to the input, and the output demonstrates consistency for the same inputs [31].

Speckle dynamics
Figure 2(a) shows an instance of the time variation of the intensity pattern measured by the above method, where the phase of the input light was modulated at 12.5 GS/s with a pseudorandom sequence and each intensity signal was measured along the x−axis shown in Fig. 1(b). The intensity pattern varies in a nonlinear manner, according to the input signal. The nonlinearity originates from the relationship Figure 1: (a) Experimental setup for generating speckle dynamics using a multimode fiber. ISO, optical isolator; Amp, electric amplifier; AWG, arbitrary waveform generator; PM, phase modulator; MMF, multimode fiber; OA, optical amplifier; PD, photodetector; OSC, digital oscilloscope. A porlarizationmaintaining single-mode fiber was coupled to the MMF through a standard FC/APC mating sleeve. The fiber probe (lensed fiber) was adjusted to measure the speckle at the end of the multimode fiber using a five-axis positioning stage. (b) Schematic of the cross-section at the end of the multimode fiber. The origin of the xy−axis was set as the center of the core layer.
between the input signal and output intensity signal, i.e., the input signal is encoded with the phase modulation, whereas the output is obtained as an intensity signal (see the Appendix). Considering that the computational capacity of the system depends on the number of linearly independent signals measured in the system [32], we measured the correlation between two signals measured at different positions, i.e., where φ λ0 (x, t) represents the intensity measured at position x at time t for input wavelength λ 0 . Here, f T = 1/T T 0 f dt denotes the mean of f over measurement time T . In addition,φ λ0 x and σ x are the mean and standard deviation of φ λ0 (x, t), respectively. The correlation C(x, x ′ ) decays for large ∆x = |x ′ − x|, as shown in Fig. 2(b). To clarify this, we calculated the mean correlation value C m (∆x) = 1/D |C(x, x + ∆x)|dx and plot C m (∆x) in Fig. 2(c). Here, C m (∆x) sufficiently decays when ∆x > 4 µm for the multimode fiber used in our experiment. We use the output signals sampled at the interval ∆x ≈ 4 µm for information processing, as discussed in Sec. 3.1.

Frequency dependence
The speckle dynamics depends on the wavelength (frequency) of the input light, and a variety of dynamics can be generated with different wavelengths in parallel. To confirm this, we varied the wavelength λ of the input light around λ 0 = 1550 nm and measured the intensity variation from the multimode fiber. In the results shown in Fig. 3(a), the input light is phase-modulated with the same pseudorandom sequence at 12.5 GS/s. As can be seen in the figure, various time variations are obtained, depending on the input wavelength λ. We also measured the correlation, where φ λ (x 0 , t) represents the intensity measured at x = x 0 and y ≈ 0 for input wavelength λ [ Fig. 3(b)]. Figure 3(c) shows the mean correlation C m (∆λ) = 1/D λ |C(λ, λ + ∆λ)|dλ. Here, C m (∆λ) sufficiently decays when ∆λ > 0.02 nm for the multimode fiber. ∆λ of 0.02 nm corresponds to lower bounds for the parallel generation of different speckles and modulation bandwidth for information encoding with input lights of different wavelengths. The correlation C m (∆λ) can more rapidly decay with an increase in the length of the multimode fiber [30], and a variety of dynamics that are more sensitive to the wavelength can be obtained. This result suggests that with multiple light sources of different wavelengths, the input information can be encoded even within the wavelength domain with a multiplexing technique, and parallel information processing is enabled, as demonstrated in Sec. 4. 3 Information processing based on speckle dynamics

Method
As demonstrated in the previous section, the use of a multimode fiber enables high-dimensional and nonlinear mapping of the input information in the speckle patterns. We use the mapping properties based on multimode fibers to process time-dependent input information in the following ways: The input light is phase-modulated with a signal, u(n) (n ∈ {1, 2, · · · }), where u(n) holds for the time interval τ . The output intensity signals, φ λ0 i (t) = |E(r i , t)| 2 , from the multimode fiber are measured at r i = (x i , y i ) with steps ∆x = ∆y ≈ 4 µm in the xy-plane [ Fig.1(b)] at the sampling time t = nτ , where i ∈ {1, 2, · · · , N } is the index of the measurement position. The readout is given by the following: where w i is a readout weight. The goal of the processing is to approximate a functional relationship between the input signal u(n) and target y(n) usingŷ(n). To this end, a finite set of training data {u(n), y(n)} Tn n=0 is utilized to calculate the optimal readout weights. This is performed through a simple ridge regression to minimize the error, Tn n |y(n) −ŷ(n)| 2 + λ N i w 2 i , where λ is the regularization parameter to avoid an ill-conditioned problem [12].

Chaotic time-series prediction
As a demonstration of a computation application using speckle dynamics, we conducted the Santa Fe time-series prediction task [33], which is a one-step ahead prediction task of chaotic data generated from a far-infrared laser [see Fig. 4(a)]. In this task, the input signal u(n) corresponds to the n-th sampling point of the chaotic waveform, and y(n) is set as the n+ 1-th sampling point, u(n+ 1). The prediction error was evaluated using the normalized mean square error (NMSE), which is given by 1/T n Tn n=1 |y(n)−ŷ(n)| 2 /σ 2 y , where σ 2 y is the variance of the target signal y(n). We used T n = 3000 steps for training and 1000 steps for testing. Figures 4(b) and 4(c) show the results of the measured intensity signals φ λ0 i (t) and readoutŷ(n) for an input rate of 1/τ = 12.5 GS/s of the signal u(n), respectively. According to Eq. (3), the readoutŷ(n) was calculated with N = 150 intensity signals, which were measured at 60×60 µm 2 near the fiber end. The NMSE of the prediction was approximately 0.09, which is comparable with the NMSEs obtained for the other photonic systems. (As an example, NMSEs of 0.106 and 0.109 have been reported for N = 388 and 124 in [13] and [16], respectively.) Unlike in other photonic information processing systems using laser devices [13,16], active devices are not involved in the processing in the proposed system; thus, the response time is not limited. We expect that faster speckle dynamics and then faster processing can be realized by increasing the input rate 1/τ . 1 (t), φ λ0 2 (t), · · · , φ λ0 N (t)], responding to the input light with wavelength λ 0 = 1550 nm, which is phase-modulated by the input signal u (see Visualization 2). These signals are used as nodes to obtain the readoutŷ [see Eq. (3)]. (c) The readout y(n) obtained using the measured signals, [φ λ0

Using wavelength-and time-multiplexing techniques
The number of nodes N required to construct a readoutŷ [Eq. (3)] is limited by the physical dimensions of the multimode fiber. However, as shown in Fig. 3(a), the signals generated for different input wavelengths are expected to be used as additional nodes in the proposed system, enabling processing on a larger-scale network. In addition, although the present experimental system does not have any memory effect, which is crucial for time-dependent signal processing, the memory effect can be introduced by adding an optical cavity in the system [22] or using the past information as a node in the time domain with a timemultiplexing method [19]. In this study, we simply use the latter method, i.e, the time-multiplexing method, and explore the availability of additional nodes in the time and wavelength domains to discuss the potential improvement of the computational performance. We here consider the following three readouts:ŷ where the output intensity φ λ0 i [(n − k)τ ] is sampled at a measurement position labeled by i for the fixed input wavelength λ 0 at time t = (n − k)τ . Here, φ λj 0 (t) is the output intensity measured at i = 0 for input wavelength λ j (j ∈ {1, 2, · · · M }) at time t. In addition,ŷ 1 (n),ŷ 2 (n), andŷ 3 (n) are the readouts in the space, wavelength, and mixed (space/wavelength) domains, respectively, and K denotes the number of past output signals used for calculating the readouts. The readout weights w i,k and w j,k were trained for the Santa Fe time-series prediction task. Figure 5 shows the result of the prediction task, which was obtained by the trained readoutsŷ l (n) (l ∈ {1, 2, 3}), where λ 0 = 1550 nm and λ j = 1549.5 + 0.02j nm. For comparison, the total number of nodes, N + M , is fixed at 50 for all cases. As K increases, the NMSEs decrease and the prediction performance improves. The readoutŷ 1 (n) shows the best NMSE of 0.06 for K = 7, which is better than that of other photonic RC [13,16]. The readoutsŷ 2 (n) andŷ 3 (n) exhibit a relatively improved performance, and the best NMSE was 0.052. This improvement may be attributed to the linear independence among the signals φ λj 0 in the wavelength domain [ Fig. 3(a)]. These results suggest that the simultaneous use of the nodes in the space, wavelength, and time domains is effective for processing time-dependent signals.

Multitasking
Finally, we discuss the simultaneous use of the above space and wavelength multiplexing techniques, which also allows a multitasking operation, i.e., parallel processing of numerous independent tasks, with a single photonic component, that is, a multimode fiber [ Fig. 6(a)]. This is enabled by the encoding of different input information at different wavelengths because there are no nonlinear interactions between the lights with different wavelengths. Figure 6(b) shows the experimental setup for the multitasking operation. We used two light sources (wavelength-tunable lasers), the wavelengths of which are λ 1 and λ 2 , respectively. The respective lights were then independently phase-modulated with different input signals, u 1 and u 2 , and injected into a multimode fiber. The output light of λ 1(2) from the multimode fiber was detected by PD1(2) after passing through a tunable optical filter (Filter 1(2)) with a bandwidth (full width at half maximum) of ∆ ≈ 1.1 nm, the center transmission wavelength of which was tuned at λ 1 (2) and used as the spatial nodes for calculating the readouts. Figure 6: (a) Schematic of multitasking based on passive photonics. Each input signal u j (j ∈ {1, 2, · · · , M }) is encoded as optical signals with wavelength λ j and is independently processed in a single photonic component (multimode fiber) to obtain the readoutỹ i . (b) Experimental setup used for multitasking operation based on a multimode fiber. Two input signals are generated from an arbitrary waveform generator (AWG). The input lights with wavelengths λ 1 and λ 2 are phase-modulated using PM1 and PM2, respectively, and injected into a multimode fiber (MMF). The outputs from the MMF are separated by two optical filters (Filter 1 and 2) and sent to photodetectors PD1 and PD2 through optical amplifiers (OA1 and OA2), respectively.
As a demonstration, we used the nonlinear channel equalization task [2]. A goal of this task is to reconstruct four digital signals {−3, −1, 1, 3} transmitted through a communication channel with nonlinear distortion. The nonlinear transformation of the communication channel is given by a model equation [20] q(n) = 0.08d(n + 2) − 0.12d(n + 1) + d(n) + 0.18d(n − 1) where d(n) ∈ {−3, −1, +1, +3} is the original signal before transmission through the communication channel, q(n) is the linear channel output, u(n) is the noisy nonlinear channel output, and v(n) is white Gaussian noise with a zero mean. In our experiment, we assumed two original signals, d 1 (n) and d 2 (n), consisting of two independent random sequences, and two different noisy communication channels with different noises, v 1 (n) and v 2 (n). The power ratio between the signal d i (n) (i = 1, 2) and noise v i (n) (i = 1, 2) was set to 30 dB. The goal of our experiment is to simultaneously recover the two original signals, d 1 (n) and d 2 (n), from the respective noisy nonlinear channel outputs u 1 (n) and u 2 (n). The computation performance was evaluated using the symbol error rate (SER). Figure 7 shows the parallel processing at a signal input rate of 12.5 GS/s for this task, where the respective readouts for channels 1 and 2 were calculated asỹ 1(2) [(n − k)τ ] with N = 100 and K = 7 in this experiment. We used T n = 3000 samples for training and 3000 samples for testing. When the difference between the two wavelengths, ∆λ = |λ 1 − λ 2 |, is smaller than the bandwidth of the optical filter ∆ ≈ 1.1 nm, the optical signal of the wavelength λ 2(1) cannot be sufficiently eliminated with the optical filter 1 (2), and their interference causes large errors in the parallel computing [ Fig. 7(a)]. The SERs were worse than that for the single operation, i.e., SER 0 = 0.022. However, when ∆λ > ∆, the SERs improved, and SER/SER 0 ≈ 1.0 was obtained, i.e., parallel processing can be performed, although the original error rate, SER 0 , was worse than in previous studies. The low SER 0 may be due to the low stability of the multimode fiber in our experiment; it is likely to improve when the multimode fiber is replaced with a multimode waveguide in a photonic chip. We consider that further parallel and multitasking processing will be conducted using optical filters with a narrower bandwidth. The number of parallel operations will not be limited in principle. This is a unique property of the proposed optical information processing system.  (2)] normalized by SER 0 as a function of the difference in wavelength |λ 1 − λ 2 |. Here, SER 0 denotes the error rate for the case of a single operation. In (a)-(c), the input signals were set to include two white Gaussian noises v 1 and v 2 with zero mean adjusted in power to yield an SNR of 30 dB.

Summary and discussion
In this study, we experimentally demonstrated that the use of multidimensional speckle dynamics in multimode fibers enables fast and parallel information processing for time-dependent signals at input rates of 12.5 GS/s. The rate of information processing is limited only by the bandwidths of the phase modulator and photodetectors in the present experiment, and can increase further when those with larger bandwidths are used. Larger-scale information processing can in principle be achieved through the simultaneous use of the output signals (network nodes) in both the space and wavelength domains. In addition, the use of wavelength multiplexing combined with spatial speckles enables a multitasking operation, as demonstrated in Sec. 4.
Although the present experimental system does not have any memory, the memory effects can be easily introduced by adding an optical cavity to store past information, and the system with a memory can be used as a reservoir computer. Figure 8 shows an example of the entire photonic architecture for reservoir computing. The multimode waveguide with a spiral geometry or interferometer structure fabricated on a silicon chip can provide a long optical path length and high NA to generate speckles sensitive to the wavelength [34,35,36]; further, it can be used to induce speckle dynamics with a shorter latency in a small footprint. Each signal from the multimode waveguide can be split by a splitter and wavelength demultiplexer and detected by the photodetector array. This type of photonic integration will offer a pathway for compact, parallel, large-scale photonic computers. Appendix: Speckle-based mapping using a multimode fiber Here, we show a mapping model based on dynamic speckle patterns in a multimode fiber for a timedependent signal u(t). In the experimental setup shown in Fig. 1(a), the laser light is phase-modulated using a voltage V (t) proportional to u(t) and sent to a multimode fiber. The input electric field is given by E in (r, t) = E 0 (r)e i(πV (t)/Vπ+ω0t) , where V π is the voltage required to produce a π-phase shift in a phase modulator, r = (x, y) is the transverse dimension of the propagation, and ω 0 is the input angular frequency. Given that the field can be expressed as the superposition of signals with angular frequency ω, i.e., E in (r, t) = Ê in (r, ω)e iωt dω, each frequency component after propagating in a multimode fiber is written as follows [30]:Ê (r, z, ω) = m A m (ω)Ψ m (r, ω)e −iβm(ω)z , where A m (ω) is the complex amplitude of the m-th guided mode which has the spatial profile Ψ m (r, ω) and propagation constant β m (ω) for input angular frequency ω, and is determined such that the boundary conditionÊ(r, 0, ω) =Ê in (r, ω) is satisfied. Thus, the entire output field at the end of a multimode fiber with length L at time t can be expressed by the Fourier transform: E(r, L, t) = Ê (r, L, ω)e iωt dω, and the intensity of the field at r i (i = 1, 2, · · · , N ) is measured as φ(r i , t) = |E(r i , L, t)| 2 . As shown in the above model, the mapping from time-dependent signal u(t) to signal φ(r i , t) is nonlinear and exhibits high dimensionality when A m (ω) = 0 for a wide range of frequencies ω and modes m, i.e., the input signal is modulated at fast rates, and many guided modes are excited by the input field.