A Model of an Oscillatory Neural Network with Multilevel Neurons for Pattern Recognition and Computing

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.

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,6,7] and neuromorphic algorithms of image processing based on synchronization [8,9,10,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.

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,21,22,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,26,27,28]. Relaxation oscillators with high order synchronization effects can be realized by using electric coupling [25,26,27,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

2.1. Oscillator Circuit and Method of ONN Organization

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.

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.

2.2. ONN Structure

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).

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:

$$S=\left[\begin{array}{ccc}{s}_{0,0}& \dots & s_{0,10}\\ \dots & \dots & \dots \\ s_{10,0}& \dots & s_{10,10}\end{array}\right]=\left[\begin{array}{ccccccccccc}0& s_{\mathrm{r}}& s_{\mathrm{r}}& s_{\mathrm{r}}& s_{\mathrm{r}}& s_{\mathrm{r}}& s_{\mathrm{r}}& s_{\mathrm{r}}& s_{\mathrm{r}}& s_{\mathrm{r}}& s_{\mathrm{r}}\\ 0& 0& s_{\mathrm{m}}& 0& s_{\mathrm{m}}& 0& 0& 0& 0& 0& s_{\mathrm{o}}\\ 0& s_{\mathrm{m}}& 0& s_m& 0& s_{\mathrm{m}}& 0& 0& 0& 0& s_{\mathrm{o}}\\ 0& 0& s_{\mathrm{m}}& 0& 0& 0& s_{\mathrm{m}}& 0& 0& 0& s_{\mathrm{o}}\\ 0& s_{\mathrm{m}}& 0& 0& 0& s_{\mathrm{m}}& 0& s_{\mathrm{m}}& 0& 0& s_{\mathrm{o}}\\ 0& 0& s_{\mathrm{m}}& 0& s_{\mathrm{m}}& 0& s_{\mathrm{m}}& 0& s_{\mathrm{m}}& 0& s_{\mathrm{o}}\\ 0& 0& 0& s_{\mathrm{m}}& 0& s_{\mathrm{m}}& 0& 0& 0& s_{\mathrm{m}}& s_{\mathrm{o}}\\ 0& 0& 0& 0& s_{\mathrm{m}}& 0& 0& 0& s_{\mathrm{m}}& 0& s_{\mathrm{o}}\\ 0& 0& 0& 0& 0& s_{\mathrm{m}}& 0& s_{\mathrm{m}}& 0& s_{\mathrm{m}}& s_{\mathrm{o}}\\ 0& 0& 0& 0& 0& 0& s_{\mathrm{m}}& 0& s_{\mathrm{m}}& 0& s_{\mathrm{o}}\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]$$

(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_10,i = 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_0,i = 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:

$$I_{\mathrm{p}\_\mathrm{i}}=\{\begin{array}{l}{I}_{\mathrm{ON}},\mathrm{if}{x}_{\mathrm{i}}=1\\ I_{\mathrm{OFF}},\mathrm{if}{x}_{\mathrm{i}}=0\end{array}$$

(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_0,10 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).

2.3. Method of Synchronization Order Definition

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

$${\mathrm{SHR}}_{\mathrm{i},\mathrm{j}}=k_{\mathrm{i}}:k_{\mathrm{j}}$$

(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_{\mathrm{S}}=k_{\mathrm{i}}\cdot F_{\mathrm{i}}^0=k_{\mathrm{j}}\cdot F_{\mathrm{j}}^0$$

(4)

where F_i⁰, F_j⁰ are frequencies of main harmonics.

As an example, Figure 4 shows the qualitative spectra of two oscillators at 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 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).

Therefore, SHR_i,j may be estimated using a phase-locking method:

$${\mathrm{SHR}}_{\mathrm{i},\mathrm{j}}=M_{\mathrm{j}}:M_{\mathrm{i}}$$

(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_i⁰ = M_j·T_j⁰, F_i⁰ = 1/T_i⁰, F_j⁰ = 1/T_j⁰).

In general, especially when a system behaves erratically, synchronization periods differ and spread in T^z_s ≠ T^z+1_s and the values of M_i and M_j may change within one oscillogram (see Figure 7).

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:

$${\mathrm{SHR}}_{\mathrm{i},\mathrm{j}}=M_{\mathrm{j}}:M_{\mathrm{i}},\mathrm{if}P=P_{\max }(M_{\mathrm{j}}:M_{\mathrm{i}})$$

(6)

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:

$$P(M_{\mathrm{j}}:M_{\mathrm{i}})=100\%\cdot NP(M_{\mathrm{j}}:M_{\mathrm{i}})\cdot M_{\mathrm{i}}/N_{\mathrm{i}}$$

(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).

For SHR_i,j, the parameter of synchronization effectiveness η is defined as the maximum probability P_max(M_j:M_i):

$$\eta =P_{\max }(M_{\mathrm{j}}:M_{\mathrm{i}})$$

(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,9,10,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

$$\mathrm{signals}\mathrm{are}\{\begin{array}{l}\mathrm{synchronized},\mathrm{if}\eta \ge {\eta }_{\mathrm{th}}\\ \mathrm{not}\mathrm{synchronized},\mathrm{if}\eta <{\eta }_{\mathrm{th}}\end{array}$$

(9)

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_0,10:

$${\mathrm{SHR}}_{0,10}=k_0:k_{10}~k_0/k_{10}$$

(10)

The same value of SHR_0,10 can be expressed in several ways: as a ratio, a simple fraction or a real number, for example, SHR_0,10 = 10:3 = 10/3 = 3.33. Later, we will use this property to present the results more vividly.

Parameter SHR_0,10 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_0, I₁₀), coupling strength constants s (s_r, s_m, s_o), noise amplitude U_n and η_th varied.

2.4. Pattern Classifier Implementation and Problem Definition

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, Figure 2 and Figure 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).

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°).

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).

Figures from one class impose the same effect on the neural network and have the same output value SHR_0,10. 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_0,10 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_0,10(X₅) = SHR_0,10(X₂₆₀) = SHR_0,10(X₃₂₀) = SHR_0,10(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:

• Synchronization of oscillators No.0 and No.10 with the corresponding value of SHR_0,10 and η > η_th, exists only for one specific class C_j with number j = m out of 102 classes:

  $${\mathrm{SHR}}_{0,10}{(\mathrm{C}}_{\mathrm{j}})=\{\begin{array}{l}{k}_0:k_{10}\mathrm{and}\eta \ge {\eta }_{\mathrm{th}}\mathrm{if}j=m\\ \mathrm{no}\mathrm{syncronysation}\mathrm{and}\eta <{\eta }_{\mathrm{th}}\mathrm{if}j\ne m\end{array}$$

  (11)
  Here we have to show the solutions of this problem with various values of m.
• 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⁽¹⁾_0,10, SHR⁽²⁾_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_0,10. By analogy with (4) the problem may be expressed as:

  $${\mathrm{SHR}}_{0,10}{(\mathrm{C}}_{\mathrm{j}})=\{\begin{array}{l}{\mathrm{SHR}}_{0,10}^{(1)}\mathrm{and}\eta \ge {\eta }_{\mathrm{th}}{\mathrm{if}\mathrm{C}}_{\mathrm{j}}\in \mathbf{C}\mathrm{and}j=Z_1,\\ {\mathrm{SHR}}_{0,10}^{(2)}\mathrm{and}\eta \ge {\eta }_{\mathrm{th}}{\mathrm{if}\mathrm{C}}_{\mathrm{j}}\in \mathbf{C}\mathrm{and}j=Z_2,\\ .............................................................................\\ {\mathrm{SHR}}_{0,10}^{\left(\mathrm{P}\right)}\mathrm{and}\eta \ge {\eta }_{\mathrm{th}}{\mathrm{if}\mathrm{C}}_{\mathrm{j}}\in \mathbf{C}\mathrm{and}j=Z_{\mathrm{P}},\\ \mathrm{no}\mathrm{syncronysation}\mathrm{and}\eta <{\eta }_{\mathrm{th}}{\mathrm{if}\mathrm{C}}_{\mathrm{j}}\notin \mathbf{C}\end{array}$$

  (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 figures belonging to three different classes takes place. Synchronization is not realized for all other classes and η < η_th.

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_0,10. Therefore, the set of input classes C = {C₁, C₂ … C₁₀₂} transfers into a set of synchronous states SHR = {SHR⁽¹⁾_0,10, SHR⁽²⁾_0,10 … SHR⁽¹⁰²⁾_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_0,10 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.

2.5. Technique for ONN Training

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_0, 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 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_0, 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:

$$U_{\mathrm{th}}-s_{\mathsf{\Sigma }}>U_{\mathrm{h}}$$

(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 ÷ 0.2 V, s_m = 0 ÷ 0.5 V, s_o = 0 ÷ 0.3 V. The variation step was chosen as 0.1% of the range magnitude.

The range 20 µV ≤ 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_0,10 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_0,10 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, NM, 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.

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_0,10, 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).

3.1. 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_0,10 = 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_0,10 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.

3.2. 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.

3.3. Solution of Problem III

The solution for problem III has not been found yet (see Section 4).

3.4. Study of the Noise Effect on the Training Results

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.

3.5. Examination of the Synchronization Threshold on the Training Result

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.

3.6. Study of the Dynamics of the Neural Network

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_0,10 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_0,10 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_0,10 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_0,10 change is reduced.

In analyzing Figure 17 and Figure 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.

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^b_0,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_0,10 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^b_0,10. The dispersion is δSHR_0,10 ~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^b_0,10 = 8:15 = 0.533, and when other parameters were varied, the elements of SHR diverge upward with δSHR_0,10 ~ 0.4. At high and similar values of current I₀ = 1058 µA and I₁₀ = 1058 µA, the base synchronization is SHR^b_0,10 = 1:1 = 1, and δSHR_0,10 ~ 0.3.

Thus, the matrix of processing layer oscillators only deviates the SHR_0,10 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.

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,41,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^b_0,10 of oscillators No.0 and No.10 and showed that the effect of oscillators’ processing layer leads to SHR_0,10 ≥ SHR^b₀. The selection of currents for the reference oscillators No.0 and No.10 determines the variation range of δSHR_0,10, 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_0,10 to ~0.3 (see Figure 19). Figure 18 demonstrates the areas of parameters where SHR_0,10 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_0,10, 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.

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_0,10 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.