A Model of an Oscillatory Neural Network with Multilevel Neurons for Pattern Recognition and Computing
Electronics , Article
A Model of an Oscillatory Neural Network with Multilevel Neurons for Pattern Recognition and Computing
Andrei Velichko *, Maksim Belyaev and Petr Boriskov
Institute of Physics and Technology, Petrozavodsk State University, 31 Lenina str., 185910 Petrozavodsk, Russia; [email protected] (M.B.); [email protected] (P.B.) * Correspondence: [email protected]; Tel.: +7-8142-63-5773 Received: 20 November 2018; Accepted: 05 January 2019; Published: date
Abstract:
The current study uses a novel method of multilevel neurons and high order synchronization effects described by a family of special metrics, for pattern recognition in an oscillatory neural network (ONN). The output oscillator (neuron) of the network has multilevel variations in its synchronization value with the reference oscillator, and allows classification of an input pattern into a set of classes. The ONN model is implemented on thermally-coupled vanadium dioxide oscillators. The ONN is trained by the simulated annealing algorithm for selection of the network parameters. The results demonstrate that ONN is capable of classifying 512 visual patterns (as a cell array 3 × 3, distributed by symmetry into 102 classes) into a set of classes with a maximum number of elements up to fourteen. The classification capability of the network depends on the interior noise level and synchronization effectiveness parameter. The model allows for designing multilevel output cascades of neural networks with high net data throughput. The presented method can be applied in ONNs with various coupling mechanisms and oscillator topology.
Keywords: oscillatory neural networks; pattern recognition; higher order synchronization; thermal coupling; vanadium dioxide; computing
1. Introduction
Hypotheses about the functional importance of synchronization for information processed by the brain were put forward long ago [1,2] and its experimental discovery [3,4] encouraged the creation of neural networks models with oscillatory dynamics [5–7] and neuromorphic algorithms of image processing based on synchronization [8–11]. Research on ONN is mainly based on phase oscillator models. Such models, primarily the Kuramoto model [12], have been very useful for studying systems of various topology with a large number ( N > 10 ) of oscillators, where not only global synchronization but the synchronization by middle field [13], mode of quasi-synchronization [14] and even chimeras synchronization [15] are feasible. Recently, the problem of pattern recognition in ONN has been intensively studied [16,17], and two main global synchronization methods have been outlined: frequency-shift [10] and phase-shift [11] keying methods. Frequency encoding, based on synchronized frequency shift [10], allows usage of an oscillator star configuration and N couplings; however, the disadvantage of this method is that only it stores a single pattern. The phase-shift keying method of encoding [11] enables storage of more than one pattern by certain combinations of phase shift at the same weight ONN matrix. However, this method has the following drawbacks: N couplings with tunable weights and a two-stage procedure of pattern recognition. lectronics , , x FOR PEER REVIEW 2 of 26 Another class of ONN is based on relaxation oscillators that generate multiple pulses of short duration and fixed amplitude. Such pulses (spikes) can code information at pulse-repetition frequency. An important distinguishing feature of a pulse type ONN from classic spiking neural networks (SNN) is a self-oscillating mode of ONN neuro-oscillators. This mode is not indicative for SNNs and real (biological) neurons, which are characterized by a forced response through generation of a single spike or their group, when the neuron threshold level is achieved. However, ONNs are fascinating due to the simplicity of the hardware implementation as available advanced micro- and nano-electronic self-oscillators ensure the compact size and energy efficiency of the circuit. Pulse type ONNs with a multi-frequent spectrum of periodic oscillation feature a specific mode of multiple frequency synchronization, or in other words, a high order synchronization effect [18] (harmonic locking [6]). We demonstrated this effect experimentally by using thermally-coupled vanadium dioxide (VO ) oscillators [19]. In relaxation oscillators with VO -based film elements, electric self-oscillations are activated by the effect of electric switching governed by metal—insulator transition (MIT) [20–23]. The processing speed of VO devices switching which amounts to ~10 ns [24] and manufacturing technology that allows switching elements to be created with high levels of nano-scalability make VO -switch based oscillators the perfect objects for research on neuro-oscillators to solve cognitive technology problems [25–28]. Relaxation oscillators with high order synchronization effects can be realized by using electric coupling [25–28] and by using not only VO -switches, but any other switching elements such as thyristors [29], tunneling diodes [30], resistive memory cells [31], and spin torque oscillator [16,32]. In this paper, we studied the ONN of thermally-coupled VO -oscillators and present a general concept of visual pattern recognition based on high order synchronization effects. We used the multiplicity of synchronous states to extract object classes by using a single output oscillator (compared to, for example, an array of oscillators at the output [10,11]). The concept of a multi-level neuron allows for using a smaller number of output neurons (oscillators) to implement the more complex cognitive functions of the neural network. We developed a set of special metrics [19,33], such as the high-order synchronization value and the synchronization effectiveness value. Compared to the neural network presented in [33], this network is able not only to memorize and classify patterns, but also to perform logical operations, computer calculations and emulate other functions of artificial intelligence systems.
2. Materials and Methods
The basic element in the studied ONN is an oscillator implemented on the circuit of a relaxation generator (Figure 1) based on a VO -switch. We described the process of fabrication and electric I–V characteristics of an electric VO -switch in detail in [19]. I-V characteristics are well approximated by a piecewise linear function (see Appendix A.1) that has two conduction states (low-resistance and high-resistance) and a region of negative differential resistance (NDR). These switches may be used in the circuit of a relaxation generator with power supply I p holding the operation point in the NDR region of the I-V characteristic and with capacity C , connected to the switch in parallel. In addition, a source of noise U n models the circuit’s interior or exterior noises, such as current noise of a switch [34]. The oscillator’s output signal is current I sw ( t ) flowing through the VO -switch, which is used to calculate the synchronization level of two oscillators (see Section 2.3). The current signal directly determines the effect of thermal coupling inside the network. lectronics , , x FOR PEER REVIEW 3 of 26 Figure 1.
The oscillator circuit and example of oscillators’ interaction via thermal coupling of VO switches. I sw ( t )—current signal in a VO switch, that leads to its Joule heating resulting in thermal impulses that spread along the substrate; U n —source of noise, s —thermal coupling strength. We applied thermal coupling to connect oscillators in a network. The presence of the thermal coupling mechanism between two VO oscillators was convincingly demonstrated by us in the experiment [19], and it is based on the mutual thermal effect of switches due to their Joule heating and the dependence of the threshold voltage U th on temperature. In the pre-threshold mode, when the VO temperature of the switch channel is close to the MIT temperature T th ~ 340 K [23], the external thermal influence can “push” it to switch, which is equivalent to lowering the threshold voltage by the value of s , the thermal interaction force. The switching causes an even higher temperature rise in the switch, which leads to the appearance of a temperature pulse in the surrounding space, which propagates as a temperature wave. In the oscillator circuit, the thermal effect of the switches on each other occur in the mode of repetitive pulses initiated by self-oscillations of oscillator currents, which ultimately leads to their synchronization [19]. The value of s can be controlled experimentally, for example, by varying the distances between the switches or by varying the parameters of the external circuit [19]. The choice of this coupling type is determined by the simplicity of a computing model when oscillators are electrically separated from each other, unlike in capacitive or resistive couplings [21,22,25,28]. An example of an oscillators’ interaction via thermal coupling of VO switches is shown in Figure 1. A detailed presentation of the mathematical model of thermally coupled relaxation oscillators is given in Appendix A.1. The heat assisted model we used completely describes the experimental data previously observed [19,35,36]. Indeed, the speed of propagation of a thermal wave can be taken into account, but this is not in the scope of the current paper because we operated at low oscillation frequencies where the time delay is insignificant. We considered the radius of thermal interaction R TC in [19], therefore, in this model we limited ourselves to interaction with the nearest oscillators. Consideration of more complex cases can be done in future publications. We used the concept of one-way thermal coupling of oscillators in the numerical model. This can be physically implemented when a resistance heater is used as a connecting element in one of the circuits instead of a switch [19]. For numerical modeling we used the following parameters as I-V characteristics ( U th = 5 V, U h = 1.5 V, U bv = 0.82 V, R off = 9.1 kΩ, R on = 615 Ω). In this circuit, capacity is a constant parameter C = 100 nF, and its value significantly determines the frequency range of oscillator operation [37] and its natural frequency F . In our case, the frequency range was 165 Hz ≤ F ≤ 1266 Hz at the range of feeding currents 550 µA ≤ I p ≤ 1061 µA. lectronics , , x FOR PEER REVIEW 4 of 26 The studied ONN consists of an input pattern, which is transferred as a 3 × 3 matrix on the layer of processing neurons consisting of 9 oscillators, and of an output layer with only one oscillator (output neuron) No.10 (see Figure 2).
Figure 2.
ONN organization circuit for pattern recognition as a matrix of 3 × 3 elements. Digits indicate the sequence numbers of oscillators.
The magnitude of the effect of the i -th oscillator on the j -th oscillator is determined by coupling strength s i,j , that is set by the following matrix: r r r r r r r r r rm m om m om m o0,0 0,10 m m m om m m m o10,0 10,10 m m m om m om m m om m o
00 0 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 m s s s s s s s s s ss s ss s s ss s ss s s s s sS s s s s ss s s s s ss s ss s s ss s s (1) The layer of processing neurons consists of a 3 × 3 oscillator matrix connected by similar couplings ( s i,j = s j,i = s m , where i , j are the numbers of neighboring oscillators). Neighboring oscillators are only connected by horizontal and vertical lines. So, the central oscillator No.5 has four couplings in the matrix, the corner oscillators have two couplings, and the oscillators in the center of the edges have three couplings (with the central oscillator and two corner ones) and they all unidirectionally affect oscillator No.10 (output neuron) in the output layer ( s i,10 = s o and s = 0, where i = 1 … 9). Importantly, there is the reference oscillator No.0 in the circuit and the synchronization order of all other oscillators is measured against this oscillator. Oscillator No.0 (Figure 2) unidirectionally affects all other oscillators ( s = s r , and s i,0 = 0, where i = 1 … 10). The input pattern (see Figure 3) is transferred to the processing layer by selection of the feeding currents of the oscillator matrix: lectronics , , x FOR PEER REVIEW 5 of 26 ON ip_i OFF i , if 1, if 0
I xI I x (2) where i = 1 … 9, x i are the coordinates of the input vector X = ( x , …, x ), that correspond to white ( x i = 0) and blue ( x i = 1) colors of pattern squares. Transfer of the pattern to the intermediate layer causes the change in the oscillators’ feeding currents, and in turn, it leads to the change in synchronization state for all oscillators (No.1–10). The synchronization order SHR between the reference oscillator No.0 and the oscillator in the output layer No.10 serves as the control value. The values of the feeding currents for these oscillators are fixed and may differ from currents in the matrix, therefore they are indicated as I and I (see Figure 3). Figure 3.
Principle of a pattern transfer from the input layer to the oscillator matrix via setting of their currents ( I OFF is for white squares, I ON is for blue squares). Separate colors are used for oscillators No.0 and No.10 and for their currents I and I , respectively. In this section, we present the method of calculating a family of metrics [19] to determine the high order synchronization of two oscillators. A family of metrics has two basic parameters: the ratio of subharmonics (SHR) and synchronization effectiveness η . In a general case, it is possible to use the concept of subharmonics ratio between oscillators i and j , which is defined [18] as a simple fraction i,j i j SHR : k k (3) where k i and k j are harmonics order of oscillators at the common frequency of their synchronization F s . In other words, Equation (3) uses spectral approaches to describe synchronization of the higher order k i : k j , when the following relation is observed: F k F k F (4) where F i0 , F j0 are frequencies of main harmonics. As an example, Figure 4 shows the qualitative spectra of two oscillators at SHR i,j = 2:7 lectronics , , x FOR PEER REVIEW 6 of 26 Figure 4.
Qualitative spectra of two coupled synchronized oscillators with parameter SHR i,j = 2:7.
In addition, Figure 4 shows the total synchronization frequency F s and subharmonics numbers at this frequency k i = 2 and k j = 7. Usage of signal spectra for estimating the magnitude of SHR i,j is not effective because the signals might not have a strictly periodic sequence and when noise is added to the system, the spectral lines broaden. Below, we will describe the mathematical procedure that estimates the value of SHR i,j without the use of spectral characteristics but using a current signal and array LE. Array LE stores information on the position of the current pulse leading edges (see Figure 5). Figure 5.
Oscillogram of oscillator current No. 10 and the corresponding array of positions of the leading edges of the current pulse LE [10] [Δt]. Δt—calculation time interval (Appendix A.1).
Figure 6 shows arrays LE[ i ] and LE[ j ] for two oscillators that correspond to the same ratio of the basic frequencies as shown in the spectra in Figure 4. The distance between two nearest phase-locked pulses is denoted as T z s —the period of synchronization (where z is a conditional number of periods T s ). oscillator jk i =2 A m p li t ude F F s oscillator ik j =7 A m p li t ude Frequency F SHR i,j = k i : k j =2:7 L E [ ][ t ] time, t -3 -3 -3 -3 I s w _10 , A time, s lectronics , , x FOR PEER REVIEW 7 of 26 Figure 6.
Arrays LE[ i ] and LE[ j ] for two oscillators. Therefore, SHR i,j may be estimated using a phase-locking method: i,j j i
SHR :
M M (5) where M i and M j are the number of signal periods falling into the synchronization periods T z s of two oscillators. Formula (5) can be easily obtained from Formula (3) taking into account Formula (4) and the following ratio ( T s = M i ·T i0 = M j ·T j0 , F i0 = 1/ T i0 , F j0 = /T j0 ). In general, especially when a system behaves erratically, synchronization periods differ and spread in T z s ≠ T z +1s and the values of M i and M j may change within one oscillogram (see Figure 7). Figure 7.
Arrays LE[ i ] and LE[ j ] for two oscillators with a non-constant period of synchronization T z s . Various values of synchronizations SHR i,j may occur within one oscillogram. To determine which SHR i,j value prevails, it is necessary to find the occurrence probabilities P ( M j : M i ) for each pair ( M i : M j ) that is present in the whole oscillogram and to select the pair with the maximum value of P = P max ( M j : M i ). Then, the final value of SHR i,j will be written as: i,j j i max j i SHR : , if ( : )
M M P P M M (6) L E [ i ] SHR i,j =2:7 M j =2 L E [ j ] time oscillator i oscillator j T zs M i =7 L E [ i ] T z+3s M j =2 M i =7 M j =2 M i =5 T z+2s M j =2 M i =9 T z+1s M j =2 L E [ j ] time oscillator i oscillator j T zs M i =7 lectronics , , x FOR PEER REVIEW 8 of 26 To find the probabilities P ( M j : M i ), we can count how many times NP ( M j : M i ) the given pair appeared within the whole oscillogram of the oscillator i , multiply by the number of periods in it ( M i ) and divide by the total number of all oscillations periods in the given signal ( N j ). Thus, for P ( M j : M i ) we obtain: j i j i i i ( : ) 100% ( : ) / P M M NP M M M N (7) where N i is the total number of periods in the oscillogram of oscillator i . It is convenient to present the probabilities P ( M j : M i ) as a histogram where the values are positioned in the descending order of the magnitude P . For example, for the oscillogram section in Figure 7 the following histogram can be built. The histogram in Figure 8 is calculated by Formula (7), when the pairs occur the following number of times, NP (2:7) = 2, NP (2:9) = 1, NP (2:5) = 1, and the total number of periods is N i = 28 (in real calculations, N i was in the range of 1000–3000 for greater accuracy, see Appendix A.2). Figure 8.
Histogram of probabilities distribution P ( M j , M i ) calculated by using Formula (11) for signals LE, shown in Figure 7. For SHR i,j , the parameter of synchronization effectiveness η is defined as the maximum probability P max ( M j : M i ): max j i ( : ) P M M (8) Therefore, we define the synchronization calculated above as SHR i,j = 2:7 with effectiveness of η = 50%. The family of metrics (SHR i,j , η ) allows sufficient determination of the synchronization state of two oscillators and calculation of the distance between the states, i.e., the difference between the synchronization degree. This property allows the use of metrics in an oscillator neural network training, pattern recognition systems and artificial intelligence [8–11]. Depending on the task, for example, the network training for data coding and pattern recognition, the problem of the presence or absence of synchronization SHR i,j can be solved by formally setting the synchronization effectiveness threshold η th , so th th synchronized, if signals are not synchronized, if (9) (2:7) (2:9) (2:5)020406080100 P ( M j : M i ) , % ( M j : M i ) =50% lectronics , , x FOR PEER REVIEW 9 of 26 In the majority of cases, we set η th = 90%, meaning the signals are synchronized if 90% of their durability have a certain synchronization pattern. For the network training, this parameter can be chosen within a selected range, and it is one of the important parameters of the network adjustment [33]. The main technical problem we faced, was the problem of defining the synchronization order between the reference oscillator No.0 and the oscillator of the output layer No.10 characterized by the value SHR : SHR : ~ / k k k k (10) The same value of SHR can be expressed in several ways: as a ratio, a simple fraction or a real number, for example, SHR = 10:3 = 10/3 = 3.33. Later, we will use this property to present the results more vividly. Parameter SHR has the properties of an output neuron while the reference neuron may be considered as a master generator to which we calculate synchronization of other network neurons. Two parameters SHR and η are used as the main metrics for evaluating the degree of two oscillators’ synchronization and are applied in the algorithm of ONN training. The current oscillograms I sw ( t ) of oscillators No.0–10 were calculated simultaneously and contained ~250,000 points with the time interval Δ t = 10µs (see Appendix A.1). Then, the oscillograms were automatically processed. The switch parameters were unchanged in numerical simulation (see Section 2.1), while current intensities I p_i ( I ON , I OFF , I I ), coupling strength constants s ( s r , s m , s o ), noise amplitude U n and η th varied. A “black-and white” pattern was used as an input test pattern presented as a 3 × 3 matrix (without gradation of gray color, 3 by 3 pixels). The form of the pattern can be unequivocally defined by the input vector X n = ( x , …, x ) where each cell may take the value x i = 0 (white color) or x i = 1 (blue color), and n is the number of the vector equal to the decimal value of the coordinates presented as a binary sequence. For example, Figures 2 and 3 show an input pattern that corresponds to the vector X = (1,1,1,1,0,1,0,0,1). The total number of patterns (vectors) n in the input layer matrix is 2 = 512 (X … X ). Presuming that the pattern processing layer together with the output layer has certain symmetry, a set of 512 vectors X n may be divided into 102 classes C j , where j is the number of classes ( j =1..102) (see Figure 9). The complete list of classes and their elements is described in Supplementary Materials (Data1.txt). Figure 9.
Symmetry-based distribution of 512 figures (vectors X n ) into 102 classes C j . The principle of figure distribution into classes is as follows: each class C j consists of a lot of figures (vectors) that have the same number of blue (white) cells and rotation-reflection axial symmetry of the 4th order (symmetry at rotation by 90°). lectronics , , x FOR PEER REVIEW 10 of 26 Mirror-rotation axial symmetry of the 4th order supposes an association of all figures in a separate class in the mirror operation in relation to the central columns (horizontally and vertically) and also at rotation by 90° (see the example of figures transformation for class C in Figure 10). Figure 10.
Mirror-rotation axial symmetry of the 4th order in class C . Figures from one class impose the same effect on the neural network and have the same output value SHR . Distribution into classes allows us to find figures that at any current values ( I ON , I OFF , I , I ) and coupling strength constants ( s r , s o , s m ), noise amplitude U n and η th , will have the same values of synchronization effectiveness η and SHR within one class of figures as a result of neural network symmetry. The initial distribution of all 512 figures into classes allows us to recognize not one specific figure, but a class (out of 102 possible ones) into which this figure is classified. For example, for class C the equality SHR (X ) = SHR (X ) = SHR (X ) = SHR (X ) will be observed. Such classification reduces the number of possible variants to sort as we may send only 102 figures to the input layer, one figure per each class. The problems this neural network is able to solve may be divided into several variants: I. Synchronization of oscillators No.0 and No.10 with the corresponding value of SHR and η > η th , exists only for one specific class C j with number j = m out of 102 classes: : and if =no syncronysation an S dHR (C )= < if k k j m j m (11) Here we have to show the solutions of this problem with various values of m . II. There is a set of classes C = {C Z1 , C Z2 … C ZP }, where Z , Z … Z P are arbitrary non-repeating indices, where the number is P < 102. When inputting this set into the oscillator system, it comes to the synchronization states corresponding to the set SHR = {SHR (1)0,10 , SHR (2)0,10 … SHR (P)0,10 }. The set
SHR does not have the same elements, i.e., each class of figures from set C corresponds to a unique synchronization order SHR . By analogy with (4) the problem may be expressed as: (1)0,10 th j 1(2)0,10 th j 20,10 j (P)0,10 th j SHR and if C and = ,SHR and if C and = ,SHR (C )= .............................................................................SHR and if C a j Zj Z
C C C
Pth j nd = ,no syncronysation and < if C j Z C (12) In fact, problem I is a subcase of problem II at P = 1. The principle of solving problem II is shown in Figure 11. Here P = 3, and set C = {C , C , C }, in this case the corresponding set is SHR = {16:15, 13:10, 14:13}. Therefore, unambiguous recognition of mirror symmetryrotational symmetry lectronics , , x FOR PEER REVIEW 11 of 26 figures belonging to three different classes takes place. Synchronization is not realized for all other classes and η < η th . Figure 11.
The principle of class recognition of neural network figures with one output neuron.
As there is a huge variety of set C variants, this problem can be reduced to the search of possible solutions with P > 1 and to the determination of the maximum value of P . III. The third variant of the problem corresponds to a fully trained network when it solves problem II for all possible input classes C j , when P = 102. In this problem, each class out of 102 possible variants corresponds to its unique value of synchronization order SHR . Therefore, the set of input classes C = {C , C … C } transfers into a set of synchronous states SHR = {SHR (1)0,10 , SHR (2)0,10 … SHR (102)0,10 }, where all the elements have non-duplicate values. Problems I-III have an increasing degree of complexity and are subcases of problem II with different parameter P . Problem I is the simplest one and it may be solved by a common neural network with one bistable output neuron (bistability means the presence or absence of synchronization), although two variants of answers in neural networks are often given by using two output neurons [38]. The output neuron in problems II and III have multilevel properties and differ from a common bistable neuron. All three problems can be applied to the circuit shown in Figure 2. The ONN circuit has only one output neuron, nevertheless, the multilevel high order synchronization SHR allows input pattern classification into P classes within set C . This is the most striking difference between the described neural network and variants presented in the literature. This increases the net data throughput of a single neuron and enables us to create multilevel output cascades of neural networks with high functionality. To solve problems I–III, network training is required. As ONN with a high order synchronization effect has not been studied before, there are no established methods for network training. One of the ways is to use a simulated annealing algorithm [39] for the network parameters selection: current ( I ON , I OFF , I I ), coupling strength ( s r , s o , s m ), noise amplitude U n and the synchronization effectiveness threshold η th . The algorithm’s main point is the random searching of lectronics , , x FOR PEER REVIEW 12 of 26 the problem II solution with the maximum value of P at some initial interval of parameters, followed by the narrowing of these intervals. Determination of step size is required for direct parameter searching. In our numerical experiment, the current range was determined by the generation of oscillation in a single oscillator. For the oscillator circuit described in Section 2.1, the generation of oscillation was within the feeding current range of 550 µA ≤ I p ≤ 1061 µA. Therefore, the currents ( I ON , I OFF , I I ) varied in this range. We determined the variation steps as δ I p_i = 1 µA. So, there were 512 current steps. For coupling strength s , the variation range was determined by the threshold values of I–V characteristics of a switch. The condition, that is, the integral thermal action on the neuron s Σ should not reduce the threshold voltage of the I–V characteristic of a switch U th below the holding voltage U h , must be met: th h U s U (13) On the basis of Formula (13), for the values of threshold voltages ( U th = 5 V, U h = 1.5 V) and coupling configurations (see Figure 2), the limits of the coupling strength variations are subject to the following conditions: ( s r + s o ·9) < 3.5 V and ( s r + s m ·4) < 3.5 V. This condition is met with the following ranges that we have chosen: s r = 0 ÷ s m = 0 ÷ s o = 0 ÷ U n ≤ 900 µV with the number of steps equal to 12 was chosen for the noise amplitude. For the synchronization effectiveness threshold η th , the variation range was 10 % ≤ η th < 100%, with the number of steps equal to 25 and minimal spacing δ η th = 0.1%. Parameter η th does not belong to the network parameters but rather to the parameters of the algorithm of synchronization order SHR calculation. The value η th strongly affects the results of synchronization and the results of pattern recognition because it is a conditionally chosen characteristic. Identifying its optimal value for the recognition problems is an important step in network tuning and training. Having a lot of network parameters, each parameter’s variant should be calculated in 102 classes together with oscillograms and synchronization values SHR at each stage. A full, direct search of all parameters’ variants would take a lot of time and computational resources. Therefore, we fixed the values U n and η th , and varied currents ( I ON , I OFF , I , I ) and couplings strength ( s r , s o , s m ). The algorithm for searching for the solution of problem II with the maximum value P included the following steps: Preparation step: Setting of constant values of U n and η th Step 1: Random searching of parameters ( I ON , I OFF , I , I , s r , s o , s m ) in the maximal range of their variations and finding the values meeting the maximum value P . The number of searching attempts is 1000. Step 2: Narrowing of the parameters ranges by 5 times with their symmetric distribution in relation to the results of the previous step. The number of searching attempts is 1000. Step 3: Narrowing of the parameters ranges by 25 times with their symmetric distribution in relation to the results of the previous step. The number of searching attempts is 1000. Noise amplitude U n is kept fixed because this is the parameter that is not often controlled in the experiment. It is determined by external and internal circuit noises, and we considered its effect individually in more detail. We used U n = 80 µV, which is close to the experimentally observed value [19]. The threshold value of synchronization effectiveness η th is fixed because it is the main parameter in the synchronization algorithm, and we considered its effect individually in more detail. For the pattern recognition experiments, we used a workstation (Intel Xeon quad core processor, Albuquerque , New Mexico, USA, 4 × 2 GHz, 8GB RAM) running 64-bit Windows Server 2008. CPU time for a single run on one core for the direct search procedure with 1000 search attempts took ~5 h. lectronics , , x FOR PEER REVIEW 13 of 26
3. Results
Following the proposed algorithm, we fixed the values of the following parameters ( U n = 80 µV, η th = 90%) and varied the currents ( I ON , I OFF , I , I ) and coupling strength ( s r , s o , s m ). The distribution of the solutions after the first, second and third steps of training is shown in the diagram (Figure 12), the values of P are on the x-axis, and the corresponding number of solutions N P are on the y-axis. The largest number of solutions corresponds to P = 0 when there is no synchronization at any input class C j . After steps 2–3. the number of solutions with a low value of P decreases and the solutions with a higher value of P appear. Step 3 gives more solutions in the range P = 6–10 in comparison with step 2, although the maximum value P = 14 is similar in both cases. The incorrect solution of training (see Figure 13a), which does not correspond to the problem I-III conditions, occurs rather often. For example, two or more classes may correspond to the same value of SHR , and the probability of such an event at the first step of training is ~55% (the calculation was based on the data in Figure 12), which results in the ambiguous recognition of figures and their classes. Another frequent case of wrong training is the absence of oscillatory synchronization for any input class ( P = 0), which has a probability of ~30% in the first step (see the values in Figure 12). Figure 12.
Diagram of the solution number N P distribution against the value of P for three subsequent steps of the network training algorithm. The total number of attempts at each step is 1000. The maximum value P = 14 was obtained at the following parameters: ( I ON = 722 µA, I OFF = 1034 µA, I = 1020 µA, I = 887 µA, s r =0.10176 V, s o = 0.29202 V, s m = 0.202 V, U n = 80 µV, η th = 90%). N p P Step 1 Step 2 Step 3 lectronics , , x FOR PEER REVIEW 14 of 26 Figure 13.
Training results at various values of current I ON : examples of incorrect solutions of training ( a ) and solutions of problem I ( b ); examples of correct solutions of problem II ( c – e ). The complete list of the model parameters is presented in Supplementary Materials (Data2.txt). . Solution of Problem I
The experiment demonstrated that the solution for problem I ( P = 1) is easily found for some classes C m . For example, Figure 13b shows a neural network that recognizes class C with the corresponding synchronization order SHR = 20:29, and there is no synchronization for all other classes. The probability of any solution with P = 1 during the first step is ~10% (here the total number of attempts is 1000 and the value of N p ~ 100, see Figure 12), and this is the maximum probability of the solution in comparison with other P > 1. Figure 14 shows the distribution of solutions at P = 1 for various values of m (4) after 60 repetitions of the first step of training at the maximum range of all parameter variations. The maximum probability (~4%) can be seen for solutions with sets C and C when all cells of the input pattern are either empty or colored, see Figure 9. Solutions for other m are much rarer, with the probability being two orders lower (~0.03%). Parameters corresponding to the same solution can differ significantly. For example, the system can recognize class C at the input, while at the output SHR would be different for different parameters, however, the problem is still considered solved. The histogram shows that the network can be trained to solve problem I with a certain, predetermined value of m . lectronics , , x FOR PEER REVIEW 15 of 26 Figure 14.
Distribution of the solutions number for problem I ( P = 1) according to the values m (4), calculated for 60 repetitions of the first step of training ( U n = 80 µV, η th =90%). Insets show the classes of input images C m for m = 1, m = 5, m = 102. . Solution of Problem II
The condition of problem II requires that a certain class from set C should correspond to the unique synchronization order from set SHR while there is no output oscillator synchronization for other classes. For example, at the first step of training, N p = 26 solutions were found for P = 2, two of them are shown in Figure 13c,d, and Figure 13e demonstrates the variant of set C and SHR for P = 4. We can find solutions with a random search, for example, for P = 2, however, these solutions would be for different satisfying condition (12). For P = 2, there are C = 5151 solution combinations (where C is the number of combinations 2 out of 102). For P = 14, C is a 17-digit number. Therefore, we indicate the maximum value of P , and do not indicate which solution was found. After all training steps, the maximum value of P reached P = 14. While ONN configuration and training algorithm can be further improved, the purpose of this work is to demonstrate that a multilevel neural network can be implemented and used for pattern recognition. The solution for problem III has not been found yet (see Section 4).
We constructed a 3D graph as shown in Figure 15, to find the dependence of maximum possible P and N p on the noise amplitude value in network U n ( η th = 90%). The values were taken from the first step of training. When the noise increases above U n = 400 µV, the number of solutions N P and the value P sharply decreases. Most of the solution variants are distributed in the range 1 ≤ P ≤ 2. The maximum value for N p corresponds to P = 1 at noise level U n = 400 µV. In this case, the integral value ∑ N P for all values P also has a maximum for U n = 400 µV. Therefore, “stochastic resonance” is present when a certain level of noise induces maximum number of solutions for problems I–II. This may be caused by two differently directed tendencies of N p reduction: the occurrence of “extra” synchronizations with decreasing U n and suppression of the number of synchronizations with increasing U n .
20 40 60 80 1000.1110100100010000 C ={X , X }, X , X C
102 511 ={X }
C ={X } N p mP =1 lectronics , , x FOR PEER REVIEW 16 of 26 Figure 15.
Dependence of solution number N P for various values of P on noise level U n ( η th = 90%). We constructed a 3D graph as shown in Figure 16, to find the dependence of maximum possible P and N p on the value of threshold synchronization effectiveness ( U n = 80 µV). A general tendency for reduction of N p with a decrease in η th below 40% can be seen, and the growth of η th up to η th = 99% on average does not change the values of P and N p . This seems to be caused by the fact that synchronization occurring during problems I and II solution has a high value of effectiveness η > 99%. In this case, reduction of η th just adds “extra” synchronous states, which do not meet the problem conditions. The fact that the maximum possible P reached P = 11 at η th = 30% is interesting, and we observed some resonance for values of P . Maximum P declines as with reduction of η th as with growth of η th , caused by reduction of in the solution number N p . Figure 16.
Dependence of solution number for various values of P on the threshold synchronization effectiveness η th ( U n = 80 µV). lectronics , , x FOR PEER REVIEW 17 of 26 To understand how the output oscillator changes its synchronization relative to the reference oscillator, and what role the layer of processing neurons plays in the process, we conducted the following model experiments. First, we turned off the effect of input oscillators on the output layer, making s o = 0 V. The result was a circuit equivalent to two coupled oscillators, when oscillator No.0 affects oscillator No.10 with a force of s r = 0.3 V. By varying the oscillator currents ( I , I ), we obtained the standard distribution pattern of SHR in the form of Arnold tongue, where the regions of equal synchronization are extended sections in the form of rays (see Figure 17a) with diagonal symmetry. A special feature is a smooth (gradient) increase in SHR from the lower right distribution corner ( I = 1061 µA, I = 500 µA) to the upper left corner ( I = 500 µA, I = 1061 µA). Similar patterns of synchronization distribution in the system of two oscillators were observed by us earlier, for example, in [33]. After adding the effect of the processing layer, when s o = 0.1 V, we get the type of synchronization distribution shown in Figure 17b. The synchronization areas are island type. The layer of input oscillators divides Arnold’s tongue into local areas while the tendency of the gradient increase in synchronization is preserved. If we fix the currents I and I , when s o = 0.29 V, and vary the currents I ON and I OFF , the pattern of the synchronization distribution changes its appearance (see Figure 18a). We observed a chaotic scatter of the synchronization magnitude over the field of parameters, with preservation of the diagonal symmetry and a wide scatter of values within 1–200 µA. As I ON and I OFF values change the frequency of many oscillators in the processing layer at once, it is difficult to predict the tendency of the synchronization distribution. Nevertheless, the range of SHR variations is much smaller than when I and I are varied, which is obviously due to the constancy of the main frequency between the oscillators whose synchronization is measured. With a decrease in the value of s o to s o = 0.1 V, the form of the distribution has a similar appearance, but the range of SHR change is reduced. In analyzing Figures 17 and 18, we can assume that the model annealing method, where we randomly select the system parameters, is searching for optimal combinations of synchronization regions. Due to the large number of regions and their unpredictable order, it is possible to find a solution to problems I–II, and it is not unique, as we observed in previous experiments. ( a ) ( b ) Figure 17.
The synchronization distribution SHR in the regions of currents I and I with ( a ) s o = 0 V and ( b ) s o = 0.1 V. For all cases, I ON = 725 µA, I OFF = 1036 µA, s r = 0.3 V, s m = 0.207 V, η th =90%, U n = 80 µV, and the class C is the input. lectronics , , x FOR PEER REVIEW 18 of 26 ( a ) ( b ) Figure 18.
The synchronization distribution SHR in the regions of currents I ON and I OFF with ( a ) s o = 0.29V and ( b ) s o = 0.1 V. For all cases, I = 1017 µA, I = 891 µA, s r = 0.1036 V, s m = 0.207 V, η th = 90%, U n = 80 µV, and the class C is on the input. The following numerical experiment shows some features of the network operation dynamics. The idea was to find out how the matrix of the processing layer affects the synchronization of oscillators No.0 and No.10. Let us consider two cases. In the first case, s o = 0 V, s r = 0.13 V, I = 650 µA, I = 950 µA, there is only one unidirectional effect of oscillator No.0 on No.10, and the base synchronization has the value SHR b0,10 = 23:12 = 1.917. In the second case, we added random variations of coupling strength ( s o , s m ) and currents ( I ON , I OFF ) at the same values of s r = 0.13, I = 650 µA, I = 950 µA. As a result, after the first step of training we obtained a set of solutions ( C and SHR ) for various P > 0 and positioned all elements of SHR on the plot (see Figure 19). It can be seen that all elements of the sets
SHR are distributed within some interval above the boundary SHR b0,10 . The dispersion is δSHR ~0.6. The same behavior is observed when the current values were interchanged ( I = 950 µA and I = 650 µA). The base synchronization had the value of SHR b0,10 = 8:15 = 0.533, and when other parameters were varied, the elements of SHR diverge upward with δSHR ~ 0.4. At high and similar values of current I = 1058 µA and I = 1058 µA, the base synchronization is SHR b0,10 = 1:1 = 1, and δSHR ~ 0.3. Thus, the matrix of processing layer oscillators only deviates the SHR value upward from the base synchronization value of oscillators No.0 and No.10. The effect of oscillators positioned in the processing layer only increases the frequency of oscillator No.10, but cannot decrease it because of the nature of thermal coupling, which can initiate switching but cannot suppress it. I =650 A, I =950 A, s r =0.13 SHR b0,10 =23:12=1.917 I =1058 A, I =1058 A, s r =0.13 SHR b0,10 =1:1=1 S HR , ( a . u . ) P I =950 A, I =650 A, s r =0.13 SHR b0,10 =8:15=0.533 lectronics , , x FOR PEER REVIEW 19 of 26 Figure 19.
Distribution of synchronization values SHR at random variations of coupling strengths ( s o, s m ) and currents ( I ON , I OFF ). Fixed parameters are shown in the plot. The base sync value (SHR b0,10 ) is shown by a solid line.
SHR b0,10 is the sync value at s o = 0.
4. Discussion
In the previous section, we demonstrated how the simulated annealing algorithm can be used for ONN training. In this algorithm, we used three steps for parameter searching, however, we did not manage to surpass the maximum value of P = 14 (see Figure 12). Therefore, we can assume that the simulated annealing method reaches a solution limit, and to search the network parameters with higher P , it might be necessary to vary other ONN properties, for example, the number of oscillators of the processing layer, the I–V characteristics of oscillators, and the matrix of couplings S . Furthermore, the values of U n could be selected more carefully on the basis of the regularities shown in Figure 15, where the maximum of the solution number N P at a certain noise level is presented. This effect is similar to stochastic resonance phenomena in spike networks [40–42], and requires additional research and analysis. The simulated annealing algorithm might not be an effective training method, and development of more efficient algorithms for this type of network is a separate global problem for future research. For proper network training, a deeper understanding of the physics of the oscillators’ synchronization process is necessary. For example, the results and regularities arising in Figure 19 assist in understanding the process of solution generation. We demonstrated the presence of a boundary for SHR values that is determined by the base synchronization value SHR b0,10 of oscillators No.0 and No.10 and showed that the effect of oscillators’ processing layer leads to SHR ≥ SHR b0 . The selection of currents for the reference oscillators No.0 and No.10 determines the variation range of δSHR , and most likely, the larger this range is, the higher the probability of finding the solutions with high P . If we set the current values I , I , for example, at the boundary of the maximum range I = 1058 µA, I = 1058 µA, then it will be more difficult for the matrix of oscillators in the processing layer to change the dynamics of the output oscillator that operates in a high frequency mode and under the high frequency effect of the reference oscillator. As a result, we observed narrowing of δSHR to ~0.3 (see Figure 19). Figure 18 demonstrates the areas of parameters where SHR practically does not change and where it changes very often, therefore, this may explain the productivity of the annealing training algorithm. In Step 1, we found areas with a large number of solutions, and then in Step 2 and 3, by thorough scanning of these areas, we found solutions with the largest P (see Figure 12). In addition to the configuration and network parameters, the selection of parameter η th , which is used in the synchronization magnitude determination method, is of great importance. As the results and regularities arising in Figure 16 show, too high ( η th ≥ 99.8) and too low ( η th ≤ 20) parameters significantly reduce the number of solutions N P . Nevertheless, we suppose that the usage of η th < 50% is not completely justified, as the parameter η determines the share of a synchronous signal in the total oscillogram [19]. Moreover, the technique of SHR determination may affect the result, and future research aimed at its improvement may lead to positive shifts in this field. The addition of alternative methods of pattern transfer into a network that vary, for example, current I or couplings magnitudes may significantly expand the range of variations of SHR , hence, the maximal attainable value P . The network structure defined by the coupling matrix S (1) was chosen in order to implement an actual VO –based device, therefore, oscillator No.0 affected all other oscillators. This can be easily implemented by arranging VO -switches on one substrate. Nevertheless, the development of an actual device requires specific setting of coupling strengths that depend on many parameters [19] and this is beyond the scope of the current work. Although problem III has not been solved, the presented results and solution for problem II for P = 14 showed that multilevel high order synchronization increases net data throughput of a single neuron and enables the creation of multilevel output cascades of neural networks with high functionality. lectronics , , x FOR PEER REVIEW 20 of 26 The duration of the processed oscillograms was ~2.5 s (for the number of pulses ~3000), providing an error less than 0.2% for calculating η (see Appendix A.2). In the experiment, such duration would significantly slow down the work with the system, since the time of the synchronization calculation would be several seconds. Nevertheless, we hope for a real implementation of this idea on oscillators with a large generation frequency (~ 1–10 MHz), which would significantly reduce the oscillogram duration and the time for the synchronization calculation.
5. Conclusions
The paper presents a new model of ONN with high order harmonics synchronization that recognizes and classifies visual patterns by unique synchronous states, i.e., in accordance with their definite symmetry class. The most outstanding feature of this neural network, in comparison with published variants, is that neurons possess not only bistable properties (the presence or absence of synchronization with the reference neuron) but exhibit multilevel synchronization. The presented model has only one output neuron, nevertheless, variation in the value of high order synchronization SHR allows its usage to classify the input pattern into P classes with the set C . The training implementation of this network for solving cognitive tasks is an interesting possibility, as the network consists only of oscillators and does not use other computational modules. The main purpose of the article was to demonstrate a new method of pattern recognition based on coupled oscillators, where cognitive functions are realized due to the high order synchronization effect. This is a universal effect, so regardless of the type of connection, thermal or any other, the cognitive functions can be realized using coupled oscillators. Thermal communication has no special role in pattern recognition, it is only necessary for organizing the synchronization of oscillators. The development of ONN models with high order synchronization effect offers significant potential for increasing the efficiency of artificial intelligence networks, and the development of their training techniques is an important direction for further research. The possibility of implementing this ONN as not only a program code but also as a separately operating device based on real micro- and nano-electronic self-oscillators would enhance the importance of the obtained results and promote further research in this field. Supplementary Materials:
The following are available online at list of classes: Data1.txt, Data2.txt.
Author Contributions: C onceptualization, A.V.; methodology, A.V., M.B. and P.B.; software, A.V.; validation, P.B.; writing—original draft preparation, A.V., M.B and P.B.; project administration, A.V.
Funding:
This research was supported by the Russian Science Foundation (grant no. 16-19-00135).
Acknowledgments:
The authors express their gratitude to O. Dobrynina and Andrei Rikkiev for the valuable comments in the course of the article translation and revision.
Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A
Appendix A.1. Model Circuit of a Coupled Oscillators-Based Neural Network
The circuit of the modelled ONN is presented in Figure 2. It consists of 11 oscillators numbered i = 0..10. The magnitude of the i -th oscillator effect on the j -th oscillator is determined by the coupling strength s i,j , which is set by matrix S (1). The model circuit of a single oscillator is given in Figure A1. It includes a source of current I p_i , capacitance C connected in parallel with a VO switch and noise source U n_i . The capacity magnitude C was constant C = 1 00 nF, but I p_i and U n_i varied within the following ranges: I p (550 µA ÷ 1061 µA), U n (20 µV ÷ 900 µV). The current through the VO switch and its voltage are defined as I sw_i and U i , respectively. Voltage-current characteristics of the VO switch are shown in Figure A2, where the experimental and model curves are presented. lectronics , , x FOR PEER REVIEW 21 of 26 Figure A1.
Circuit of a single oscillator of a VO -based structure. I is the number of an oscillator; I p_i is a source of current, C is a capacitance, U n_i is a source of noise, I sw_i is the current through the VO switch, U i is the voltage at the switch. Figure A2.
Typical experimental I–V characteristics of a separate switch (Experiment curve) and its model curve (Model curve).
All model switches without coupling have the same I–V characteristics, with stationary natural parameters U th = 5 V (threshold voltage), U h = 1.5 V (holder voltage), U cf = 0.82 V (cutoff voltage). The model curve I sw_i = f ( U i ), which has high-resistance (OFF) and low-resistance (ON) segments with corresponding dynamic resistances R off_i = 9100 Ω and R on_i = 615 Ω, can be presented by the following formula: i ioff_isw_i i i cf_i ion_i , if 1 (OFF)( ) ( ) , if 0 (ON) U flagRI f U U U flagR (A1) where i = 0 … 10 is the oscillator number, flag i denotes the VO switch state 1 (OFF—high resistance), 0 (ON—low resistance). Transitions from one state into another can be expressed by the following algorithm: -3 -3 -3 -3 State=ON
Experiment Model I sw = f ( U ) C u rr en t ( A ) Voltage (V) ND R State=OFF U cf U h U th I th I h lectronics , , x FOR PEER REVIEW 22 of 26 i i h_ii i th_i i flag U Uflag flag U U (A2) where U th_i and U h_i are threshold turn-on and holder voltages of switches (see Figure A3): As mentioned in Section 2.1., the model of thermal coupling is based on reduction of threshold voltage U th by the value s due to thermal effect of other switches. The value s is the thermal coupling strength. The physics of thermal interaction is caused by the generation of a heat wave when the switch is turned on, which propagates and acts (heats) on the surrounding switching structures. In [19], we showed that it is possible to introduce an interaction radius R TC beyond which the induced temperature is less than 0.2 K. Therefore, in the model, we assume that each switch only interacts with surrounding switches that are within the radius of R TC . The heating of structures leads to a decrease in the threshold voltage [19], this change characterizes the value of s , which we call the coupling force. The higher the s , the stronger the effect of one switch on the other. Coupling of oscillators results in changing their threshold turn-on voltages U th_i in regard to the switch state ( flag i ) with which they interact. In general, the magnitude U th i can be presented by the following formula: j th_i th j,i(if 0 (ON)) j flag U U s (A3) where j runs though all values, at which the switch is in turn-on state ( flag i = 0), and U th is the natural turn-on voltage without coupling. This means that all effects of s j,i on i -switch are summed up from all the other j -th switches at the turn-on state. The differential equation that describes the circuit operation (Figure A2) is written as: c_i p_i sw_i sw_i c_i n_i ( ) ( ), were ( ) ( ( ) ( )) dU tC I I t I t f U t U tdt (A4) where I is the number of an oscillator, C is capacitance, U c_i is the capacitor voltage, I sw_i is the current through the switch, U n_i is the noise in the switch, and f ( U ) is the I–V characteristic function (A1). Equations (A4) were numerically calculated with respect to time at regular intervals Δ t = 10 −5 s using implicit Euler method and discrete noise U n ( t ) was generated according to the algorithm U n ( t ) = U n0 ·randn( t ), where U n0 —noise amplitude and randn( t )—normal random numbers generated with zero mean and dispersion equal to 1, realized through the algorithm of uniformly distributed random value transformation [43]. As a result, by setting oscillator current values I p_i , noise amplitude U n0 , and coupling strengths ( s r , s m , s o ) we could calculate oscillograms of current I sw_i ( t ), voltage at the capacity U c_i ( t ) and voltage at switches: i c_i n_i ( ) ( ) ( ) U t U t U t (A5) lectronics , , x FOR PEER REVIEW 23 of 26 Figure A3.
Examples of calculated oscillograms sections (1000 points), voltage U and current I sw_0 , for the oscillator with i = 0 at I p_0 = 1200 µA. Here we do not show other circuit parameters because other oscillators do not affect this one. Figure A4.
Spectrum of a current signal for the oscillogram shown in Figure A3.
Examples of sections (1000 points) of calculated oscillograms for voltage U and current I sw_0 are shown in Figure A3. The spectrum of the current signal has a large number of harmonics, see Figure A4. The main frequency F = 1269 Hz determines the period of current oscillations T ~ 788 µs (because T = 1/ F ). The interaction between the oscillators occurs according to the current signal as it is shown in the model, and has also been observed in an actual experiment, see [19,35], because the current signal is transformed into thermal pulses that spread along the substrate. As a result, the high order synchronization effect can be observed between the coupled oscillators. Appendix A.2. Dependence of SHR and η on the Number of Pulses in the Oscillogram
An important issue in calculating SHR and η according to the method in Section 2.3, with the desired accuracy, is the determination of the number of pulses used for the calculation. In fact, it is necessary to set the minimum duration of the processed oscillogram. This is also important for determining the minimum calculation time for a family of metrics. For example, in a real experiment -3 -3 -3 -3 I s w _0 , A t,s U , V t, s -4 -4 -4 -4 -3 k =7 A m p li t ude , A F, Hz k =1 F =1,269 Hz lectronics , , x FOR PEER REVIEW 24 of 26 to calculate SHR and η , it is required to record the oscillogram, and then make the calculation. In a model experiment, this time is determined by the oscillogram simulation time. Figure A5a shows the dependence of η on the number of oscillogram points with the following system parameters ( I = 1017 µA, I = 891 µA, I ON = 725 µA, I OFF = 1035 µA s r = 0.1036 V, s m = 0.207 V, s o = 0.29298 V, η th =90%, U n = 80 µV, the input image is X ). With an increase in the number of points, the number of pulses in the oscillograms of oscillators No.0 and No.10 grows (see Figure A5b.) ( a ) ( b ) Figure A5.
The dependence of η ( a ) and the number of pulses ( b ) on the number of points on the oscillogram. The dashed line corresponds to 250,000 oscillogram points. With an increase in the number of points (pulses) in the processed oscillogram, the value of η comes to saturation, which is the true value of η . The synchronization value is defined as SHR = 38:37 with the number of points more than 60,000 (when η ≥ η th at η th = 90%); with a smaller number of points, the synchronization is not detected ( η < η th ). In this work, we used 250,000 points (the number of pulses ~3000), which gives an error less than 0.2% in determining η . The duration of the oscillogram was ~2.5 s (Δ t = 10 −5 s). References Freeman, W.J. Spatial properties of an EEG event in the olfactory bulb and cortex.
Electroencephalogr. Clin. Neurophysiol. , , 586–605, doi:10.1016/0013-4694(78)90126-8. 2. von der Malsburg, C. The Correlation Theory of Brain Function. In Models of Neural Networks ; Springer: New York, NY, USA, 1994; pp. 95–119. 3.
Gray, C.M.; Singer, W. Stimulus-specific neuronal oscillations in orientation columns of cat visual cortex.
Proc. Natl. Acad. Sci. USA , , 1698–1702, doi:10.1073/PNAS.86.5.1698. 4. Eckhorn, R.; Bauer, R.; Jordan, W.; Brosch, M.; Kruse, W.; Munk, M.; Reitboeck, H.J. Coherent oscillations: A mechanism of feature linking in the visual cortex?
Biol. Cybern. , , 121–130, doi:10.1007/BF00202899. 5. Burton, S.D.; Ermentrout, G.B.; Urban, N.N. Intrinsic heterogeneity in oscillatory dynamics limits correlation-induced neural synchronization.
J. Neurophysiol. , , 2115–2133, doi:10.1152/jn.00362.2012. 6. White, J.A.; Chow, C.C.; Rit, J.; Soto-Treviño, C.; Kopell, N. Synchronization and Oscillatory Dynamics in Heterogeneous, Mutually Inhibited Neurons.
J. Comput. Neurosci. , , 5–16, doi:10.1023/A:1008841325921. 7. Akam, T.; Kullmann, D.M. Oscillations and Filtering Networks Support Flexible Routing of Information.
Neuron , , 308–320, doi:10.1016/j.neuron.2010.06.019. 8. Kuzmina, M.; Manykin, E.; Surina, I. Oscillatory network with self-organized dynamical connections for synchronization-based image segmentation.
Biosystems , , 43–53, doi:10.1016/J.BIOSYSTEMS.2004.05.005. 9. Corinto, F.; Bonnin, M.; Gilli, M. Weakly Connected Oscillatory Network Models for Associative and Dynamic Memories.
Int. J. Bifurc. Chaos , , 4365–4379, doi:10.1142/S0218127407020014. , % number of points oscillator No.0 oscillator No.10 nu m be r o f pu l s e s number of points lectronics , , x FOR PEER REVIEW 25 of 26 Nikonov, D.E.; Csaba, G.; Porod, W.; Shibata, T.; Voils, D.; Hammerstrom, D.; Young, I.A.; Bourianoff, G.I. Coupled-Oscillator Associative Memory Array Operation for Pattern Recognition.
IEEE J. Explor. Solid-State Comput. Devices Circuits , , 85–93, doi:10.1109/JXCDC.2015.2504049. 11. Hoppensteadt, F.C.; Izhikevich, E.M. Pattern recognition via synchronization in phase-locked loop neural networks.
IEEE Trans. Neural Netw. , , 734–738, doi:10.1109/72.846744. 12. Acebrón, J.A.; Bonilla, L.L.; Pérez Vicente, C.J.; Ritort, F.; Spigler, R. The Kuramoto model: A simple paradigm for synchronization phenomena.
Rev. Mod. Phys. , , 137–185, doi:10.1103/RevModPhys.77.137. 13. Klinshov, V.V.; Nekorkin, V.I. Synchronization of delay-coupled oscillator networks.
Physics-Uspekhi , , 1217–1229, doi:10.3367/UFNe.0183.201312c.1323. 14. Vassilieva, E.; Pinto, G.; Acacio de Barros, J.; Suppes, P. Learning Pattern Recognition Through Quasi-Synchronization of Phase Oscillators.
IEEE Trans. Neural Netw. , , 84–95, doi:10.1109/TNN.2010.2086476. 15. Hagerstrom, A.M.; Murphy, T.E.; Roy, R.; Hövel, P.; Omelchenko, I.; Schöll, E. Experimental observation of chimeras in coupled-map lattices.
Nat. Phys. , , 658–661, doi:10.1038/nphys2372. 16. Vodenicarevic, D.; Locatelli, N.; Abreu Araujo, F.; Grollier, J.; Querlioz, D. A Nanotechnology-Ready Computing Scheme based on a Weakly Coupled Oscillator Network.
Sci. Rep. , , 44772, doi:10.1038/srep44772. 17. Romera, M.; Talatchian, P.; Tsunegi, S.; Abreu Araujo, F.; Cros, V.; Bortolotti, P.; Trastoy, J.; Yakushiji, K.; Fukushima, A.; Kubota, H.; et al. Vowel recognition with four coupled spin-torque nano-oscillators.
Nature , , 230–234, doi:10.1038/s41586-018-0632-y. 18. Pikovsky, A.; Rosenblum, M.; Kurths, J.
Synchronization: A Universal Concept in Nonlinear Sciences ; Cambridge University Press: Cambridge, UK, 2001; ISBN 9780521533522. 19.
Velichko, A.; Belyaev, M.; Putrolaynen, V.; Perminov, V.; Pergament, A. Thermal coupling and effect of subharmonic synchronization in a system of two VO based oscillators. Solid-State Electron. , , 40–49, doi:10.1016/j.sse.2017.12.003. 20. Sakai, J. High-efficiency voltage oscillation in VO planer-type junctions with infinite negative differential resistance. J. Appl. Phys. , , 103708, doi:10.1063/1.2930959. 21. Shukla, N.; Parihar, A.; Freeman, E.; Paik, H.; Stone, G.; Narayanan, V.; Wen, H.; Cai, Z.; Gopalan, V.; Engel-Herbert, R.; et al. Synchronized charge oscillations in correlated electron systems.
Sci. Rep. , , 4964, doi:10.1038/srep04964. 22. Maffezzoni, P.; Daniel, L.; Shukla, N.; Datta, S.; Raychowdhury, A. Modeling and Simulation of Vanadium Dioxide Relaxation Oscillators.
IEEE Trans. Circuits Syst. I Regul. Pap. , , 2207–2215, doi:10.1109/TCSI.2015.2452332. 23. Belyaev, M.A.; Boriskov, P.P.; Velichko, A.A.; Pergament, A.L.; Putrolainen, V.V.; Ryabokon, D.V.; Stefanovich, G.B.; Sysun, V.I.; Khanin, S.D. Switching Channel Development Dynamics in Planar Structures on the Basis of Vanadium Dioxide.
Phys. Solid State , , 447–456, doi:10.1134/s1063783418030046. 24. Boriskov, P.P.; Velichko, A.A.; Pergament, A.L.; Stefanovich, G.B.; Stefanovich, D.G. The effect of electric field on metal-insulator phase transition in vanadium dioxide.
Tech. Phys. Lett. , , 406–408, doi:10.1134/1.1482750. 25. Datta, S.; Shukla, N.; Cotter, M.; Parihar, A.; Raychowdhury, A. Neuro Inspired Computing with Coupled Relaxation Oscillators. In
Proceedings of the 51st Annual Design Automation Conference on Design Automation Conference—DAC ’14 ; ACM Press: New York, NY, USA, 2014; pp. 1–6. 26.
Parihar, A.; Shukla, N.; Datta, S.; Raychowdhury, A. Exploiting Synchronization Properties of Correlated Electron Devices in a Non-Boolean Computing Fabric for Template Matching.
IEEE J. Emerg. Sel. Top. Circuits Syst. , , 450–459, doi:10.1109/JETCAS.2014.2361069. 27. Shukla, N.; Parihar, A.; Cotter, M.; Barth, M.; Li, X.; Chandramoorthy, N.; Paik, H.; Schlom, D.G.; Narayanan, V.; Raychowdhury, A.; et al. Pairwise coupled hybrid vanadium dioxide-MOSFET (HVFET) oscillators for non-boolean associative computing. In ; IEEE: San Francisco, CA, USA, 2014; p. 28.7.1–28.7.4. 28.
Velichko, A.; Belyaev, M.; Putrolaynen, V.; Pergament, A.; Perminov, V. Switching dynamics of single and coupled VO -based oscillators as elements of neural networks. Int. J. Mod. Phys. B , , 1650261, doi:10.1142/S0217979216502611. lectronics , , x FOR PEER REVIEW 26 of 26 Ghosh, S. Generation of high-frequency power oscillation by astable mode arcing with SCR switched inductor.
IEEE J. Solid-State Circuits , , 269–271, doi:10.1109/JSSC.1984.1052131. 30. Chen, C.L.; Mathews, R.H.; Mahoney, L.J.; Calawa, S.D.; Sage, J.P.; Molvar, K.M.; Parker, C.D.; Maki, P.A.; Sollner, T.C.L.G. Resonant-tunneling-diode relaxation oscillator.
Solid-State Electron. , , 1853–1856, doi:10.1016/S0038-1101(00)00105-2. 31. Sharma, A.A.; Bain, J.A.; Weldon, J.A. Phase Coupling and Control of Oxide-Based Oscillators for Neuromorphic Computing.
IEEE J. Explor. Solid-State Comput. Devices Circuits , , 58–66, doi:10.1109/JXCDC.2015.2448417. 32. Locatelli, N.; Cros, V.; Grollier, J. Spin-torque building blocks.
Nat. Mater. , , 11–20, doi:10.1038/nmat3823. 33. Velichko, A.; Belyaev, M.; Putrolaynen, V.; Boriskov, P. New Method of the Pattern Storage and Recognition in Oscillatory Neural Networks Based on Resistive Switches.
Electronics , , 266, doi:10.3390/electronics7100266. 34. Velichko, A.A.; Stefanovich, G.B.; Pergament, A.L.; Boriskov, P.P. Deterministic noise in vanadium dioxide based structures.
Tech. Phys. Lett. , , 435–437, doi:10.1134/1.1579818. 35. Velichko, A.; Belyaev, M.; Putrolaynen, V.; Perminov, V.; Pergament, A. Modeling of thermal coupling in VO -based oscillatory neural networks. Solid-State Electron. , , 8–14, doi:10.1016/j.sse.2017.09.014. 36. Velichko, A.; Putrolaynen, V.; Belyaev, M. Effects of Higher Order and Long-Range Synchronizations for Classification and Computing in Oscillator-Based Spiking Neural Networks. arXiv , arXiv:1804.03395. 37. Belyaev, M.; Velichko, A.; Putrolaynen, V.; Perminov, V.; Pergament, A. Electron beam modification of vanadium dioxide oscillators.
Phys. Status Solidi Curr. Top. Solid State Phys. , , doi:10.1002/pssc.201600236. 38. Reljan-Delaney, M.; Wall, J. Solving the linearly inseparable XOR problem with spiking neural networks. In ; IEEE: London, UK, 2017; pp. 701–705. 39.
Callan, R.
The Essence of Neural Networks ; Prentice Hall Europe: Upper Saddle River, NJ, USA, 1999; ISBN 013908732X. 40.
Collins, J.J.; Chow, C.C.; Imhoff, T.T. Stochastic resonance without tuning.
Nature , , 236–238, doi:10.1038/376236a0. 41. Kawaguchi, M.; Mino, H.; Durand, D.M. Stochastic Resonance Can Enhance Information Transmission in Neural Networks.
IEEE Trans. Biomed. Eng. , , 1950–1958, doi:10.1109/TBME.2011.2126571. 42. Uzuntarla, M.; Cressman, J.R.; Ozer, M.; Barreto, E. Dynamical structure underlying inverse stochastic resonance and its implications.
Phys. Rev. E , , 42712, doi:10.1103/PhysRevE.88.042712. 43. Parzen, E.