materials Article Self-Organizing Neural Networks Based on OxRAM Devices under a Fully Unsupervised Training Scheme Marta Pedró 1, * , Javier Martín-Martínez 1 , Marcos Maestro-Izquierdo 2 , Rosana Rodríguez 1 and Montserrat Nafría 1 1 Departament d’Enginyeria Electrònica, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain; javier.martin.martinez@uab.es (J.M.-M.); rosana.rodriguez@uab.es (R.R.); montse.nafria@uab.es (M.N.) 2 Institut de Microelectrònica de Barcelona, IMB-CNM, CSIC, 08193 Barcelona, Spain; marcos.maestro@imb-cnm.csic.es * Correspondence: marta.pedro@uab.es  Received: 9 September 2019; Accepted: 22 October 2019; Published: 24 October 2019  Abstract: A fully-unsupervised learning algorithm for reaching self-organization in neuromorphic architectures is provided in this work. We experimentally demonstrate spike-timing dependent plasticity (STDP) in Oxide-based Resistive Random Access Memory (OxRAM) devices, and propose a set of waveforms in order to induce symmetric conductivity changes. An empirical model is used to describe the observed plasticity. A neuromorphic system based on the tested devices is simulated, where the developed learning algorithm is tested, involving STDP as the local learning rule. The design of the system and learning scheme permits to concatenate multiple neuromorphic layers, where autonomous hierarchical computing can be performed. Keywords: memristors; neuromorphic engineering; OxRAM; self-organization maps; synaptic device 1. Introduction The implementation of electronic synapses is nowadays one of the challenges of hardware-based neuromorphic engineering, which aims to design electronic circuits with a similar architecture and behavior to the one found in biological brains. Within this context, the conductivity of an electronic device with memristive characteristics is identified as the weight or strength of a connection between two neurons (Figure 1), usually within a crossbar array which implements the synaptic matrix layer of an electronic neural network (Figure 2a). An analog behavior of the electronic synapse is desirable to improve the robustness of the network [1–3], showing a large window between its higher and lower conductivities and displaying many accessible conductivity levels in between. The conductivity of an analog synaptic device is then finely tuned according to some learning rule during the training stage of a learning algorithm. Among the different technologies that have been proved to be suitable for synaptic applications, the oxide-based resistive random access memory (OxRAM) technology stands out when analog conductivity changes are required [1–7]. Materials 2019, 12, 3482; doi:10.3390/ma12213482 www.mdpi.com/journal/materials Materials 2019, 12, x FOR PEER REVIEW 2 of 18 Materials 2019, 12, 3482 2 of 18 Materials 2019, 12, x FOR PEER REVIEW 2 of 18 Figure 1. Electronic synapses can be implemented with memristive devices. The electronic neural network Figure 1.isElectronic implemented in a synaptic synapses matrix layer. with memristive devices. The electronic neural can be implemented Figure 1. Electronic synapses can be implemented with memristive devices. The electronic neural network is implemented in a synaptic matrix layer. network is implemented 1.1. Spike-Timing in a synaptic Dependent Plasticity matrix (STDP) inlayer. Memristive Electronic Synapses 1.1. Spike-Timing Dependent Plasticity (STDP) in Memristive Electronic Synapses The synapticDependent 1.1. Spike-Timing weight updating Plasticityprocess (STDP)isintherefore Memristive theElectronic basis for Synapses the application of any learning The synaptic algorithm in a neural weight updating network, and isprocess relatedistotherefore the capabilitythe basis of theforsynapse the application to adapt its of any learning conductivity The synaptic algorithmexperience, weight in a neuralnamely updating network, process is therefore the basis for the application of any learning through itsand is related property to the capability of plasticity. In the case of theofsynapse electronic to synapses, adapt its conductivity this feature algorithm in a through experience, neural network, namely and is related its conductivity to property of plasticity. the capability of the synapse to adapt its conductivity involves the modulation of the (G) of an In the casedevice, electronic of electronic wheresynapses, changes (ΔG) this feature can be through theexperience, namely its conductivity property of plasticity. In the casedevice,of electronic synapses, this feature induced by applying the appropriate voltage drop between its two terminals (Figure 2b). can involves modulation of the (G) of an electronic where changes (∆G) Thesebe involves induced by the modulation applying of the the appropriate conductivity voltage (G) of an electronic device, where changes (ΔG) can be updates in the conductivity of the device aredrop between applied its two to according terminals the recent (Figure 2b). of activity These updates the neurons induced in by applying the conductivity thedevice of the appropriate voltage arecorrelations applied drop between according to the recent its two activity terminals (Figure it2b). of theactivity neurons These connects. it connects. For instance, temporal and causality between the recent of the input updates in For instance, the conductivity temporalcan of correlationsthe device are applied according to the recent activity of the neurons and output neurons determineand thecausality magnitude between the recent of and direction activity of the input the relative synapticand weight output it connects. neurons can For instance,the determine temporal magnitude correlations and and causality direction of the betweensynaptic relative the recent activity weight of the∆G/G. change, input change, ΔG/G. The so-called spike-timing dependent plasticity (STDP) has been reported in andso-called The output neurons spike-timing can dependent determine plasticity the magnitude(STDP)and has direction been of the reported relative synaptic in biological weight systems neural [8–16], biological systems [8–16], and is a popular bio-inspired learning rule implemented in artificial change, ΔG/G. and is a popular The so-called bio-inspired learning spike-timing dependent rule implemented plasticity in artificial (STDP) has been reported in networks and computational neuroscience [10–15], where ΔG/G isneural described networks and computational as a function of the time biological neuroscience systems [10–15],[8–16], whereand is ∆G/G a popular bio-inspired learning rule implemented in artificial ∆t between (Figure neural delay Δt between the pre (input) andispost-synaptic described as(output)a function of thespike neurons time firing, delay respectively the pre networks (input) and computational neuroscience [10–15], where ΔG/G is described as a function of the oftime 2c). Theand nature post-synaptic of the synaptic (output) change neurons is whatspike depends firing, on respectively the causality (Figure between2c). theTheinput nature and output the delay synapticΔt between change isthe whatpre (input) depends and on post-synaptic the causality (output) between neurons the input spike and firing, output respectively neurons (Figure activations. neurons activations. 2c). The nature of the synaptic change is what depends on the causality between the input and output neurons activations. Figure 2. (a) Neuromorphic memristive array. Each node within the crossbar corresponds to the weighted Figure 2. (a) connection Neuromorphic (synapse) between memristive two neurons,array. Each node implemented withwithin the crossbar a memristor. (b) Thecorresponds conductivityto of the weighted connection (synapse) between two neurons, implemented with a memristor. device can be changed according to the activity of the neurons it connects. (c) STDP function. (b) The conductivity of the device can be changed according to the activity of the neurons it connects. the Figure 2. (a) Neuromorphic memristive array. Each node within the crossbar corresponds to (c) weighted In STDPorder connection to function. (synapse) demonstrate between two the plasticity neurons, property implementeddevices, of memristive with a memristor. (b) The the input and output activities are assumed to be in the form of voltage pulses, and the significant change in ∆G/G (c) conductivity of the device can be changed according to the activity of the neurons it connects. occurs whenInSTDP these function. orderpulses meet at the terminals to demonstrate of the synaptic the plasticity propertydevice, overlapping of memristive in time, devices, thecausing input aand significant output voltage drop activities are (Figure assumed 2b). to In bethis case, in the the of form STDP function voltage (Figure pulses, and 2c) the can be tunedchange significant by changing in ΔG/Gtheoccurs shape In order of the input to demonstrate the plasticity property of memristive devices, the input and output when these and pulsesoutput meetneuron at thepulses [12–15]. terminals Thesynaptic of the most popular device,shape of the STDP, overlapping resembling in time, causingthea activities one reportedarein assumed to beculture a biological in the byform and of voltage pulses, and shows the significant change inexperimental ΔG/G occurs significant voltage drop (Figure 2b). Bi In thisPoo [16] case, (Figure the STDP2cfunction the average (Figure 2c)ofcan the be tuned by when has data), these pulsesbeen already meet at the in reported terminalselectronic of the synaptic device, overlapping in time,possibility causing a changing the shape of the input andmany output neuron devices [12,13,17–21]. pulses [12–15]. The most However, popularthe shape of the significant to voltage tuneresembling the STDP drop (Figure concerning function 2b). In thisthecase, the STDP function (Figure 2c) can be tuned by STDP, the one shape, reported in a biological electronic culture synapse by Bi and electrical Poo [16]characteristics, (Figure 2c shows is often the changing the shape of the input and output neuron pulses [12–15]. The most popular shape of the STDP, resembling the one reported in a biological culture by Bi and Poo [16] (Figure 2c shows the Materials 2019, 12, 3482 3 of 18 skipped. Variety in STDP functions appears in biological synapses [8,12–15]. This variety extends the application of the STDP as a local learning rule in artificial intelligence learning schemes [12], especially in those based on unsupervised techniques. 1.2. Unsupervised Learning and Self-Organizing Neural Networks Unsupervised learning involves a methodology where the training stage does not require the calculation of any error made by the system for a certain input, in order to improve its performance. That is, both user and system are not meant to know the actual solution of the problem entered to the network, nor detailed information about the input dataset properties, in contraposition with supervised learning techniques. Unsupervised learning implementation would be beneficial for neuromorphic architectures, since on the one hand, it does not rely on the error computation and correction as the supervised learning techniques do, so extra circuitry specialized for this purpose could be avoided. Furthermore, unsupervised learning models, such as the above mentioned STDP learning rule, are considered to be biologically plausible. By reverse engineering simple and primitive biological nervous systems as a first approach, the neuromorphic community would take advantage and inspiration because of the simplicity of their design, compared to the ones found in artificial deep learning neural networks, which present an extremely high density of synapses, neural layers and complex pathways and dynamics. Applications of unsupervised learning algorithms are related to classification, symbolic representation, and associative tasks, usually by extracting the relevant statistical features of the input dataset. Examples of bio-inspired unsupervised learning implementations based on memristive devices for image recognition tasks can be found in [17–21]. However, the hardware-based implementation of other unsupervised learning applications remains unexplored. A particular example of bio-inspired unsupervised learning is the self-organizing map (SOM), also called Kohonen network [22]. Applications of SOM extend to financial predictions, medical diagnosis, or data mining, among others [22–24]. The aim of this learning algorithm consists in mapping the input dataset onto a regular and usually two-dimensional grid, which corresponds to the output layer, under an unsupervised and competitive learning scheme. A diagram of a Kohonen network is depicted in Figure 3a. In here, the input layer is unidimensional and consists of three nodes (input neurons). The output layer is bidimensional, and each node corresponds to an output neuron. Output neurons can communicate to their immediate neighbors. All of the input nodes have a weighted connection (synapse) with every output node. The weight of the synapse determines how strong an output neuron responds to the activation of a particular input. These neural networks are inspired in the topological maps found in the sensory-processing areas of the brain (Figure 3b), where neurons that respond to similar inputs are spatially located very close. The key of this algorithm consists in evaluating the similarity between the set of weights of an output neuron and the input data, which is fed to the system as a vector. The original software algorithm consists in the sequential execution of the following steps, parting from a network with randomly initialized weights. For randomly chosen input from a particular dataset, the Euclidean distance between the input and the weights of every output neuron must be computed, in order to determine which is the output neuron whose weights are closer to the input. This element is identified as the best matching unit (BMU) and its weights are updated in order to slightly reduce its distance with the input data. Once trained, these networks present topographical organization such as the one found in sensory processing areas of the brain, such as the tonotopic map found in the primary auditory cortex, in charge of processing sound (Figure 3b). In here, the neurons that respond to similar sound frequencies are grouped in clusters, which appear in a frequency-ordered fashion. In this way, similar inputs activate neurons in the output layer which are found close to each other, whereas dissimilar ones affect distant regions [22,25–27]. The output layer neurons in a trained SOM appear organized in clusters, whose relative size and location provides statistical information of its corresponding input data item characteristics. It is actually the presence or absence of an active response of an output Materials 2019, 12, 3482 4 of 18 Materials 2019, 12, x FOR PEER REVIEW 4 of 18 neuron cluster, and not so much the exact input–output signal transformation or magnitude of the neuron cluster, and not so much the exact input–output signal transformation or magnitude of the response, that provides an interpretation of the input information [22–24]. response, that provides an interpretation of the input information [22–24]. Many methods are derived from the SOM algorithm, where the neural system is built with SOMs Many methods are derived from the SOM algorithm, where the neural system is built with SOMs as basic blocks or layers, such as the multi-layer or hierarchical SOM (HSOM) [22]. In the latter as basic blocks or layers, such as the multi-layer or hierarchical SOM (HSOM) [22]. In the latter case, case, the network is constituted by concatenating SOMS in a feed-forward way (cascade), where one the network is constituted by concatenating SOMS in a feed-forward way (cascade), where one SOM SOM layer is trained by receiving as input the outputs of another previous SOM. The advantage of layer is trained by receiving as input the outputs of another previous SOM. The advantage of HSOMs HSOMs is that they require less computational effort than a standard SOM to perform certain tasks or is that they require less computational effort than a standard SOM to perform certain tasks or problems that present a hierarchical or thematic structure, and moreover, HSOMs provide a simpler problems that present a hierarchical or thematic structure, and moreover, HSOMs provide a simpler representation of the results, which leads to an easier interpretation because they allow the user to representation of the results, which leads to an easier interpretation because they allow the user to check what clustering has been performed at each level of the hierarchy. check what clustering has been performed at each level of the hierarchy. Figure 3. Figure 3. (a) (a) Example Example of of aa self-organizing self-organizing map. map. (b) An example (b) An example of of aa topological topological map map in in the the human human brain, corresponding corresponding to to the the tonotopic tonotopic map map of of the the primary primary auditory auditorycortex. cortex. In this work, work, we we propose proposean anunsupervised unsupervisedhardware hardwareadaptation adaptationofof thethe SOMSOM algorithm algorithm to to be be implemented implemented in on-line in an an on-line learning learning neuromorphic neuromorphic OxRAM-based OxRAM-based crossbar crossbar array,array, by means by means of of bio- bio-inspired unsupervised inspired unsupervised learning learning methods, methods, being being thethe first first ofofitsitskind, kind,totothe thebest bestof of our our knowledge. There There isis another another work work related related toto the the electronic electronic implementation implementation of the SOM algorithm: [28] is also a simulation work, and is based upon the previous calculation of the desired synaptic weight update, hence not being being an an unsupervised unsupervised learning learning algorithm. algorithm. InIn contrast, contrast, in in our work we provide a fully unsupervised learning algorithm, in which the weight updating process relies on the STDP property of the the employed employed memristive devices. For the sake of simplicity, a very simple input dataset is used used as an example, for which a color identification task is is provided. provided. First of all, a model from a previous work work [29] [29] is is used used to to analyze analyze thethe plasticity plasticity property property ofof the the tested tested devices, devices, which which isis further further verified verified experimentally. experimentally. A methodology for tuning the STDP function, which is a key element to control the learning process, is proposed. proposed. The obtained STDP curves are used as the local learning learning rules within within the adapted adapted SOMSOM algorithm, algorithm, forfor which which aa fundamental fundamental application application isis demonstrated. demonstrated. The learning learning mechanisms mechanismsintroduced introducedin in thisthis work can concatenate work multiple can concatenate SOMs without multiple extra circuitry, SOMs without extra providing circuitry, aproviding step towards a stepthetowards implementation of hardware-based the implementation hierarchicalhierarchical of hardware-based computing computing systems. systems. 2. Materials and 2. Materials and Methods Methods 2.1. Electrical Characterization and Device Modeling 2.1. Electrical Characterization and Device Modeling The devices employed in this study are TiN/Ti-HfO2 -W metal–insulator–metal (MIM) structures. The devices employed in this study are TiN/Ti-HfO2-W metal–insulator–metal (MIM) structures. They were fabricated on silicon wafers either with an oxide isolation scheme or as a single crossbar They were fabricated on silicon wafers either with an oxide isolation scheme or as a single crossbar on a thermally grown 200 nm-thick silicon dioxide. The 10 nm-thick HfO2 layer was deposited by on a thermally grown 200 nm-thick silicon dioxide. The 10 nm-thick HfO2 layer was deposited by atomic layer deposition at 225 ◦ C using TDMAH and H2 O as precursors, and N2 as carrier and purge atomic layer deposition at 225 °C using TDMAH and H2O as precursors, and N2 as carrier and purge gas. The top and bottom metal electrodes were deposited by magnetron sputtering and patterned by gas. The top and bottom metal electrodes were deposited by magnetron sputtering and patterned by photolithography. The bottom electrode (BE) consists of a W layer and the top electrode (TE) of TiN on photolithography. The bottom electrode (BE) consists of a W layer and the top electrode (TE) of TiN on a 10 nm-Ti layer acting as oxygen getter material. The fabricated devices are square cells with an area of 5 × 5 μm2. Figure 4a shows a scanning electron microscope (SEM) (IMB-CNM (CSIC), Materials 2019, 12, 3482 5 of 18 a 10 nm-Ti layer acting as oxygen getter material. The fabricated devices are square cells with an area Materials 2019, 12, x FOR PEER REVIEW 5 of 18 of 5 × 5 µm2 . Figure 4a shows a scanning electron microscope (SEM) (IMB-CNM (CSIC), Barcelona, Spain) image Spain) Barcelona, of the tested imagestructures, where of the tested the TE and structures, BE the where are TE indicated. and BE More details onMore are indicated. the electrical details on behavior and fabrication process of these samples can be found in [30,31]. the electrical behavior and fabrication process of these samples can be found in [30, 31]. 5 µm Figure 4. (a) SEM image of the tested structure (b) Experimental I-V characteristics of the analyzed Figure[29]. devices 4. (a)ASEM image voltage of for limit thethe tested structure RESET (b)was process Experimental I-Vwhereas set to −1.6V, characteristics for theof theprocess, SET analyzed thedevices [29]. A voltage limitparameter conductivity-controlling for the RESET wasprocess was set tocurrent the maximum −1.6V, whereas driving for thethe SET process, device (currentthe conductivity-controlling parameter compliance, Ic ) set by the user. was the maximum current driving the device (current compliance, Ic) set by the user. In Figure 4b, a few examples of experimental I–V curves are shown, where it can be noted that In Figure the tested devices4b, a few examples display of experimental a bipolar resistive switching I–Vbehavior, curves are shown, where consisting it can be in transitions noted from highthat the tested (HRS) to lowdevices display a bipolar (LRS) resistance states andresistive switching vice versa. These behavior, consisting transitions in transitions are identified as the SET fromandhigh (HRS)processes, RESET to low (LRS) resistance respectively. Thestates mainand viceofversa. results Thesework a previous transitions [30] showarethat identified as the SET small changes and in the RESET processes, conductivity at the lowrespectively. resistanceThestatemain (LRS) results can be ofinduced a previous work [30] by means show that small of controlling changes in the maximum the conductivity current at the during driving the devices low resistance state (LRS) the SET process, provingcan be plasticity their induced property, by means andofthus controlling indicatingthe maximum that the testedcurrent devices driving the devices are suitable to playduring the SETrole the synaptic process, proving theircrossbar-array. in a neuromorphic plasticity property, and In [29], thus indicating that setup a pulse-programming the tested devices are was proposed, suitable with the aim to of play the synaptic analyzing in whichroleways in a fine neuromorphic changes in crossbar-array. the conductivity In of [29], thea device pulse-programming can be induced setup wasapplication by the proposed, with the aim of single of analyzing pulses. in which The proposed setup allowed obtaining the experimental G–V characteristics of the tested devices, by means ofpulses. ways fine changes in the conductivity of the device can be induced by the application of single the The proposed application setup allowed of increasing obtaining amplitude and decreasing the experimental G–V characteristics single pulses of the testedover with a fixed pulse-width devices, time.by meansfrom Results of the application [29] are shownof inincreasing and decreasing Figure 5, where amplitude and the pulse amplitude single thepulses with a measured conductivity fixed pulse- width after everyover time. single Results applied from pulse (in[29] are shown Go units, beingin = 77.5 5, GoFigure µSwhere the pulse the quantum amplitude and of conductance unit)the areconductivity plotted againstmeasured after every the number singlepulses. of applied appliedThe pulse (in Go units, conductivity being state Go = measured G was 77.5 μS theafterquantum the of conductance application of everyunit) pulseare(Figure plotted5a,against the number red pulses), of applied by means pulses. of applying The (Figure 50mV conductivity 5b, gray state G was pulses) andmeasured reading theaftercurrent the application of every flowing through thepulse (Figure device. 5a,analyzed In the red pulses), voltageby means of applying 50mV range, conductivity can (Figure take values 5b, gray pulses) between ~10 Go and andreading 30 Go.the current flowing through the device. In the analyzed voltage range, conductivity By means can take values of representing between the obtained ~10 Go and conductivity experimental 30 Go. as a function of the applied voltage, the experimental the G–V characteristics can be fitted according to the compact model of [32]. In here, the so-called hysteron function is used to describe a time-independent conductivity window as a function of the applied voltage in non-linear memristive devices. An example of an ideal hysteron function of a non-linear memristive device is depicted in Figure 6a. The normalized internal state λ is represented as a function of voltage drop at the memristor. The top and bottom boundaries are identified as the maximum (gmax ) and minimum (gmin ) conductivity states. In order to increase (decrease) the conductivity state of the device, a positive (negative) voltage has to be applied so that λ shifts towards gmax (gmin ), describing the Γ+ (Γ− ) trajectories. The pair of logistic ridge functions Γ+ and Γ− can be modeled with two cumulative distribution functions (cdf) [29], related to the pulse amplitudes applied to the non-linear memristive device, being V+ , σ+ and V− , σ− the average and standard deviation values of the cdf related to Γ+ (for dV/dt > 0) and Γ− (dV/dt < 0) curves, respectively. Figure 5. (a) Stair-case pulse-programming scheme used in [29] for obtaining the G–V characteristics of a memristive device. (b) The pulse amplitude was increased and decreased over time to change the device conductivity. (c) Conductivity G as a function of the applied pulses. ways fine changes in the conductivity of the device can be induced by the application of single pulses. The proposed setup allowed obtaining the experimental G–V characteristics of the tested devices, by means of the application of increasing and decreasing amplitude single pulses with a fixed pulse- width over time. Results from [29] are shown in Figure 5, where the pulse amplitude and the conductivity measured Materials 2019, 12, 3482 after every single applied pulse (in Go units, being Go = 77.5 μS the quantum 6 of 18 of conductance unit) are plotted against the number of applied pulses. The conductivity state G was measured after the application of every pulse (Figure 5a, red pulses), by means of applying 50mV Both of them define the boundaries of the possible conductivities of the device within a range limited (Figure2019, Materials 5b, gray pulses) 12, x FOR PEERand reading the current flowing through the device. In the analyzed voltage REVIEW 6 of 18 by the minimum and maximum conductivity states, gmin and gmax , respectively. range, conductivity can take values between ~10 Go and 30 Go. By means of representing the obtained experimental conductivity as a function of the applied voltage, the experimental the G–V characteristics can be fitted according to the compact model of [32]. In here, the so-called hysteron function is used to describe a time-independent conductivity window as a function of the applied voltage in non-linear memristive devices. An example of an ideal hysteron function of a non-linear memristive device is depicted in Figure 6a. The normalized internal state λ is represented as a function of voltage drop at the memristor. The top and bottom boundaries are identified as the maximum (gmax) and minimum (gmin) conductivity states. In order to increase (decrease) the conductivity state of the device, a positive (negative) voltage has to be applied so that λ shifts towards gmax (gmin), describing the Γ+ (Γ-) trajectories. The pair of logistic ridge functions Γ+ and Γ- can be modeled with two cumulative distribution functions (cdf) [29], related to the pulse amplitudes applied to the non-linear memristive device, being V+, σ+ and V−, σ− the average and Figure standard 5. (a) Stair-case deviation valuespulse-programming scheme of the cdf related to Γ+ (forused in [29] dV/dt > 0)for andobtaining the