The theoretical foundations of CNNs were laid in 1988 by L. Chua (Chua and Yang, 1988a; Chua and Yang, 1988b). In the most general case, a CNN consists of a spatially discrete collection of locally coupled k ^th-order continuous-time processing elements, called cells, arranged at regular positions within a l-dimensional lattice. The architecture of a small two-dimensional CNN with M = 6 rows and N = 6 columns is presented in Figure 1A, under the assumption that each cell C ( i , j ) – i ∈ { 1 , … , M } , j ∈ { 1 , … , N } – is physically coupled to its 8 adjacent neighbors only ³ . Each cell is assigned a state, an input, an output, as well as a threshold. The rate of change of the state of a cell is influenced by the inputs and outputs of its 8 adjacent neighbors, as well as by its own input and output, as respectively sketched through eight brown directed segments and through one magenta directed loop in Figure 1A for the processing element located where the 3^rd row crosses the 4^th column. The block diagram in plot (b) of Figure 1 illustrates once more the key factors affecting the dynamical behaviour of a CNN cell. The neighbors’ inputs (outputs) are modulated by feedforward (feedback) synaptic weights before accessing the cell to affect the time evolution of its state, similarly as it occurs in biological neural networks. The cell C ( i , j ) of a standard time- and space-invariant two-dimensional CNN (Chua and Roska, 2002) is implemented by the circuit of Figure 2, where the computing core is mathematically described by ⁴ ( i ∈ { 1 , … , M } , j ∈ { 1 , … , N } ) d v x i , j d t = − v x i , j C x i , j ⋅ R x i , j + z ⋅ I C x i , j + 1 C x i , j ⋅ ∑ k = − r k = r ∑ l = − r l = r ( i a k , l + i b k , l ) , (1) in case it is physically coupled to its 8 adjacent neighbors only ⁵ i.e., r = 1. With reference to the circuit of Figure 2, the main physical quantity within the input stage, the computing core, and the output stage respectively are the input voltage v u i , j , the voltage v x i , j across a capacitor with capacitance C x i , j ℝ > 0 , expressing the state, and the output voltage v y i , j . Focusing on the output stage, the latter physical quantity is defined as v y i , j = R y i , j ⋅ i f i , j , (2) where R y i , j ℝ > 0 is the resistance of a passive linear resistor, whereas i f i , j is the current of a voltage-controlled source, featuring the piecewise linear expression i f i , j ≜ i f ( v x i , j ) = g l i n ⋅ | v x i , j + v s a t | − | v x i , j − v s a t | 2 , (3) generally known as standard nonlinearity, where g l i n and v s a t are positive parameters with units Ω − 1 and V, respectively. Importantly, the piecewise-linear standard nonlinearity identifies three domains, specifically the negative saturation region v x i , j < − v s a t , the linear region | v x i , j | ≤ v s a t , and the positive saturation region v x i , j > + v s a t , where the cell output voltage v y i , j attains the negative saturation voltage − v s a t , is a linear function of the state v x i , j , and attains the positive saturation voltage v s a t , respectively. Inspecting now the computing core in the cell circuit of Figure 2, i z , defined as i z ≜ z ⋅ I , (4) where z is a dimensionless parameter referred to as threshold in CNN theory (Chua and Roska, 2002), while I is a coefficient with positive 1 A value, symbolizes the threshold current. Further, R x i , j ℝ > 0 stands for the resistance of a passive linear resistor ⁶ , while, importantly, i a k , l ( i b k , l ), defined as i a k , l ≜ a k , l ⋅ v y i + k , j + l (5) i b k , l ≜ b k , l ⋅ v u i + k , j + l , (6) where a k , l ( b k , l ), with k , l ∈ { − 1,0,1 } , is known as feedback (feedforward) synaptic weight in CNN theory (Chua and Roska, 2002), represents the feedback (feedforward) synaptic current due to the neighboring cell C ( i + k , j + l ) , and constituting one of the 18 contributions to the cell capacitor current i x i , j enclosed within the round brackets to the right of the double sum in Eq. 1. It is worth mentioning that, the CNN mathematical description reported in Eq. 1 is known as Chua-Yang model (Chua and Yang, 1988a; Chua and Yang, 1988b). Despite an alternative mathematical description, known as Full-Range model (Vázquez et al., 1993), facilitates the hardware implementation of the CNN paradigm by restricting the set of allowable values for the cells’ states, this paper adopts the original Chua-Yang mathematical description for pedagogical purposes.

Typically, the right hand side of the CNN cell state Eq. 1 may be recast as d v x i , j d t = i g i , j + i w i , j C x , (7) where i g i , j , the so-called Internal Driving Point (DP) Component is a function of the cell state, being defined as i g i , j ≜ i g ( v x i , j ) = { − v x i , j R x − a 0,0 ⋅ R y ⋅ g l i n ⋅ v s a t if v x i , j < − v s a t , (8) ( a 0,0 ⋅ R y ⋅ g l i n − 1 R x ) ⋅ v x i , j if | v x i , j | ≤ v s a t , (9) − v x i , j R x + a 0,0 ⋅ R y ⋅ g l i n ⋅ v s a t if v x i , j > + v s a t . (10) while i w i , j , known as offset current, mostly accounts for the coupling effects, being expressed by i w i , j = i w ( { v u i + k , j + l } , { v y i + k , j + l } ) ≜ z ⋅ I + b 0,0 ⋅ v u i , j + ∑ k , l = − 1 ( k , l ) ≠ ( 0,0 ) ( a k , l ⋅ v y i + k , j + l + b k , l ⋅ v u i + k , j + l ) , (11) where { v u i + k , j + l } ( { v y i + k , j + l } ) denotes the set of input (output) voltages of all the processing elements in the 3 × 3 neighborhood of the cell C ( i , j ) . Nineteen are the key parameters, which define the image processing operation, which a CNN may accomplish, for a predefined input/initial condition combination, and under suitable boundary conditions, specifically the threshold z, the feedback synaptic weights { a − 1 , − 1 , a − 1,0 , a − 1 , + 1 , a 0 , − 1 , a 0,0 , a 0 , + 1 , a + 1 , − 1 , a + 1,0 , a + 1 , + 1 } , and the+feedforward synaptic weights { b − 1 , − 1 , b − 1,0 , b − 1 , + 1 , b 0 , − 1 , b 0,0 , b 0 , + 1 , b + 1 , − 1 , b + 1,0 , b + 1 , + 1 } . Given the crucial role that this 19-number set plays on the dynamical evolution of the cells’ states, it is generally known as gene in CNN theory. A gene defines the set of rules, which apply concurrently in the neighborhood of each cell, allowing the standard space-invariant CNN to carry out a given computation.

Remark 1. A CNN may be used to carry out any computing task. The calculations are based upon the analogue dynamics of the cell states. As transients vanish, depending on cell parameter settings, inputs and initial conditions, each capacitor voltage may tend toward an equilibrium, or converge to an oscillatory waveform, which may be periodic, quasi-periodic, or even chaotic. A CNN may then process the information, inserted as cell inputs and/or initial conditions, in two different forms i.e., producing computation results in the form of equilibria or waves, hence the names CNN equilibria-based computing or CNN wave computing attributed to the respective operating principle. In this manuscript the attention is focused on CNNs computing via equilibria ⁷ .

 Let us elucidate the classical method to design a CNN, so that it may execute a fundamental image processing task ⁸ , carrying out the result of the computation in the cell equilibria.

 A rigorous technique to synthesize the gene of a standard CNN, so as to allow the successful execution of a given equilibria-based computing task, even in the presence of deviations of some cell circuit parameters from their nominal values, was introduced in (Zarándy, 2003), and comprehensively presented in (Itoh and Chua, 2003). Before summarizing the line-of-thought at the basis of this classical methodology, let us provide a brief overview of the nonlinear dynamics of a CNN processing element.

The first aspect to consider for the synthesis of a suitable CNN gene is the choice of a proper value for the self-feedback synaptic weight a 0,0 . As will be clarified through qualitative sketches below, this parameter crucially influences the directed Internal DP Characteristic, consisting of the arrowed locus of the rate of change of the state v ˙ x i , j vs. the state v x i , j itself under i w i , j = 0   A. As anticipated in section 2.2 in the context of memristors, this type of graphical representation is typically referred to as State Dynamic Route (SDR). In the upper (lower) half of the v x i , j – v ˙ x i , j plane, where v ˙ x i , j > ( < ) 0   V ⋅ s⁻¹, arrows, supeimposed over the Internal DP Characteristic, point toward the east (west) to indicate an increase (a decrease) in the state v x i , j over time. The intersections between this v ˙ x i , j – v x i , j locus and the horizontal v x i , j -axis, marked as filled (hollow) circles, denote the stable (unstable) equilibria of the cell state equation under null offset current. All in all, fixing the value for a 0,0 unequivocally determines the dynamical behaviour of the cell state v x i , j from any initial condition v x i , j , 0 under zero offset current. With reference to ⁹ Figure 3, plot (a) shows how a 0,0 affects the shape of the locus of the state rate of change v ˙ x i , j vs. the state v x i , j itself under i w i , j = 0   A . As anticipated in section 2.2 in the context of memristors, a family of directed loci of this kind, one for each value of a control parameter (in this case a 0,0 ), is called DRM, here more specifically referred to as cell DRM. Varying a 0,0 from − ∞ to + ∞ , the CNN processing element may exhibit three distinct stability characters:

• If a 0,0 < R x − 1 ⋅ g l i n − 1 ⋅ R y − 1 (see the green and brown curves, associated to non-null negative and null self-feedback values, respectively) the cell is monostable with one and only one GAS equilibrium at v ¯ x i , j = 0   V.

The ordinates of the two breakpoints of the three-segment ¹⁰ piecewise-linear Internal DP Characteristic, located at ( − v s a t , − ( a 0,0 ⋅ R y ⋅ g l i n − R x − 1 ) ⋅ v s a t ) , and at ( + v s a t , + ( a 0,0 ⋅ R y ⋅ g l i n − R x − 1 ) ⋅ v s a t ) , respectively, are of significant importance in the analysis of the Shifted ¹¹ DP Characteristic (Chua and Roska, 2002), i.e. the locus of the state rate of change v ˙ x i , j vs. the state v x i , j itself under non-null offset current. First of all, it is important to point out that, under specific hypotheses, including fixed values for all input voltages and thresholds, a standard space-invariant ¹² CNN is completely stable (Chua and Roska, 2002), in the sense that the state v x i , j of each cell C ( i , j ) converges asymptotically toward an equilibrium point v ¯ x i , j , which, in general, depends upon the initial conditon v x i , j , 0 . Moreover, for a completely stable standard space-invariant CNN, according to the bistability criterion (Chua and Roska, 2002), in case the condition a 0,0 > R x − 1 ⋅ g l i n − 1 ⋅ R y − 1 (12) holds true, the slope of the Internal DP Characteristic in the linear region–refer to Eq. 9 – is positive, and, consequently, the stable equilibria of each cell under i w i , j ≠ 0   A are found to lie in either of the two saturation regions of the standard nonlinearity of Eq. 3, as will be elucidated, shortly, through the analysis of the resulting Family of Shifted DP Characteristics, implying that, given the form of the standard nonlinearity Eq. 3, the final value for the output voltage of any processing element is one of the two saturation levels in the set { − v s a t , v s a t } . This is highly desirable for image processing, where, as will be shown through an illustrative example shortly, a CNN equilibria-based computation is typically visualized in the form of a binary output image, coding the final values of the cell output voltages. Importantly, the output voltage of each processing element will attain its final value i.e., one of the two possible saturation levels, already at a finite time instant, let us denote it as t i , j ( s ) , at which its state v x i , j enters the saturation region, which hosts the equilibrium point v ¯ x i , j . For all t ≥ t ( s ) ≜ max 1 ≤ i ≤ M , 1 ≤ j ≤ N { t i , j ( s ) } the CNN may be considered at steady state as far as its cell output voltages are concerned ¹³ . It follows that, irrespective of the location of the CNN cell C ( i , j ) under analysis, the offset current i w i , j , which, in general, changes over time during transients, keeps a fixed value for all t ≥ t ( s ) . This is of significant importance, since from this time instant onward, the Shifted DP Characteristic will no longer bounce, as, on the other hand, may be the case during transients, facilitating the study of the dynamic behaviour of the state v x i , j from any initial condition v x i , j , 0 .

With reference to the viewgraphs of Figure 3B, neglecting the marginal cases, three are the possible equilibria configurations, which a cell C ( i , j ) may feature under the bistable criterion hypothesis for i w i , j ≠ 0   A.

• If i w i , j < − ( a 0,0 ⋅ R y ⋅ g l i n − R x − 1 ) ⋅ v s a t (see the blue curve) the cell turns into a monostable dynamical system, featuring one and only one globally asymptotically stable (GAS) equilibrium in the negative saturation region at v ¯ x i , j = − R x ⋅ ( a 0,0 ⋅ R y ⋅ g l i n ⋅ v s a t − i w i , j ) .

The standard method (Zarándy, 2003; Itoh and Chua, 2003) to synthesize a suitable CNN gene for the execution of a given computing task is based upon the set up and later solution of an ad-hoc set of inequalities, expressed in terms of unknown gene parameters, ensuring a robust accomplishment of the computing task of interest. The functionalities of a given uncoupled ¹⁴ standard space-invariant completely stable CNN, satisfying the hypotheses of the bistability criterion, are dictated by a set of local rules, establishing the asymptotic value for the state v x i , j ( ∞ ) ≡ v ¯ x i , j , and, correspondingly, the steady-state output voltage v y i , j ( t i , j ( s ) ) of each cell C ( i , j ) , depending upon inputs, and initial conditions of all the processing elements within its 3 × 3 neighbourhood.

For the reasons motivated above, in CNN designs for image processing applications, it is common to set the numerical value for the self-feedback synaptic weight a 0,0 so as to guarantee the satisfaction of inequality Eq. 12. Typically, to facilitate the CNN design process, the structure of the gene under synthesis is simplified as much as possible, given the degree of complexity of the computing task, which the cellular array is expected to execute, and/or the values of some of the elements from the 19-number parameter set are assumed to be known. The family of all the possible invariable arrowed Shifted DP Characteristics, which a cell may admit for all t > t ( s ) for any value, which the offset current may assume, then, under the hypothesis of each of the rules, is then examined, so as to identify the worst-case scenario, where deviations of cell circuit parameters from their nominal values are most likely to induce a change in the cell equilibria configuration ¹⁵ , and, consequently, the emergence of CNN computation errors. The next step is to write down an inequality, establishing a constraint for the offset current, and ensuring that, in such critical scenario, v ˙ x i , j is negative (positive) at the specified initial condition v x i , j , 0 , so that the state v x i , j would approach a desired equilibrium v ¯ x i , j in the negative (positive) saturation region. Repeating this procedure for each rule allows to derive an inequality set (IS), whose solutions may be determined through numerical techniques, or, in case the number of unknowns is small, via a geometry-based approach. The particular values assigned to the unknowns, allowing to program the CNN with a suitable gene, are then selected to endow the computation with the highest degree of tolerance against parameter variation.

The aim of this section is to synthesize the gene of a standard space-invariant non-autonomous uncoupled CNN so that the resulting nonlinear dynamic array is able to extract the edges from an input binary image. Next, the classical CNN design methodology, briefly described earlier, will be applied to the cell ODE (1) in order to achieve this purpose. The three local rules, which each cell should obey, are reported in Table 1, where n B defines how many of the 8 adjacent neighbors feature a positive one V-valued input voltage. The first rule establishes that, if the cell C ( i , j ) features an input voltage v u i , j equal to − 1   V , its output voltage v y i , j ( ∞ ) at equilibrium is found to attain the negative saturation voltage − v s a t irrespective of the value of n B . Plot (a) in Figure 4 depicts a possible 3 × 3 pattern around a white pixel at row i and column j in the input binary image under rule 1. Here 3 of the 8 neighboring pixels are black. The pixel in the position ( i , j ) of the output binary image is white at equilibrium, as depicted on the bottom of the input pattern. The second (third) rule imposes that, in case the cell C ( i , j ) features an input voltage v u i , j equal to + 1   V , its output voltage v y i , j ( ∞ ) at equilibrium is found to attain the negative (positive) saturation voltage − v s a t in case n B is exactly equal to 8 (is less or equal to 7). Plot (b) ((c)) in Figure 4 depicts the only (a) possible 3 × 3 pattern around a black pixel at row i and column j in the input binary image under rule 2 (3). Here all (4) of the 8 neighboring pixels are black. The pixel in the position ( i , j ) of the output binary image is white (black) at equilibrium, as depicted on the bottom of the input pattern.

As clarified by Figure 5A, the CNN under design is expected to extract the edges from an input binary image, visualizing them in the output binary image at steady state. Since the CNN is meant to be uncoupled, the offset current from Eq. 11 reduces to ¹⁶ i w i , j = z ⋅ I + b 0,0 ⋅ v u i , j + ∑ k , l = − 1 ( k , l ) ≠ ( 0,0 ) 1 b k , l ⋅ v u i + k , j + 1 . (13)

Assuming that all the feedforward synaptic weights, with the exclusion of b 0,0 , are identical one to the other, namely b k , l = b for all k , l ∈ { − 1,0 , + 1 } such that ( k , l ) ≠ ( 0,0 ) , indicating how many, among the 8 neighbours of the cell C ( i , j ) , feature a negative one V-valued input voltage through the variable n W , and noting that n B + n W = 8 , the formula Eq. 13 for i w i , j reduces to i w i , j ( v u i , j , n B ) = z ⋅ I + b 0,0 ⋅ v u i , j + b ⋅ ( 2 ⋅ n B − 8 ) V , (14) where the argument reveals the dependence of the offset current upon v u i , j and n B . Assuming that a 0,0 ∈ ℝ > 0 and b ∈ ℝ < 0 are given parameters, the only two unknowns for the specification of a suitable gene are then b 0,0 and z. Figure 5B shows the directed Internal DP Characteristic ¹⁷ . The state Eq. 1 under i w i , j = 0   A admits two locally stable equilibria, located one in the negative saturation region, namely v ¯ x i , j = − a 0,0 ⋅ R y ⋅ g l i n ⋅ R x ⋅ v s a t , and one in the positive saturation region, specifically v ¯ x i , j = a 0,0 ⋅ R y ⋅ g l i n ⋅ R x ⋅ v s a t , and separated by an unstable equilibrium, i.e. v ¯ x i , j = 0   V, positioned in the linear region.

Let us set the initial condition v x i , j ( 0 ) of the cell ODE (1) to + 1   V. Following the line-of-thought inspiring the classical CNN design methodologies, discussed in the seminal papers (Zarándy, 2003) and (Itoh and Chua, 2003), and briefly reviewed above, let us now examine the Family of arrowed Shifted DP Characteristics, which a processing element may admit in all scenarios, which may possibly emerge under the hypothesis of each of the three rules from Table 1.

In order to fulfill rule 1, where v u i , j = − 1   V , a condition should be enforced to ensure that the maximum value, which the offset current may ever attain i.e., max 0 ≤ n B ≤ 8 { i w i , j ( − 1   V , n B ) } = i w i , j ( − 1   V , 0 ) , is smaller than the ordinate − ( a 0,0 ⋅ R y ⋅ g l i n − G x ) ⋅ v s a t of the left breakpoint of the v ˙ x i , j – v x i , j piecewise linear characteristic of Figure 5B. This would guarantee a negative sign for the ordinate of the right breakpoint of the resulting v ˙ x i , j – v x i , j piecewise-linear characteristic, as is the case for the arrowed blue locus in Figure 6A, illustrating the dynamic route followed by the cell state for i w i , j ( v u i , j , n B ) = − 2 ⋅ ( a 0,0 ⋅ R y ⋅ g l i n − G x ) ⋅ v s a t , where v u i , j = − 1   V, and n B = 0 , under the parameter setting, reported in the caption of Figure 5B. As a result, for all possible n B values in { 0,1,2,3,4,5,6,7,8 } , the CNN cell would be monostable, and v x i , j would decrease monotonically over time from the initial condition v x i , j ( 0 ) toward an equilibrium, i.e. v ¯ x i , j = ( − a 0,0 ⋅ R y ⋅ g l i n ⋅ v s a t + i w i , j ( − 1   V , n B ) ) ⋅ R x , located in the negative saturation region, as established by rule 1. Therefore, the first EDGE CNN design constraint sets an upper bound for the maximum offset current, according to i w i , j ( − 1   V , 0 ) < − ( a 0,0 ⋅ R y ⋅ g l i n − G x ) ⋅ v s a t . (15) It is worth pointing out that the farther away from the horizontal axis, within the plane lower half, would be positioned the right breakpoint of the v ˙ x i , j – v x i , j piecewise-linear characteristic in the worst-case scenario from rule 1, and the more robust would be the EDGE CNN design ¹⁸ .

Let us now derive the condition allowing the CNN to apply rule 2 from Table 1 in the sphere of influence of any processing cell C ( i , j ) , which features, as each of its eight neigbours, a positive one V-valued input voltage. Since the expected cell steady-state output voltage v y i , j ( t i , j ( s ) ) is once again − 1   V , as in rule 1, Figure 6A can be reused to work out a suitable inequality for rule 2 under v u i , j = + 1   V and n B = 8 . The second EDGE CNN design condition, ensuring that the state v x i , j of a cell C ( i , j ) with v u i , j = + 1   V and n B = 8 would asymptotically approach an equilibrium, specifically v ¯ x i , j = ( − a 0,0 ⋅ R y ⋅ g l i n ⋅ v s a t + i w i , j ( + 1   V , 8 ) ) ⋅ R x , located in the negative saturation region, is then similar to the inequality Eq. 5, reading as i w i , j ( + 1   V , 8 ) < − ( a 0,0 ⋅ R y ⋅ g l i n − G x ) ⋅ v s a t . (16) In order for the CNN under design to apply rule 3 from Table 1 in the 3 × 3 neighbourhood of each processing element, which features a positive one V-valued input voltage, and is physically coupled to at least one neighbour with a negative one V-valued input voltage, the minimum value, which the offset current may ever attain, namely min 1 ≤ n B ≤ 7 { i w i , j ( + 1   V , n B ) } = i w i , j ( + 1   V , 7 ) should be larger than the ordinate − ( a 0,0 ⋅ R y ⋅ g l i n − G x ) ⋅ v s a t of the left breakpoint of the v ˙ x i , j – v x i , j piecewise-linear characteristic of Figure 5B. This would ensure a positive sign for the ordinate of the right breakpoint of the resulting v ˙ x i , j vs. v x i , j piece-wize linear characteristic, as is the case for the arrowed blue locus in Figure 6B, illustrating the dynamic route of the cell state for i w i , j ( v u i , j , n B ) = − 0.5 ⋅ ( a 0,0 ⋅ R y ⋅ g l i n − G x ) ⋅ v s a t , where v u i , j = + 1   V, and n B = 7 , under the parameter setting, reported in the caption of Figure 5B. Consequently, for all admissible n B values in { 0,1,2,3,4,5,6,7 } , the cell state v x i , j would monotonically increase over time toward an equilibrium, i.e. v ¯ x i , j = ( a 0,0 ⋅ R y ⋅ g l i n ⋅ v s a t + i w i , j ( v u i , j , n B ) ) ⋅ R x , located in the positive saturation region, as required in rule 3. The third EDGE CNN design inequality is then establishing a lower bound for the minimum offset current, i.e. ¹⁹ i w i , j ( + 1   V , 7 ) > − ( a 0,0 ⋅ R y ⋅ g l i n − G x ) ⋅ v s a t . (17) For a robust CNN design the right breakpoint of the v ˙ x i , j – v x i , j piecewise-linear characteristic of a cell C ( i , j ) in the worst-case scenario from rule 3 should be positioned as farther away as possible from the horizontal axis within the plane upper half ²⁰ .

For the parameter setting reported in the caption of Figure 5, the three inequalities Eqs. 15–17 are solved through a geometric approach on the z– b 0,0 parameter plane, as shown in Figure 7, where the green region visualizes the set of admissible solutions. For the specification of a suitable gene, guaranteeing the expected EDGE CNN functionality even in the presence of some small deviation of either of the two parameters z and b 0,0 from their nominal values, it is adviceable to choose a particular solution ( z * , b 0,0 * ) , whose graphical point-based representation on the parameter plane features an adequate distance from the boundaries of the green region, as indicated by means of an asterisk marker in Figure 7. The gene, synthesized in this section, allows the CNN to extract edges from an input binary image, as displayed in plot (a) of Figure 5.

Since each of their processing elements interacts simultaneously with the respective neighbors, CNNs may process multi-variate signals in a massively parallel fashion, as crucially necessary in time-critical application fields, such as industry process control, electronic surveillance, medical augmented reality, and IoT smart sensing. In order to harness more efficiently the bio-inspired operating principles of these nonlinear dynamic arrays, which make them a suitable mathematical framework for modeling neural systems, Chua and Roska proposed an innovative computer, called CNN Universal Machine (CNN-UM) (Roska, 1993), to implement their signal processing paradigm. The CNN-UM, fabricated in various forms over the years through the well-established CMOS technology ²¹ (Vázquez et al., 2018), consists of an array of locally coupled computing units, each of which is endowed with data storage blocks, which allow to distribute the memory across the cellular array, endowing the computing machine with a truly non-von Neumann architecture, and to reconfigure the array so as to solve any computation problem. Thanks to their massively parallel computing power, CNN-UM hardware realizations (Vázquez et al., 2018) may process images at rates as high as 30,000 frames per second. Considering that, furthermore, a universal cellular array may be physically realized within the IC area of a single chip (Vázquez et al., 2018), CNNs are particularly suitable for the development of miniaturized IoT technical systems, in which the integration between a matrix of sensing elements and a network of locally coupled computing units with local stored programmability on board enables information processing at the same location, where data detection takes place. A major problem, which prevents to widen the applicability scope of this class of sensor-processor arrays, is the limited degree of complexity of the dynamical phenomena, which may possibly emerge within their physical media, due to the simplicity of the input-output behaviours of the electrical components employed in the CNN-UM constitutive blocks. Thanks to their extremely rich dynamics, memristors may be adopted in novel designs of cellular computing arrays so as to extend significantly the spectrum of asymptotic spatio-temporal behaviours, which purely CMOS CNNs may currently exhibit. Another critical issue, which affects the performance of technical systems, combining sensing and processing functionalities on the same physical platform, is due to the rather low spatial resolution of state-of-the-art CNN-UM hardware realizations, originating from the presence of spacious data storage units within their computing units, as discussed earlier. This limits the maximum number of sensing and processing elements, which may be paired ²² one-to-one within the available IC area of these IoT commercial products (Toshiba Ltd., 2012). The adoption of memristive devices, endowed with memprocessing capabilities, may allow to obviate the inclusion of additional memory banks within the IC design of each CNN-UM computing unit, allowing to shrink considerably its size, and enabling the future realization of sensor-processor arrays with unprecedented spatial resolution, of great appeal to the IoT industry, nowadays. In this respect, it is timely to commence investigations aimed to explore the functionalities of Memristor CNNs (M-CNNs). In general, introducing memristors in the circuit implementation of a CNN processing element ²³ increases the order of its ODE model, calling for the development of a new theory to investigate the operating principles of the resulting nonlinear dynamic array, and to program its gene to allow the accomplishment of a pre-defined memcomputing task. The theoretical foundations of M-CNNs shall be discussed in the section to follow.

The M-CNN cell C ( i , j ) of Figure 8 may be described by the following pair of first-order coupled ODEs ²⁵ ( i ∈ { 1 , … , M } , j ∈ { 1 , … , N } ): d x m i , j d t = g ( x m i , j , v x i , j ) , and (18) d v x i , j d t = i ˜ g i , j + i w i , j C x . (19) The first ODE Eq. 18 governs the time evolution of the state x m i , j of the first-order nonvolatile resistance switching memory ℳ x , which the M-CNN cell C ( i , j ) acccommodates, according to an enhanced variant of a voltage-controlled memristor model, originally formulated by Pershin and Di Ventra (Pershin et al., 2009), and capable to capture the switching kinetics of real memristor devices (Jo et al., 2009), as discussed in Pershin and Di Ventra (2011). The model of the resistance switching memory in the cell C ( i , j ) is a first-order element from the class of generic memristors, defined by the DAE set d x m i , j d t = g ( x m i , j , v m i , j ) , (20) i m i , j = G ( x m i , j ) ⋅ v m i , j . (21) Note that, within the processing element C ( i , j ) , the memristor voltage v m i , j coincides with the capacitor voltage v x i , j , thus the expression for the memristor current i m i , j in Eq. 21 reduces to i m i , j = G ( x m i , j ) ⋅ v x i , j . (22) The state evolution function g ( x m i , j , v m i , j ) and the memductance function G ( x m i , j ) in the Pershin and Di Ventra model of the memristor in the M-CNN cell C ( i , j ) are respectively expressed by g ( x m i , j , v m i , j ) = κ ( v x i , j ) ⋅ ( step ( v m i , j ) ⋅ f + ( p ) ( x m i , j ) , + step ( − v m i , j ) ⋅ f − ( p ) ( x m i , j ) ) , and (23) G ( x m i , j ) = 1 x m i , j . (24) The memristor state x m i , j , representing the device memristance, is constrained to lie at all times within the closed set D ≜ [ x o n , x o f f ] , where x o n and x o f f denote the lowest and highest possible device resistances, respectively. With reference to Eq. 23, step ( ⋅ ) stands for the Heaviside function, while κ ( v m i , j ) is a piecewise-linear nonlinearity of the form κ ( v m i , j ) = − β ⋅ v m i , j + β − α 2 ⋅ ( | v m i , j + V t | − | v m i , j − V t | ) , (25) where α ∈ ℝ > 0 and β ∈ ℝ > 0 are coefficients, measured in units Ω ⋅ V − 1 ⋅ s − 1 , denoting the smaller and larger slopes of the characteristic for | v m i , j | ≤ V t and | v m i , j | > V t , respectively, where V t ∈ ℝ > 0 represents the memristor switching threshold voltage. Figure 9A depicts the k ( v m i , j ) – v m i , j chapacteristic for the parameter setting reported in its caption.

Since the memristor state existence domain D is finite, the state evolution function in Eq. 23 is endowed with boundary conditions, which ensure that x m i , j never decreases below (increases above) its lowest (largest) possible value under v m i , j > ( < )   0   V . In order to facilitate the numerical simulation of the memristor DAE set, we reformulate the boundary conditions as compared to their original definition ²⁶ in the Pershin and Di Ventra model (Pershin et al., 2009), adopting continuous and differentiable functions, inspired to Biolek’s window (Biolek et al., 2009), and reading as f + ( p ) ( x m i , j ) = 1 − ( x m i , j − x o n x o f f − x o n − 1 ) 2 ⋅ p ,   and (26) f − ( p ) ( x m i , j ) = 1 − ( x m i , j − x o n x o f f − x o n ) 2 ⋅ p , (27) where p ∈ ℕ > 0 controls the decay rate of the window function Eqs. 26, 27 as x m i , j approaches x o n ( x o f f ). As graphically illustrated in plot (b) ((c)) of Figure 9 for the parameter configuration provided in its caption, the window function f +   ( − ) ( p ) ( x m i , j ) in Eqs. 26, 27 enforces the memory state evolution function Eq. 23 to feature a zero, and, consequently, the memristor ODE Eq. 18 to admit an equilibrium at x ¯ m i , j = x o n   ( o f f ) under positive (negative) values of the capacitor voltage. Since the memory state ODE Eq. 18, with evolution function expressed by Eq. 23, is of first-order, the classical DRM graphical tool (Chua, 2018a) may be applied to investigate the memristor nonlinear dynamics. The DRM of the modified Pershin and Di Ventra memristor model is illustrated in Figure 10A for the parameter arrangement defined in its caption. The DC value V m i , j , assigned to the voltage falling across the resistance switching memory, parametrizes the family of memristor SDRs. Within the family of x ˙ m i , j vs. x m i , j loci, the characteristic obtained for V m i , j = 0   V, known as POP, provides hints on the nonvolatile memory capability of the circuit element. On the basis of the memristor model under focus, the POP is a segment of the x m i , j axis lying between x o n and x o f f . Each of the points on this segment–shown in black in Figure 10A–represents a stable but not asymptotically stable equilibrium (Strogatz, 2000) for the ODE Eq. 18 with state evolution function Eq. 23. Particularly, the existence of a continuum of equilibria, namely x ¯ m i , j ∈ D , for the memristor state equation under zero input clearly reveals that the resistance switching device is an analogue non-volatile memory. Any value for x m i , j within its existence domain D is a possible state, which the memristor may store, from the time at which the power is turned off, till the time at which a new voltage stimulus is applied between its terminals. With regard to the x ˙ m i , j – x m i , j loci, associated to nonzero values for V m i , j , in Figure 10A, the device asymptotically approaches the fully off (fully on) resistive equilibrium state x ¯ m i , j = x o f f   ( o n ) in case any negative (positive) DC voltage is applied continually between its terminals, as indicated by the arrow superimposed on each blue (red) characteristic, which dictates a memory state rate of change increasing monotonically with | V m i , j | . Irrespective of the negative (positive) DC value assigned to the memristor voltage, the upper (lower) bound in the memristor state existence domain D is found to be a globally asymptotically stable equilibrium for the ODE (18) with state evolution function Eq. 23. For the very same parameter setting, Figure 10B demonstrates now the smooth periodic change, which the state x m i , j of the memristor in the cell C ( i , j ) undergoes over each cycle of a sinusoidal voltage appearing between its terminals, and mathematically expressed by v m i , j = v ^ m i , j ⋅ sin ( 2 ⋅ π ⋅ f m i , j ⋅ t ) , where v ^ m i , j = 2   V and f m i , j = 100   Hz. Clearly, at any given time instant, the cell is effectively a second-order dynamical system with degrees of freedom provided one by the memristor state and one by the capacitor voltage, which is also illustrated in plot (b) of Figure 10. Visualizing the memristor current flowing through the memristor as a result of the capacitor voltage from plot (b) vs. the capacitor voltage itself, the resulting pinched hysteresis loop, shown in Figure 10C, gives further evidence for the analogue dynamic behaviour of the cell memristor. The second M-CNN cell ODE Eq. 19 governs the time evolution of the cell capacitor voltage v x i , j within the memcomputing core of the circuit of Figure 8. Its right hand side is identical as in the ODE Eq. 1 dictating the rate of change of the capacitor voltage within the computing core of the cell of the standard time- and space-invariant two-dimensional CNN discussed in section 2.1, except for the presence of an additional addend, resulting from the current through the memristor. It follows that the expression for the offset current i w i , j of the memristive processing element of Figure 8 is still given by Eq. 11, while, using Eq. 22 to express the current through the memristor, the formula for the cell Internal DP Component g ˜ i , j features the new form i ˜ g i , j ≜ i ˜ g ( x m i , j , v x i , j ) = − v x i , j x m i , j − v x i , j R x − a 0,0 ⋅ R y ⋅ g l i n ⋅ v s a t if v x i , j < − v s a t , (28) ( a 0,0 ⋅ R y ⋅ g l i n − 1 R x − 1 x m i , j ) ⋅ v x i , j if | v x i , j | ≤ v s a t , (29) − v x i , j x m i , j − v x i , j R x + a 0,0 ⋅ R y ⋅ g l i n ⋅ v s a t if v x i , j > + v s a t . (30) in which Eqs. 22, 24 were employed to model the cell memristor current i m i , j and the memductance function G ( x m i , j ) , respectively. It is worth to note that the number of variables in the argument of i g i , j is a signature for the order of the cell, as can be inferred by comparings Eqs. 8–10 and Eqs. 28–30. The classical cell DRM technique (Chua, 2018a), reviewed in section 2.2, and adopted for the analysis and synthesis of standard CNNs with first-order processing elements, is applicable to dynamical systems with one degree of freedom only. As a result, the development of a systematic procedure to investigate and design M-CNNs with second-order memristive processing elements calls for a preliminary generalization of the DRM graphic tool. Drawing inspiration from the phase portrait concept from the theory of nonlinear dynamics (Strogatz, 2000), the next section introduces a new system-theoretic notion, which we name Second-Order DRM (DRM₂), enabling the investigation of the memcomputing capabilities of cellular nonlinear arrays with second-order memristive cells.