2010 12th International Workshop on Cellular Nanoscale
Networks and their Applications (CNNA)

Leaning Theory of Cellular Neural Networks
based on Covariance Structural Analysis
Mamoru Tanaka

Hisashi Aomori

Yoshifumi Nishio

Sophia University

Tokyo University of Science

The University of Tokushima

Kumiko Oshima and Martin Hasler
Ecole Polytechnique Federale de Lausanne(EPFL)

Abstract—This paper describes a learning theory of the CNN
based on the covariance structure analysis using new numericalintegral methods. In general, a Cellular Neural Network (CNN)
is defined as a local connected circuit which has continuous statevariables x ∈ Rn . The importance is in that the piece-wise linear
function of the CNN has a linear region |x| ≤ 1 for x ∈ x
because the learning method can be constructed only in linear
state and measurement equations, and because the linear region
can be quantized from the continuous variable x to the multilevel quantized variablef (x) by each 1-bit Σδ modulator whichis corresponding to a spiking neuron model. That is, our purpose
is to determine the weight parameters θ in the connection
matrices A, B, C, D, T and e by the machine learning method
for equilibrium points of the CNN states equation ẋ = 0. The
covariance structure for the equilibrium point to the linear region
will be constructed based on extended Chua’s CNN theorem tohave symmetric edges for aij = aji and asymmetric one-wayedge aij = 0 for aji = 0 for A-matrix A = [aij ].

I. I NTRODUCTION
Recently, neural machine learning methods have been usedas data mining to acquire the important information frommassive amount of data and to predict future. That is, thedata mining is to construct ’model’ for classification, prediction, I/O mapping and association by using machine learningalgorithm. However, the data mining has not been combined
to a circuit model based on the CNN [2], because the CNN
learning methods have not been proposed. Also, from circuit
viewpoints, it is very important to construct various spiking
neuron circuits because the connections will be dynamically ona large reconfigurable FPGA device. The circuit will be used in
the field of neural prosthesis. Though various spiking neuronmodels [3] which are simplified Hodgkin-Huxley models havebeen proposed, the structure has not been built from thegeneral learning. If each cell in the CNN can be changed to aΣδ modulator [1], the CNN with learning will be a superior
spiking neuron model to a brain computer. The problem is
how we construct a generalized Σδ CNN from data. This
paper describes a machine learning of the generalized CNN
based on the covariance structure analysis which includes
factor analysis. The explicit or implicit numerical integral
methods are used for the learning. In the latter case, ateach learning iteration step for quasi-Newton method, theinverse Hessian matrix approximation for the matrix including
second-order partial derivatives is updated by using DavidonFletcher-Powell (DFP) or Groyden-Fletcher-Goldfad-Shanno

978-1-4244-6678-8/10/$26.00 ©2010 IEEE

(BFGS) method. By using both of Backward Euler method
and DFP(or BFGS) method, the solution for stiff systems in
the parameter space can be obtained. That is, our purpose is todetermine the weight parameters θ in the connection matrices
A, B, C, D, T and e by the machine learning method forequilibrium points in linear region of the CNN states equationẋ = 0. The Backward Euler method, one of methods forsolving a nonlinear differential equation, is used in CNN
parameter space. The important point is that this method
can solve stiff systems which have large difference among
eigenvalues. The resulted parameters are used as weights onedges in the cellular signal flow graph (SFG) which works as
a CNN circuit model with V-I transformation devices and as aprediction model for unknown input data. Our simulation forthe model shows good result. This machine learning algorithm
will be in SDP(Sophia Dynamics Program) which has beenused as nonlinear circuit simulator based on Relaxation-basedalgorithm with Fast Automatic Differentiation. The SFG will
be constructed from netlist of the SDP in future.
II. CNN S TATE AND M EASUREMENT E QUATIONS
The observed variable y is a visible information data
which has been obtained from real human behavior, naturalenvironment and so on. The mean μ and the real covariancematrix S = E((y − μ)(y − μ) ) ∈ Rl×l are calculated byusing the observed variables y.
Let x ∈ Rn ,u ∈ Rm be the state(endogenous latent) andinput(exogenous latent) vectors respectively, then the cellular
structural equation is expressed byẋ = −x + Af (x) + Bu + T(1)

where A ∈ Rn×n ,B ∈ Rn×m are coefficient weight matrices
which express connections between the state variables, andT ∈ Rn is the dsturbance variables for the state variables.
Each nonlinear function of f (x) ∈ Rn is a piece-wise linear
function as
f (x) =

1(|x + 1| − |x − 1|).
2

(2)

The cellular measurement equation should be used to express the casual and/or resulted relations between the observed(endogenous and/or exogenous) variables y ∈ Rl and thelatent variables x (u). The measurement equation is givenbyy = μ + Cf (x) + Du + e(3)
where C ∈ Rl×n , D ∈ Rl×m are coefficient weight matrices
which express connections between the latent variables x,uand the observed variables y for its mean μ, and e ∈ Rl isthe error variables for the observed variables.
For example, the matrix A includes a weight element aij on
the edge from a causal cell Cj to the corresponding resulted
cell Ci . Generally, all matrices are sparse.
Our purpose is to make a general program which candetermine the parameters of A, B, C, D, T and e by machinelearning method for equilibrium points of the CNN statesequation ẋ = 0. It is very important that the model constructedby the proposed learning method becomes a CNN circuit
because of high speed prediction and association. Also, it isvery important that the multiplication in Af (x) and Cf (x)becomes integer summation by the use of Σδ modulations ineach of which the nonlinear function f (x) ∈ f (x) can beconstructed by the sum of local dynamical outputs from 1bit quantizer (spiking neuron). The nonlinear function basedon the Σδ modulation type of CNN is described equivalently
from theoretical viewpoints as the following multi-quantizedfunction:⎧
1 (x ≥ ξ0 )⎪⎨ 1
g(x) (−ξ0 ≤ x ≤ ξ0 )(4)
f (x) =⎪⎩ ξ
−1 (x ≤ −ξ0 )
where

⎧
⎪
⎨ Δ([ x + 1 ]) (l is odd)
Δ 2
(5)g(x) =x1⎪⎩ Δ([ ] + ) (l is even)
Δ
2
2ξ
(6)
Δ=
l−1
The parameters Δ and l are quantized step size and thenumber of quantized levels, respectively. The parameter ξ0is the value x to give the definition area of multi-steps
function 1ξ g(x) and the parameter ξ is the value x to satisfy|f (x)| = 1.
The important problems are how to learn the matrices suchthat we have equilibrium point in the linear region |x| ≤ 1(ξ =1) corresponding to Chua’s theorem aii < 1 ( aii =  or 0)
for multi-level AD conversion region and how to construct thegeneralized CNN structure.
III. O PTIMIZATION
A. SDP Program: Minimization of Fit Function
Let z be the model standardized vector, then it is given byz = y − μ.

(7)

If the diagonal elements aii of the matrix A for CNN with
Lyapunov energy function are ,the attractor of dynamics for
the CNN state equation will be in equilibrium point on thelinear region and the equations are expressed for the l × l,n × n or (l + n) × (l + n) unit matrix I as follows:


z
=
x
⎛
⎞


 

e
0 C
z
I 0 D ⎝T ⎠ (8)
+
0 A
f (x)
0 I B
u
The network is not symmetric. Therefore, for asymmetric(causal and resulted) connection, if aij = 0, we set aji = 0 for
its equilibrium. There is no equilibrium problem for symmetric(covariance) connections according to Chua’s theorem[2]. Andif aii =   1, it seems that each state variable x is in
the equilibrium point of linear region as |x| ≤ 1 to become
f (x) = x. This will be proved. In the linear region, if we
define as
E=

0

−1
0 C
A0 = I −
0 A


I 0 D
B0 =
,
0 I B
it is derived that

I


⎛

⎞
e
z = EA0 B 0 ⎝ T ⎠ ,u(9)where (e, T , u) is an exogenous(source) vector.
Let Σ ∈ Rl×l be the covariance matrix E(zz  ) for the statevariables in linear region of the piece-wise linear function,then it is derived theoretically from the cellular structural andmeasurement equations as follow:
Σ = EA0 B 0 Φ0 B 0 A0 E (10)where Φ0 defined as the covariance matrix of the vector(e, T , u) is a (l + n + m) × (l + n + m) diagonal block
matrix by assumption as
⎞
⎛
0
0
ee
0 ⎠.
Φ0 = ⎝ 0 T T 
0
0
uu


0 C
To use in (14), we define that A! =.
0 A
By using the real covariance matrix S ∈ Rl×l and aparameter element vector θ ∈ Rp , the fit function is definedas1
fc (θ) = tr((S − Σ)S −1 )2 .
(11)
2
where tr(X) means the trace of a matrix. In order to minimizec (θ )
the fit function, the function g(θ) = ∂f∂θ
is changed toi∂fc (θ) ∂Σ(θ)∂fc (θ)
= tr
.
(12)
∂θi
∂Σ(θ) ∂θi

In this equation, we derive that

B. SDP Program: Nolinear Solution∂fc= S −1 (Σ − S)S −1 .
∂ΣFor the element φij of the matrix Φ0 , we derive as∂f ∂EA0 B 0 Φ0 B 0 A0 E 
∂fc
= tr
∂φij
∂Σ
∂φij


∂f
= tr
EA0 B 0 (I ij + I ji )B 0 A0 E 
∂Σ


∂f
EA0 B 0 (I ij + I ji )= tr B 0 A0 E 
∂Σ


∂f



EA0 B 0×2=
B 0 A0 E∂Σij

The function (11) should be minimized such that the parameters at that time are determined. So we want to solve theequation:g(θ) = 0(15)where g(θ) = ∂fc .
∂θHowever, the convergence depends on initial value when
the nonlinear equation (15) is solved by an iterative solution
method. In order to escape the initial problem, we solve thefollowing equation by using Backward Euler method.
θ̇ = −g(θ)

where I ij represents a matrix in which only (i, j)-element is
1 and others are 0. For the matrix Φ0 ,∂fc
∂fc
EA0 B 0 .
= 2B 0 A0 E ∂Φ0∂ΣFor the element b0ij of the matrix B 0 , we derive as
∂fc ∂EA0 B 0 Φ0 B 0 A0 E 
∂fc
= tr
∂b0ij
∂Σ
∂b0ij


∂f
∂B 0 Φ0 B 0  
= tr
EA0
A0 E
∂Σ
∂b0ij


∂f

EA0 (Iij Φ0 B 0 + B 0 Φ0 Iji= tr A0 E 
)
∂Σ


∂f
EA0 B 0 Φ0× 2.
=
A0 E ∂Σij
For the matirx B 0 ,∂fc∂f
EA0 B 0 Φ0 .
= 2A0 E ∂B 0∂Σ
For the element αij of the matrix A! , considering

(13)

dA0
dA!
= −A0 (−
)A0 , tr (P Q) = tr (QP ) ,
dαij
dαij
we derive as
∂fc
∂αij




∂fc ∂EA0 B 0 Φ0 B 0 A0 E 
∂Σ
∂αij


∂fc
∂A0 B 0 Φ0 B 0 A0 
= tr
E(
)E
∂Σ
∂αij
∂fc
= tr(
E(A0 Iij A0 B 0 Φ0 B 0 A0
∂Σ

= tr

The Backward Euler method is an implicit numerical integral
method for solving stiff systems. This method approximates
the solution at virtual time tk+1 = tk + h by solving theimplicit equation:
θ k+1 = θ k − hg(θ k+1 )=∂fc
EA0 B 0 Φ0 B 0 A0A0 E then the Newton method is described as
(n+1)

∂Σ

θ k+1 = (n) θ k+1 − (∂F((n) θ k+1 ) −1 (n)) F( θ k+1 ). (19)∂ (n) θ k+1

The inverse matrix computation of the Jacobi matrix (I +(n)h ∂g((n) θ k+1 ) )−1 is not available or expensive. Then the ap∂θ k+1
proximation technique is used. The Hessian matrix H(θ n ) ∈Rp×p which is correspoding to the Jacobi matrix is a secondorder partial derivatives of the function fc with respect toθ n . It is important to use the inverse of the Hessian matrix .
However, since the Hessian leads to algorithmic and computational complexities, an approximation technique of the inverse
Hessian is often used. We use Davidon-Fletcher-Powell (DFP)or Broyden-Fletcher-Goldfab-Shanno (BFGS) method which
is one of quasi-Newton methods. The update formula is asfollows:
DFP method:
Hn+1 = Hn +



(17)

where the gradient vector g(θ k ) is evaluated at θ k .
Since this equation(17) may be nonlinear, solving it ingeneral requires an iterative solution method. An implicit
method requires the solution of the nonlinear equation at eachtime step. For each step of the Backward Euler method, we
use the quasi-Newton method.
Let
F(θ k+1 ) = θ k+1 − θ k + hg(θ k+1 ),(18)+A0 B 0 Φ0 B 0 A0 Iji A0 )E  )(16)H dd Hn
rr
− n

rd
d Hn d(20)

BFGS method:ij

For the matrix A! ,∂fc
∂fc
EA0 B 0 Φ0 B 0 A0 .
= 2A0 E 
∂A!
∂Σ

× 2.

Hn+1 = Hn + (1 +rd Hn + Hn drd Hn d rr
)−(21)r d r d
d r

where
(14)

(n)

(n+1)

(n)

r = −αH n F (θ k+1 ) d = F (θ k+1 ) − F (θ k+1 ).

+

IV. C IRCUIT M ODEL
A circuit model which is constructed by the learning algorithm becomes a general CNN circuit described asState Equationdxij (t)1C=−
xij (t) +A(i, j; k, l)ykl (t)dtRx
C(k,l)∈Nr (i,j)+

1  dyij dxij (t)yij (t)Rx
dxij dt
(i,j)

−



B(i, j; k, l)(i,j) (k,l)

(i,j)

(26)
we can change the A-element AΩ (i, j; k, l) to the original
A-element. Then we can prove the Lyapunov stability as

B(i, j; k, l)ukl + T

C(k,l)∈Nr (i,j)(22)
Output1(|xij + 1| − |xij − 1|)(23)2
where a cell C(i, j) is placed in pixel P (i, j) in 2D-plane.
Definition 1: A connection between a cell C(k, l) and a cellC(i, j) in the CNN is called symmetry for A1 (i, j; k, l) =A1 (k, l; i, j) and one-directional asymmetry for that if
A0 (i, j; k, l) = 0 then A0 (k, l; i, j) = 0.
For the type of each branch connection, we define eachelement of A-template by using a binary parameter Ω suchthat:

2A(i, j; k, l)if Ω = 0(24)
AΩ (i, j; k, l) =A(i, j; k, l)if Ω = 1

dE(t)
dt

= −

 dyij dxij (t)(i,j)

yij (t) =where Ω = 0 and Ω = 1 mean asymmetrical and symmetricalconnections respectively.
Definition 2: if any connection is Ω = 0 or Ω = 1, the CNN
is called a data-mining CNN.
All Cell have different templates in the data-mining CNN.
The weights depends on the learning for the structure determined by a data-mining analyzer. The structure is determinedto satisfy the conditions of the definitions. The data-mining
CNN must give a stable state in linear region [-1,1] of thepiece-wise linear function in parallelism because it is a circuit
model for the data mining. The dynamics will be converged
to an equilibrium point according to the following theorem.
Theorem 1: The data-mining CNN energy function which
is defined as
 1{Ω + (1 − Ω)}AΩ (i, j; k, l)yij (t)ykl (t)E(t) = −2(i,j) (k,l)+

1 yij (t)2
2Rx
(i,j)

−



B(i, j; k, l)yij (t)ukl −(i,j) (k,l)



T yij (t).

(25)

(i,j)

Taking it consideration as
dE(t)
=
dt

1
dyij dxij (t)−ykl (t){Ω + (1 − Ω) }AΩ (i, j; k, l)2dxij dt(i,j) (k,l)

 dyij dxij (t)dyij dxij (t)ukl (t) −,Tdxij dt
dxij dt

dxij

dt


·[

A(i, j; k, l)ykl (t)−C(k,l)∈Nr (i,j)

+1yij (t)Rx

B(i, j; k, l)ukl +T ]C(k,l)∈Nr (i,j)=−

 dxij (t)·dt|xij |<1

·[A(i, j; k, l)ykl (t) −C(k,l)∈Nr (i,j)+

1xij (t)Rx

B(i, j; k, l)ukl + T ]C(k,l)∈Nr (i,j)

dE(t)
=−
dt




C

|xij |<1

dxij (t)
dt

(27)

2
≤0

(28)

Thorem 2: If
1
(29)
Rx
then, each cell converges to the linear region:
The proof should be done such as that in the paper[2].
A(i, j; i, j) <V. C ONCLUSION
A novel cellular learning analysis method of covariancestructure was proposed. It is very important that we build a- ζ in which
structure such that the degree of freedom l(l+1)
2
ζ is the number of variables to be obtained must be positive.
ACKNOWLEDGMENT
This research by corporation of Sophia, Tokushima andEPFL Universities is supported by the fund of Open ResearchCenter Project from MEXT of Japanese Government (20072011).
R EFERENCES[1] H. Aomori, T. Otake, N. Takahashi, and M. Tanaka,” Sigma-detal cellularneural network for 2D modulation,” Neural Networks, vol. 21, no. 2?3,
pp. 349-357, Mar. 2008.
[2] L.O.Chua and L.Yang, ”Cellular neural networks:thorery,” IEEE Trans.
Circuits Syst., vol, no.10, pp. 1257-1272, Oct. 1988.
[3] E.M. Izhikevich,”Which model to use for cortical spiking neurons?,” IEEE
Trans. on Neural Networks, vol. 15, pp. 1063-1070,2004.