Reservoir Computing with 3D Nanowire Networks

The most basic dynamical element within a NWN is the junction that results from two wires making contact. We model the internal dynamics of the memristive junction as a conductive ‘hillock’ (Figure 1b) that grows through the insulating layer between wires under the influence of an electric field, creating a non-equilibrium structure. Surface energy effects will attempt to reduce the size of this protrusion (Wang et al., 2019). The dynamical equation governing the growth of the hillock is then

where $z$ is the height of the hillock and $D$ is the distance between wire cores as in Figure 1b. $V$ is the voltage across the junction, and $\mu$ and $\kappa$ are the parameters which govern growth and relaxation respectively ($\mu=0.346$nm²V^-1 and $\kappa=0.038$s^-1 – see discussion below). Eq (1) can then be numerically integrated using the Euler method to obtain a discretized change in hillock height at each computational timestep. The conductance of each junction is governed by quantum mechanical tunneling, and so is a function of the instantaneous size of each tunnel gap, $D-z$:

Hence each junction has a nonlinear response to a voltage input. Note that the response of the network results from the interaction of all the junctions.

Using Eq (1), we can identify a critical voltage $V_{c}$. For $V<V_{c}$, the hillock grows towards a stable equilibrium value. For $V>V_{c}$, the hillock will grow until $z=D$, forming a conductive bridge between the two wires (Figure 1c). We focus here on the low voltage regime where no hillock completely bridges a tunnel gap.

We now show that the junction model given in Eq (1) is equivalent to that of an ideal voltage-driven memristor (Strukov et al., 2008) parameterized using an internal variable $w$ and governed by the dynamical equation

with the resistance across the gap of size $D$, given by

where $R_{ON}$ and $R_{OFF}$ are the resistance across the memristor in the high conductive state and the low conductive state respectively. If we assume (Caravelli, 2019) a fixed ratio of $r=R_{OFF}/R_{ON}$, we can combine Eq (3) and (4) to get

where $\chi=(R_{OFF}-R_{ON})/R_{OFF}$. The ratio $r$ is typically on the order of $10^{3}-10^{5}$ (Kuncic and Nakayama, 2021), so when $r$ is large, $\chi\to 1$. The resulting equation is then

When compared with Eq (1), we can see that they are equivalent up to a factor of $r$ in the first term, which has the effect of renormalizing the relevant range of $\mu$ (note that $r$ is a parameter of the memristor model, but not our simulations). Hence the hillock growth model is equivalent to models of diffusive unipolar memristive junctions operating in a low voltage regime, and the results presented in this paper are applicable to the more general case of nanowire networks with memristive junctions.

We note that other authors have considered alternate regimes and junction models, including bipolar switching models. For example Hochstetter et al. (2021) and Loeffler et al. (2021) consider a regime where each hillock can extend close to the opposite side of the junction gap, leading to a dramatic increase in conductance which then redistributes the voltage across other junctions. In this regime, conductive bridges can continuously evolve throughout the network, leading to more complex dynamical behaviour. We emphasize that we consider the low voltage regime where the hillocks never completely bridge the junction gaps.

We constructed two different types of NWNs, 2D and Q3D. The $x-$ and $y-$coordinates of the centres of the wires were chosen randomly from a uniform distribution in the range $[0,\Lambda]$, where $\Lambda$ is the width of the square deposition area. The orientation of the wires with respect to the $x$-axis were drawn randomly from a uniform distribution in the range $[-\pi/2,\pi/2]$. Every point where two wires intercept is considered to be a junction in the network. The connectivity is then stored in an adjacency matrix.

In the 2D networks, the wires are rigid 1D lines, and hence all lie within the same plane. However, in the Q3D networks, the wires are rigid 3D volumes in space. The positioning of a new wire therefore depends upon those that have been previously deposited. For a new wire, the algorithm determines all potential intercept points, and in the general case, the point of intercept with the highest $z$ value acts as the first contact point and as a pivot around which the new wire must rotate. The center of mass of the new wire then determines the direction of rotation and hence the second contact point. This results in a stacking of the wires in the vertical direction, leading to very different topological properties of the 2D and Q3D networks (for more details, see Daniels and Brown (2021)).

The input and output electrodes are defined by placing additional nanowires in the deposition area to act as contacts (they are placed before the deposition of the NWN). The input electrode is deposited on the left edge of the deposition area, and the output electrodes are evenly distributed along the top, bottom, and right edges (as in Figure 1a). The ‘contact’ wires are then treated in the same way as all other wires, with additional rows and columns added to the adjacency matrix representation of the network. The input signal is a voltage that is applied to the input electrode, and the output electrodes are grounded. We then solve Kirchhoff’s circuit laws for the network, and the current is recorded from each of the output electrodes. All simulations were implemented in Python v3.8.5.

The memory capacity task aims to reconstruct a delayed version of the input signal from the measured outputs (Jaeger, 2002a). In other words, the task determines how precisely the reservoir can reconstruct the input of previous time-steps.

Figure 2a shows a schematic of the MC task. The input to the network, $\mathbf{U}$, is a random sequence of voltage values drawn uniformly from the range $[0,V_{max}]$. The input is applied at the single input electrode and the output signal is collected from the readouts (which can be either the electrodes, as illustrated in Figure 2a, or every wire). The measured readout sequences are then combined into a predictor matrix, $\mathbf{X}$. $MC_{k}$ quantifies the performance of the network in reconstructing a version of the input that is delayed by a time $k$. In other words the target $\mathbf{{y}}$ at time $t$ is the same as the input at time $t-k$. The reconstructed signal $\mathbf{\hat{y}}$ (i.e. the prediction of the reservoir) is obtained from a linear combination of the readouts at time $t$. For the delay $k$ the required weight vectors, $\mathbf{w}_{k}=\textbf{X}^{\dagger}\mathbf{y}$, are obtained from the Moore-Penrose Pseudo inverse $\mathbf{X}^{\dagger}$, which minimizes the root mean square error (RMSE)

Here $\mathbf{\hat{y}}=\mathbf{w}_{k}\mathbf{X}$ and $L$ is the length of the input sequence. After training, the memory capacity for each $k$ is obtained from the predicted output for a different test input sequence, and is defined as the squared correlation coefficient between the target and predicted signals,

The total memory capacity is then a sum over all possible delays

The NLT task demonstrates the ability of the network to perform a non-linear mapping using the higher harmonics generated by the internal dynamics of the network (Sillin et al., 2013; Jaeger, 2002c). As illustrated schematically in Figure 2, an input sine wave $\mathbf{U}$ with signal level between $V_{min}$ and $V_{max}$ is fed into the network at the input electrode. The process of training and testing the network is essentially identical to that of the MC task, except the target function $\mathbf{y}$ is a square wave in the range $[-1,1]$ (Figure 2b). The normalized root mean square error (NRMSE) between the target function and the transformed function $\mathbf{\hat{y}}$,

is used as the performance metric.

In many machine learning models, adding noise to the input signals acts as a type of regularization (Bishop, 1995) - i.e. it prevents the weights from taking extreme values (overtraining). The result is a model that generalizes better and is more resilient to noisy test data (Bishop, 2006; Goodfellow et al., 2016). Therefore, in order to investigate the generalizability of NWN to performing the nonlinear transformation, a degree of random noise is added to the input training signals. The networks are then tested using the regression weights obtained from training.

Figure 3a shows the impact of the number of electrodes on the memory capacity in E-mode $MC_{E}$ and in N-mode $MC_{N}$. Perhaps surprisingly, given the significant differences in network structure (Daniels and Brown, 2021), $MC_{E}$ for the 2D (red) and Q3D (blue) networks are strikingly similar. The same is also true of $MC_{N}$ (yellow for 2D, green for Q3D). In both types of network, as $E$ is increased, $MC_{E}$ increases gradually but gains in performance for $E>12$ become small. $MC_{N}$ is not affected by the number of electrodes, and is higher than $MC_{E}$, consistent with the expectation that readout from all nodes should provide an upper limit on performance.

In order to determine the impact of $N$ on the memory capacity, we set $E=$ 24 and varied the number of wires. Figure 3b shows that $MC_{E}$ is always similar for the 2D and Q3D networks (compare red and blue curves) and the difference becomes smaller as $N$ increases. This is surprising, since as $N$ is increased the stacking in the Q3D network becomes more pronounced. $MC_{N}$ for the 2D NWNs begins to drop off gradually and becomes comparable with that of $MC_{E}$, whereas $MC_{N}$ for the Q3D NWNs remains relatively constant. Hence, for large $N$ the upper limit on performance of the Q3D networks is significantly higher than for the 2D networks. As $N$ increases, more wires become shorted together in the 2D networks (decreasing the average degree), and this effectively reduces the number of wires with unique outputs and the results become similar to those obtained from using a small sample of the wires (and the performance using electrode currents). Recent work by Han et al. (2021) has also shown that the MC of directed acyclic networks decreases rapidly as the density of nodes in the network increases. These results are also consistent with the findings of Loeffler et al. (2021) and Zhu et al. (2020b), which show MC declining with an increase in the mean degree of the networks.

Figure 3c shows the impact of the input voltage on the memory capacity for a fixed $E=$ 24. The 2D and Q3D networks were subjected to the same voltage input sequence, but for each realization of the networks, the range of $V_{max}$ was constrained by the requirement that all junctions remain in the low voltage regime (i.e. that no hillock formed a complete bridge across any junction). Figure 3c shows the change in MC with $V_{max}$ (the values of $V_{max}$ on the $x$-axis are min-max normalized which allows us to compare performance relative to the critical voltage, $V_{C}$). Initially, an increase in input voltage causes a rapid increase in both $MC_{N}$ and $MC_{E}$, as a greater degree of nonlinearity in the junctions can be exploited. However, as $V_{max}$ increases further there is a saturation and then a slight decline. This effect is also seen in directed acyclic networks: as the input signal is increased, MC decreases (Han et al., 2021).

The inset of Figure 3c shows $MC_{E}$ of each realization of the system as a function of the absolute voltage. Nearly twice the input voltage can be applied to the Q3D networks compared to the 2D networks. This is consistent with previous work (Daniels and Brown, 2021) showing that the Q3D networks have longer path lengths – as voltage is distributed across many junctions in series, the longer path between electrodes ensures a lower voltage across each junction.

Figure 4 compares the performance in the NLT for different network types. The input sequence used was a 1V sinusoidal signal of 1Hz frequency for 21 periods. The network contains $N=500$ nanowires. In Figure 4a the red and blue curves show the performance in E-mode (NRMSE_E) for the 2D and Q3D networks respectively. As $E$ is increased, NRMSE_E drops (i.e. performance improves), with the 2D networks slightly outperforming the Q3D networks. Figure 4 also shows the theoretical maximum performance, NRMSE_N, of the networks in $N$-mode for the 2D and Q3D networks (yellow and green respectively). NRMSE_E approaches NRMSE_N as $E$ is increased, and for $E>12$ there is very little performance gain, as was also observed for the MC task. When using N-mode to perform the NLT, the testing and training NRMSE_N are similar to one another for the Q3D networks, yet the testing NRMSE_N is significantly higher than the training NRMSE_N, suggesting that the 2D network weights are overfitting the training data (not shown).

Figure 4b shows the impact of $N$ on the nonlinear tranformation performance. As with the MC task, $E$ is fixed at 24, and $N$ ranges from 400 to 1500. In E-mode, performance improves slightly with increasing $N$, until $N\sim$1000. When using N-mode, performance remains constant for $N>500$ wires, with very little variance across different realizations in both networks.

Figure 4c shows that as $V_{max}$ increases, there is a gradual improvement in performance, which then plateaus at 3-4V for both networks. There is very little difference in performance between the 2D and Q3D networks. This is once again surprising given the different topological structure of the 2D and Q3D networks. However, as shown in the inset of Figure 4c, the Q3D network can be subjected to nearly twice the input amplitude compared to the 2D networks, achieving consistent NRMSE scores over a wider range of voltage inputs.

In order to successfully transform the input signal into the target, the network must be capable of higher harmonic generation (HHG). We therefore also examined the higher harmonics produced by each network.

Figure 5 shows the HHG for an example 2D and an example Q3D network both in a 24 electrode configuration, showing a highly nonlinear network response to a sinusoidal input, and the generation of a wide range of higher harmomics. The generation of the higher harmonics can be largely attributed to the nonlinear response of the network. Hence the network is able to perform a nonlinear mapping of the input features to a higher-dimensional space. Note that in Figure 5, we can see that the Q3D network has a smaller degree of hysteresis. This is likely due to the lower voltage across each junction caused by the shorter mean path lengths within the Q3D network.

In order to test the resilience of the network performance, we trained the output weights for each of the networks in the same manner as for the NLT task, but added input noise ranging from 0 to 5% of the input signal amplitude. We then exposed the trained networks to testing inputs with different amounts of noise.

When the networks receive training input with small amounts of noise ($<1\%$), the Q3D networks achieve a lower NRMSE on the testing signal than the 2D networks. In fact, when the networks are trained in the absence of noise, the test performance for the 2D network is poor when even the slightest amount of noise is added to the test input (see the nearly vertical red line close to the vertical axis of Figure 6a). When the 2D networks are trained on input signals with greater amounts of noise, the testing performance significantly improves.

In general, the 2D networks have a slightly better performance (lower NRMSE) than the Q3D networks. However, the Q3D networks show a significant jump in performance when trained on input with noise above $\sim$2%. The average performance (bold dashed lines) of the two networks begins to converge for higer levels of noise. Figure 6b shows the evolution of the weights as more input noise is added. When trained with low noise, the distribution of weights in both networks is narrow. As the noise on the training input is increased, the weight distribution for the 2D networks remains stable, but the weights for the Q3D network undergo a sharp transition and become much more widely distributed.

Figure 6(c-f) show examples of NLT performance for a 2D and Q3D network exposed to a low amount of input noise (0.25%). The 2D network achieves a better training performance, but the Q3D network achieves a significantly better test NRMSE. The high values of the weights when trained on low noise and the comparatively poor performance of the networks on the testing data suggest that the regression of the output from the 2D network is overfitting. Adding noise acts to regularize the output weights, allowing it to generalize better on the more noisy test data, as seen by the flatter curves in Figure 6a for the Q3D networks. The slightly narrower range of weights for the Q3D network trained on low noise avoids over fitting and therefore makes them more resilient.