Online training for high-performance analogue readout layers in photonic reservoir computers

A reservoir computer, schematised in figure 1, consists of a recurrent network of internal variables, usually called “nodes” or “neurons”, from its biological origins lukovsevivcius2009survey . These $N$ variables, denoted by $x_{i}(n)$, with $i=0,\ldots,N-1$, evolve in discrete time $n\in\mathbb{Z}$, as follows

where $f$ is a nonlinear function, $u(n)$ is the input signal that is injected into the system, $\alpha$ and $\beta$ are feedback and input gains, respectively, used to adjust the dynamics of the system, and $M_{i}$ is the input mask, drawn from a uniform distribution over the interval $[-1,+1]$, as in e.g. rodan2011minimum ; paquot2012optoelectronic ; duport2012all . In our implementation, we use a sine function $f=\sin(x)$ as nonlinearity and a ring topology appeltant2011information ; rodan2011minimum to simplify the interconnection matrix of the network, so that only the first neighbour nodes are connected. Both choices are dictated by the proposed hardware setup of the opto-electronic reservoir. As will be discussed in Sec. 2.4, we use a light intensity modulator with a sine transfer function as the nonlinear node, and a delay system with time-multiplexing of the reservoir states. A detailed discussion of these experimental aspects can be found in the Supplementary Material of paquot2012optoelectronic .

The reservoir computer produces an output signal $y(n)$, given by a linear combination of the states of its internal variables

where $w_{i}(n)$ are the readout weights, trained either offline (using standard linear regression methods), in which case they are time independent, or online, as described in the next section, in order to minimise a task-dependent error function, which will be introduced alongside the benchmark tasks further in this paper. As discussed below in detail, in hardware implementations of reservoir computing the readout layer can be more complex than Eq. (2).

Contrary to offline, or batch learning, where the entire training dataset is used at once to compute the best readout weights $w_{i}$, online training approach handles the data sequentially in order to optimise the performance step by step. As discussed in Sec. 1, this allows the Reservoir Computer to be optimised without accurate knowledge of the underlying hardware, which is exactly what is required for an analogue readout layer. This approach can be realised with various algorithms, and in this work we chose to work with the simple gradient descent algorithm, for ease of implementation on the FPGA board.

The gradient, or steepest, descent method is an algorithm for finding a local minimum of a function using its gradient arfken1985mathematical . For the task considered here (see the next section) we update the readout weights using the procedure given in bishop2006pattern

where $\lambda$ is the step size, used to control the learning rate, and $d(n)$ is the task-dependent target signal (see sections 2.3.1 and 2.3.2). The origin of this procedure is that if the error at time $n$ is given by $(d(n)-y(n))^{2}$, then the derivative of the error with respect to $w_{i}$ gives $(d(n)-y(n))x_{i}(n)$, i.e. the right-hand side of Eq. (3). At high values of $\lambda$, the weights get close to the optimal values very quickly (in a few steps), but keep oscillating around these values. At low values, the weights converge slowly to the optimal values. In practice, we start with a high value $\lambda=\lambda_{0}$, and then gradually decrease it during the training phase until a minimum value $\lambda_{min}$ is reached, according to the equation

with $\lambda(0)=\lambda_{0}$ and $m=\left\lfloor n/k\right\rfloor$, where $\gamma<1$ is the decay rate and $k$ is the update rate for the parameter $\lambda$. Previous work has shown that setting $\lambda_{0}=0.4$, $\lambda_{min}=0$ and $\gamma=0.999$ is a reasonable choice for good performance antonik2016online . The update rate $k$ defines the convergence speed and the resulting error rate : higher $k$ requires a longer training dataset but offers better results. More complex optimisation methods could be used here, such as simulated annealing press1986numerical or stochastic gradient descent bottou2004stochastic . However, the above technique is very simple and provides sufficiently good results for this application.

We tested our system on two benchmark tasks commonly used by the reservoir computing community: wireless channel equalisation and emulation of a 10-th order Nonlinear Auto Regressive Moving Average system (NARMA10). Note that, to the best of our knowledge, the latter has never been tested on an online-trained reservoir computer.

Fig. 2 depicts the proposed experimental setup that we have investigated using numerical simulations. This section overviews the three main components: the optoelectronic reservoir, the analogue readout layer and the FPGA board.

All numerical experiments were performed in Matlab. We used a custom model of our reservoir computer, based on previous investigations paquot2012optoelectronic ; antonik2016online , that has been shown to emulate very well the dynamics of the real system. The simulations were performed in discrete time, and took into account the internal structure of the Reservoir Computer described above, such as the ring-like topology, sine nonlinearity and the analogue readout layer with an RC filter. The simulations allow to try out different configurations and to scan various experimental parameters, including values that are impossible to achieve experimentally or imposed by the hardware. All simulations were performed on a dedicated high-performance workstation with 12-core CPU and $64$ Gb RAM. Since the convergence of the gradient descent algorithm is quite slow, we limited our investigations to a fast update rate $k=10$ (see Eq. (4)), so that each simulation lasted about 24 hours.

The principal goal of the simulations was to check how th online learning approach would cope with experimental difficulties encountered in previous works smerieri2012analog ; duport2016fully . To that end, we gathered a list of possible issues and scanned the corresponding experimental parameters in order to check the system performance. In particular, we investigated the following parameters:

• •

  The RC integrator ratio $\rho$. This is the most important parameter of the analogue readout layer. While its accurate measure is not required in our setup – since we do not correct the readout weights $w_{i}$ – it defines the integration span of the filter, and thus the reservoir states that contribute to the output signal. It can thus significantly influence the results. Another question of importance is how dependent the system performance is on the exact value of $\rho$.

• •

  The MZ modulator bias. Mach-Zehnder modulators need to be applied a constant voltage to maintain their transfer function symmetric around zero. The devices we were using up to now are subject to slight drifts over time, often resulting in a non-perfectly symmetric response. We thus checked in simulations whether such an offset would impact the results.

• •

  The DAC resolution. The precision of the DACs on the FMC151 daughtercard is limited to 16 bits. Numerical investigations have shown that the precision of readout weights has a significant impact on the performance, see e.g. soriano2013optoelectronic ; soriano2015delay ; antonik2016towards . We thus checked whether the resolution available is enough for this experiment.

Besides these potentially problematic parameters, we also scanned the input and feedback gain parameters (denoted by $\beta$ and $\alpha$ in Eq. (1)) in order to find the optimal dynamics of the reservoir for each task.

In a separate set of simulations, we investigated the applicability of the proposed method to nonlinear readout layers. That is, we checked whether the simple gradient descent method would still work with a nonlinear response of the analogue readout layer with respect to the reservoir states $x_{i}(n)$ (see Eq. (12)). We picked two “saturation” functions of sigmoid shape. This choice arises from the transfer function of common light detectors that are linear at lower intensities and become nonlinear at higher intensities. We used the following functions: a logistic function, given by

and the hyperbolic tangent function, given by

These functions, $g_{\text{lg}}$ and $g_{\text{ht}}$, do not model any particular photodiode, but are two simple examples that allow us to address the above question. Both functions are plotted in figure 5, together with a linear response, for comparison.

We investigated two possible nonlinearities in the output layer. In the first case, the readout photodiode ($\text{P}_{\text{r}}$ in figure 2) produces a nonlinear response, while the balanced photodiode ($\text{P}_{\text{b}}$ in figure 2) remains linear. This scenario, that we shall refer to as “nonlinear readout”, allows one to investigate what happens when the reservoir states $x_{i}$ used to compute the output signal $y(n)$ (see Eq. (2)) differ from those employed to update the readout weights (see Eq. (3)). Thus in this case the update rule (Eq. (3)) for the output weights becomes

where $g$ is given by either Eq. (13) or Eq. (14), while the output layer is given by Eq. (12).

In the second case, called “nonlinear output”, the readout photodiode is linear, but the balanced photodiode exhibits a saturable behaviour. In this case the update rule Eq. (3) for the output weights is unchanged, but the output layer Eq. (12) becomes

Note that we have only considered cases with just one nonlinear photodiode, so as to check whether the difference between the reservoir states used for training and those to compute the readout (see Eqs. (3) and (2), respectively) would degrade the performance of the system. The scenario with both nonlinear photodiodes is hence more simple, as the reservoir states are the same in both equations. One could consider the case with two photodiodes exhibiting different nonlinear behaviours. In that situation, similar to the results we will show in Sec. 3, we expect the algorithm to cope with the difference up to a certain point, before running into troubles. For this reason, we leave that scenario for future investigations.

In this section we present our numerical results and answer the questions raised in the previous section.

For each of the two tasks considered here, we performed three kinds of simulations: we scanned the RC integrator ratio $\rho=\theta/\tau$ in the first simulation, the MZ bias in the second, and the resolution of the DAC in the third. Furthermore, since different values of these parameters may work better with different dynamics of the reservoir, we also scanned the input gain $\beta$ and the feedback gain $\alpha$ in all three simulations independently, and applied the optimal values in each case.

For both tasks, we used a network with $N=50$ neurons, as in most previous experimental works paquot2012optoelectronic ; duport2012all ; vinckier2015high ; duport2016fully . The reservoir was trained on 83000 inputs, with an update rate $k=10$, and then tested over $10^{5}$ symbols for the channel equalisation task and $10^{4}$ inputs for NARMA10 task. For statistical purposes, we ran each simulation 10 times, with different random input masks. In the following figures, averaged results over the masks are plotted, while the error bars give the standard deviation over the different input masks. Results related to the channel equalisation task are plotted with solid lines, while dashed lines correspond to those for NARMA10.

For the channel equalisation task, our system yields SERs between $10^{-4}$ and $10^{-3}$ depending on the input mask, as summarised in Tab. 1 (first line). This is comparable to previous experiments with the same opto-electronic reservoir: error rates of order of $10^{-4}$ were reported in paquot2012optoelectronic using a digital readout and in duport2016fully with an analogue readout, using an RLC filter. The first experimental analogue system, using a simple RC circuit, as we did in this work, performed significantly worse, with SER of order of $10^{-2}$ smerieri2012analog . That is, online learning does not outperform other methods, but allows to obtain significantly better results with a simpler setup.

As for the NARMA10 task, we obtain a NMSE of $0.20\pm 0.02$. Previous experiments with a digital readout layer produced $0.168\pm 0.015$ paquot2012optoelectronic and $0.107\pm 0.012$ vinckier2015high . With an analogue readout layer, the best NMSE reported was $0.230\pm 0.023$ duport2016fully . Our system thus slightly outperforms the analogue approach, and gets close to the digital one, except for the very good result obtained with a different reservoir, based on a passive cavity vinckier2015high . Again, our results were obtained with a simple setup and no modelling of the readout, contrary to duport2016fully .

Furthermore, the error rates obtained here can be significantly lowered with more training, as has been demonstrated numerically and experimentally in antonik2016online . To keep reasonable simulation times (about 24 hours per simulation), we limited the training to 83000 input values, with an update rate $k=10$. Higher update rates can be used experimentally, since running the opto-electronic setup is much faster than simulating it. We thus expect to obtain lower error rates experimentally with longer training sets and update rates up to $k=200$. To illustrate this point with results reported in antonik2016online , short training sets with $k=10$ yielded SERs of order of $10^{-4}$ for the channel equalisation task. Increasing $k$ up to $200$ allowed to decrease the error rate down to $5.7\times 10^{-6}$.

Figures 3 and 3 show the influence of input and feedback gain parameters on the performance of the system. All curves present a pronounced minimum, except for the input gain $\beta$ for the NARMA10 task, where values above 0.6 seem to produce comparable results. Note that the channel equalisation task requires a low input signal with $\beta=0.2$, while NARMA10 works best with stronger input and $\beta=0.8$. As for the feedback gain, NARMA10 shifts the system close to the chaotic regime with $\alpha=0.95$, while channel equalisation works better with $\alpha=0.8$.

Figure 4 shows the results of the scan of the RC integrator ratio $\rho$. Both tasks work well on a relatively wide range of values, with NARMA10 much less sensitive to $\rho$ than channel equalisation. In particular, the channel is equalised best with $\rho=0.03$. With $N=50$, this corresponds to $\tau=T/0.03N=5.29\;\text{\textmu s}$, which is shorter than the roundtrip time $T=7.94\;\text{\textmu s}$. On the other hand, NARMA10 output is best reproduced with $\rho=0.003$, which yields $\tau=T/0.003N=52.93\;\text{\textmu s}$. This is significantly longer than the roundtrip time $T$, meaning that reservoir states from previous time steps are also taken into account for computation of an output value. This is not surprising, since NARMA10 function has a long memory (see Eq. (7)). However, this memory effect in the readout layer is not crucial, as the system performs equally well with higher $\rho$ and thus lower $\tau$. All in all, these results are very encouraging for upcoming experiments, as they show that an accurate choice of capacitor is not crucial for the performance of the system.

Figure 4 illustrates the impact of the bias of the readout Mach-Zehnder modulator on the reservoir computer performance. NARMA10 task is clearly more affected by this offset, as the NMSE grows quickly from a bias of roughly $0.06$. The SER, on the other hand, stays low until $0.1$. For a MZ modulator with $V_{\pi}=4.5\;\text{V}$ this corresponds to a tolerance of roughly $0.1\;\text{V}$, which is superior to expected experimental deviations. The Hameg power supply that we use to bias the modulator (see Sec. 2.4) has a resolution of $0.001\;\text{V}$.

Figure 5 shows that the 16-bit DAC resolution is not an issue for this experiment, as the minimal precision required for good performance is 8 bits, for both tasks.

Table 1 sums up the results obtained with a nonlinear readout layer. We used optimal experimental parameters, as described above, and generated new sets of data for the training and test phases. We investigated two scenarios and used two functions of sigmoid shape, $x\rightarrow g_{\text{lg}}(x)$ and $x\rightarrow g_{\text{ht}}(x)$, as described in section 2.5. The system was trained over 83000 inputs, with an update rate $k=10$, and tested over $10^{5}$ symbols for the channel equalisation task and $10^{4}$ inputs for NARMA10. We report error values averaged over 10 trials with different random input masks, as well as the standard deviations. The figures show that the performance deterioration is more manifest with the hyperbolic tangent function $g_{\text{ht}}$, as it is much more nonlinear than the logistic function $g_{\text{lg}}$. Overall, the added nonlinearity does not have a significant influence on the results in both cases. The SER roughly doubles, at most, for the channel equalisation task. The impact on NARMA10 is barely noticeable, as the error increase of $5\%$ is smaller than the standard deviation. Using offline training on the same system (i.e. with nonlinear output) we observed an increase of the SER by one order of magnitude for the channel equalisation task, and a 30% increase of the NMSE with the NARMA task. These results show that online training is very well suited for experimental analogue layers, as it can cope with realistic components that do not have a perfectly linear response.