Application of Gradient Optimization Methodsin Defining Neural Dynamics

1. Introduction and background

Recurrent neural networks (RNNs) are important class of algorithms for computing matrix (generalized) inverses. These algorithms are used to find the solutions of the matrix equations or to minimize certain nonlinear matrix functions. RNNs are divided into 2 subgroups: Gradient Neural Networks (GNN) and Zhang Neural Networks (ZNN). The GNN design is explicit and mostly applicable to time-invariant problems, which means that coefficients of the equations which are addressed are constant matrices. ZNN models can be implicit and capable to solve time-varying problems where coefficients of the equations depend on the variable $t\in \mathbb{R},t>0,$ representing the time [31].

GNN neural design models for computing the inverse or the Moore-Penrose inverse and linear matrix equations were proposed in [6,25,26,27]. Further, various dynamical systems aimed to approximating the pseudo-inverse of rank-deficient matrices were originated in [6]. Wei in [28] proposed three RNN models for approximation of the weighted Moore-Penrose inverse. Cichocki in [29] proposed a feed-forward neural network for approximating the Drazin inverse. Online matrix inversion in complex matrix case was considered in [30].

Applications of this type of inverses can be found in important areas such as modeling of electrical circuits [10], estimation of DNA sequences [22,23], balancing chemical equations [11,12] and in other important research domains related to robotics [13] and statistics [14].

In the following sections we will focus on GNN and ZNN dynamic systems based on gradient of the objective function and their implementation. The main goal of this research is the analysis of convergence and the study of analytic solutions. In this study we are concerned with solving the matrix equation $AXB=D$ [20,21] in real time using the GNN model, denoted by GNN$(A,B,D)$ [3,4,5,6,7,8,9], and the novel gradient-based GGNN model, termed as GGNN$(A,B,D)$. The proposed GGNN model is defined evolving the standard GNN dynamics along the gradient of the standard error matrix. The convergence analysis reveals the global asymptotic convergence of GGNN$(A,B,D)$ without restrictions, while the output belongs to the set of general solutions to the matrix equation $AXB=D$. The implementation is performed in MATLAB Simulink and numerical experiments are developed with simulations of GNN and GGNN models.

The GNN used to solve the general linear matrix equation $AXB=D$ is defined over the error matrix $E\left(t\right)=D-AV\left(t\right)B$, where $t\in [0,+\infty )$ is the time and $V\left(t\right)$ is an unknown state-variable matrix that approximates the unknown matrix X in $AXB=D$. The goal function is $\epsilon \left(t\right)=\left|\right|D-{AV\left(t\right)B\left|\right|}_F^2/2$ and its gradient is equal to

$$\frac{\partial \epsilon \left(t\right)}{\partial V}=\nabla \epsilon =\frac{1}{2}\frac{\partial \left|\right|D-{AV\left(t\right)B\left|\right|}_F^2}{\partial V}=-A^T(D-AV\left(t\right)B)B^T.$$

The GNN evolutionary design is defined by the dynamic system

$$\dot{V}\left(t\right)=\frac{\mathrm{d}V\left(t\right)}{\mathrm{d}t}=-\gamma \frac{\partial \epsilon \left(t\right)}{\partial V},V\left(0\right)=V_0,$$

(1)

where $\gamma >0$ is real parameter used to speed up the convergence. Thus, the linear GNN aimed at solving $AXB=D$ is given by the following dynamics:

$$\dot{V}\left(t\right)=\gamma A^T(D-AV\left(t\right)B)B^T.$$

(2)

The dynamical flow (2) is denoted as as GNN$(A,B,D)$. The nonlinear GNN$(A,B,D)$ for solving $AXB=D$ is defined by

$$\dot{V}\left(t\right)=\gamma A^T\mathcal{F}(D-AV\left(t\right)B)B^T.$$

(3)

The function array $\mathcal{F}\left(C\right)$ is based on appropriate odd and monotonically increasing activation function, which is applicable to elements of a real matrix $C=\left(c_{ij}\right)\in {\mathbb{R}}^{m\times n}$, i.e. $\mathcal{F}\left(C\right)=\left(f\left(c_{ij}\right)\right),i=1,\dots ,m,j=1,\dots ,n,$.

Proposition 1 restates restrictions about the solvability of $AXB=D$ and its general solution.

Proposition 1.

[1,18] If $A\in {\mathbb{R}}^{m\times n},B\in {\mathbb{R}}^{p\times q}$ and $D\in {\mathbb{R}}^{m\times q}$ then fulfillment of the condition

$$A{A}^{†}D{B}^{†}B=D$$

(4)

is necessary and sufficient for solvability of the linear matrix equation $AXB=D$. In this case, the set of all solutions is given by

$$X=\left\{A^{†}D{B}^{†}+Y-A^{†}AYB{B}^{†}|Y\in {\mathbb{R}}^{n\times p}\right\}.$$

(5)

The following results from [15] describe the conditions of convergence and the limit of the unknown matrix $V\left(t\right)$ from (3) as $t\to +\infty$.

Proposition 2.

[15] Suppose the matrices $A\in {\mathbb{R}}^{m\times n},B\in {\mathbb{R}}^{p\times q}$ and $D\in {\mathbb{R}}^{m\times q}$ satisfy (4). Then the unknown matrix $V\left(t\right)$ from (3) converges as $t\to +\infty$ with the equilibrium state

$$V\left(t\right)\to \tilde{V}=A^{†}D{B}^{†}+V\left(0\right)-A^{†}AV\left(0\right)B{B}^{†}$$

(6)

for any initial state variable matrix $V\left(0\right)\in {\mathbb{R}}^{n\times p}.$

The research in [24] investigated various ZNN models based on optimization methods. The goal of current research is to develop the GNN model based on the gradient $E_G\left(t\right)$ of ${\parallel E\left(t\right)\parallel }_F^2$ instead of the original goal function $E\left(t\right)$.

Obtained results are summarized as follows.

• A novel error function $E_G\left(t\right)$ is proposed for development of the GNN dynamical evolution.
• GNN design evolved upon the error function $E_G\left(t\right)$ is developed and analyzed theoretically and numerically.
• A hybridization of GNN and ZNN dynamical systems based on the error matrix $E_G$ is proposed and investigated.

Global organization of section is as follows. Motivation and derivation of GGNN and GZNN models are presented in Section 2. Section 3 is aimed to convergence analysis of the GGNN dynamics. Numerical comparison on GNN and GGNN dynamics are given in Section 4, in which SubSection 4.1 investigates a practical application of GGNN to electrical networks. Neural dynamics based on the hybridization of GGNN and GZNN model for solving matrix equations are considered in Section 6. Numerical examples on hybrid models model are analysed in Section 6. Finally, the last section presents some concluding remarks and vision of further research.

2. Motivation and derivation of the GGNN and GZNN models

The standard GNN design (2) solves the GLME $AXB=D$ under the constraint (4). Our goal is to resolve this restriction and propose dynamic evolutions based on the error functions that tend to zero without restrictions.

Our goal is to define the GNN design for solving the GLME $AXB=D$ based on the error function

$$E_G\left(t\right):=\nabla \epsilon \left(t\right)=A^T\left(D-AV\left(t\right)B\right)B^T=A^{T}E\left(t\right)B^T.$$

(7)

The equilibria points of (7) satisfy

$$E_G\left(t\right):=\nabla \epsilon \left(t\right)=0.$$

We continue investigation from [24]. More precisely, we will develop the GNN model based on the error function $E_G\left(t\right)$ instead of the error function $E\left(t\right)$. In this way, new neural dynamics is aimed to force the gradient $E_G$ to zero instead of the standard goal function $E\left(t\right)$. It is reasonable to call such RNN model as gradient-based GNN (GGNN shortly).

Proposition 3 gives conditions for solvability of the matrix equations $E\left(t\right)=0$ and $E_G\left(t\right)=0$ and general solutions to these systems.

Proposition 3.

[24] Consider arbitrary matrices $A\in {\mathbb{R}}^{m\times n}$, $B\in {\mathbb{R}}^{k\times h}$ and $D\in {\mathbb{R}}^{m\times h}$. The next statements are true.

(a) The equation $E\left(t\right)=0$ is solvable if and only if (4) is satisfied and the general solution to $E\left(t\right)=0$ is given by (5).

(b) The equation $E_G\left(t\right)=0$ is always solvable and its general solution coincides with (5).

In this way, the matrix equation $E\left(t\right)=0$ is solvable under the condition (4), while the equation $E_G\left(t\right)=0$ is always consistent. In addition, the general solutions to equations $E\left(t\right)=0$ and $E_G\left(t\right)=0$ are identical [24].

The next step is to define the GGNN dynamics using the error matrix $E_G\left(t\right)$. Let us define the objective function ${\epsilon }_G=\left|\right|E_G{\left|\right|}_F^2/2$, whose gradient is equal to

$$\frac{\partial {\epsilon }_G\left(V\left(t\right)\right)}{\partial V}=\frac{\partial \left|\right|A^T(D-AV\left(t\right)B)B^T{\left|\right|}_F^2}{\partial V}=-A^{T}A\left(A^T(D-AV\left(t\right)B)B^T\right)B{B}^T.$$

The dynamical system for the GGNN formula is obtained applying the GNN evolution along the gradient of ${\epsilon }_G\left(V\left(t\right)\right)$ based on $E_G\left(t\right)$, as follows

$$\begin{array}{cc}\hfill \dot{V}\left(t\right)& =-\gamma \frac{\partial {\epsilon }_G}{\partial V}\hfill \\ \hfill & =\gamma A^{T}A\left(A^T\left(D-AV\left(t\right)B\right)B^T\right)B{B}^T.\hfill \end{array}$$

(8)

The nonlinear GGNN dynamics is defined as

$$\dot{V}\left(t\right)=\gamma A^{T}A\mathcal{F}\left(A^T\left(D-AV\left(t\right)B\right)B^T\right)B{B}^T,$$

(9)

in which $\mathcal{F}\left(C\right)$ denotes an odd and monotonically increasing function array, as mentioned in previous section for the GNN model (3). An arbitrary monotonically increasing odd activation function $f(·)$ is used for the construction of the GNN neural design. The model (9) is termed as GGNN$(A,B,D)$. Figure 1 represents the Simulink implementation of GGNN$(A,B,D)$ dynamics (9).

On the other hand, the GZNN model defined upon the Zhangian matrix $E_G\left(t\right)$ is defined in [24] by the general evolution design

$${\dot{E}}_G\left(t\right)=\frac{\mathrm{d}{E}_G\left(t\right)}{\mathrm{d}t}=-\gamma \mathcal{F}\left(E_G\left(t\right)\right).$$

(10)

3. Convergence analysis of GGNN dynamics

In this section, we will analyze convergence properties of GGNN model given by dynamics (9).

Theorem 1.

Consider matrices $A\in {\mathbb{R}}^{m\times n},B\in {\mathbb{R}}^{p\times q}$ and $D\in {\mathbb{R}}^{m\times q}$. If an odd and monotonically increasing array activation function $\mathcal{F}(·)$ based on an elementwise function $f(·)$ is used, then the neural state matrix $V\left(t\right)\in {\mathbb{R}}^{n\times p}$ of the $\mathrm{GGNN}(A,B,D)$ model (9) asymptotically converges to the solution of the matrix equation $AXB=D$, i.e., $A^{T}AV\left(t\right)B{B}^T\to A^{T}D{B}^T\mathrm{as}t\to +\infty$, for an arbitrary initial state matrix $V\left(0\right)$.

Proof. 

From the statement b) of Proposition 3, the solvability of $A^{T}AVB{B}^T=A^{T}D{B}^T$ is ensured. The substitution $V\left(t\right)=\overline{V}\left(t\right)+A^{†}D{B}^{†}$ transforms the dynamics (9) into

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\overline{V}\left(t\right)}{\mathrm{d}t}& =\frac{\mathrm{d}V\left(t\right)}{\mathrm{d}t}=\gamma A^{T}A\mathcal{F}\left(A^T\left(D-AV\left(t\right)B\right)B^T\right)B{B}^T\hfill \\ \hfill & =\gamma A^{T}A\mathcal{F}\left(A^T\left(D-A\overline{V}\left(t\right)B-A{A}^{†}D{B}^{†}B\right)B^T\right)B{B}^T\hfill \\ \hfill & \stackrel{\left(4\right)}{=}\gamma A^{T}A\mathcal{F}\left(A^T\left(D-A\overline{V}\left(t\right)B-D\right)B^T\right)B{B}^T\hfill \\ \hfill & =-\gamma A^{T}A\mathcal{F}\left(A^{T}A\overline{V}\left(t\right)B{B}^T\right)B{B}^T.\hfill \end{array}$$

(11)

Lyapunov function candidate which measures the convergence performance is defined by

$$L\left(\overline{V}\left(t\right),t\right)=\frac{1}{2}\left|\right|\overline{V}\left(t\right){\left|\right|}_F^2=\frac{1}{2}\mathrm{Tr}\left(\overline{V}{\left(t\right)}^T\overline{V}\left(t\right)\right).$$

(12)

The conclusion is $L(\overline{V}\left(t\right),t)\ge 0$. According to (12), assuming (11) and using $\mathrm{d}\mathrm{Tr}\left(X^{T}X\right)=2\mathrm{Tr}\left(X^T\mathrm{d}X\right)$ in conjunction with basic properties of the matrix trace function, one can express the time derivative of $L(\overline{V}\left(t\right),t)$ as in the following

$$\begin{array}{cc}\hfill \frac{\mathrm{d}L(\overline{V}\left(t\right),t)}{\mathrm{d}t}& =\frac{1}{2}\frac{\mathrm{d}\mathrm{Tr}\left(\overline{V}{\left(t\right)}^T\overline{V}\left(t\right)\right)}{\mathrm{d}t}\hfill \\ \hfill & =\frac{1}{2}·2·\mathrm{Tr}\left(\overline{V}{\left(t\right)}^T\frac{\mathrm{d}\overline{V}\left(t\right)}{\mathrm{d}t}\right)\hfill \\ \hfill & =\mathrm{Tr}\left[\overline{V}{\left(t\right)}^T\left(-\gamma A^{T}A\mathcal{F}\left(A^{T}A\overline{V}\left(t\right)B{B}^T\right)B{B}^T\right)\right]\hfill \\ \hfill & =-\gamma \mathrm{Tr}\left[\overline{V}{\left(t\right)}^T{A}^{T}A\mathcal{F}\left(A^{T}A\overline{V}\left(t\right)B{B}^T\right)B{B}^T\right]\hfill \\ \hfill & =-\gamma \mathrm{Tr}\left[B{B}^T\overline{V}{\left(t\right)}^T{A}^{T}A\mathcal{F}\left(A^{T}A\overline{V}\left(t\right)B{B}^T\right)\right]\hfill \\ \hfill & =-\gamma \mathrm{Tr}\left[{\left(A^{T}A\overline{V}\left(t\right)B{B}^T\right)}^T\mathcal{F}\left(A^{T}A\overline{V}\left(t\right)B{B}^T\right)\right].\hfill \end{array}$$

(13)

Since the scalar-valued function $f(·)$ is an odd and monotonically increasing, it follow for $W\left(t\right)=A^{T}A\overline{V}\left(t\right)B{B}^T$

$$\begin{array}{cc}\hfill \frac{\mathrm{d}L(\overline{V}\left(t\right),t)}{\mathrm{d}t}& =-\gamma \mathrm{Tr}\left[\left(W^T\mathcal{F}\left(W\right)\right)\right]\hfill \\ \hfill & =-\gamma \sum _{i=1}^m\sum _{j=1}^n{w}_{ij}f\left(w_{ij}\right)\left\{\begin{array}{ccc}\hfill <0& \mathrm{if}& W\left(t\right):=A^{T}A\overline{V}\left(t\right)B{B}^T\ne 0\hfill \\ \hfill =0& \mathrm{if}& W\left(t\right):=A^{T}A\overline{V}\left(t\right)B{B}^T=0,\hfill \end{array}\right.\hfill \end{array}$$

(14)

which implies

$$\frac{\mathrm{d}L(\overline{V}\left(t\right),t)}{\mathrm{d}t}\left\{\begin{array}{ccc}\hfill <0& \mathrm{if}& W\left(t\right)\ne 0\hfill \\ \hfill =0& \mathrm{if}& W\left(t\right)=0.\hfill \end{array}\right.$$

(15)

Observing the identity

$$\begin{array}{cc}\hfill W\left(t\right)& =A^{T}A\overline{V}\left(t\right)B{B}^T\hfill \\ \hfill & =A^{T}A\left(V\left(t\right)-A^{†}D{B}^{†}\right)B{B}^T\hfill \\ \hfill & =A^{T}AV\left(t\right)B{B}^T-A^{T}D{B}^T\hfill \\ \hfill & =A^T\left(AV\left(t\right)B-D\right)B^T,\hfill \end{array}$$

and using the Lyapunov stability theory, $W\left(t\right):=A^T\left(AV\left(t\right)B-D\right)B^T$ globally converges to the zero matrix, from arbitrary initial value $V\left(0\right)$. □

Theorem 2.

The activation state variables matrix $V\left(t\right)$ of the model $\mathrm{GGNN}(A,B,D)$, defined by (9), is convergent as $t\to +\infty$ and its equilibrium state is

$$V\left(t\right)\to \tilde{V}\left(t\right)=A^{†}D{B}^{†}+V\left(0\right)-A^{†}AV\left(0\right)B{B}^{†}$$

(16)

for every initial state matrix $V\left(0\right)\in {\mathbb{R}}^{n\times p}$.

Proof. 

From (9), the matrix $V_1\left(t\right)={\left(A^{T}A\right)}^{†}{A}^{T}AV\left(t\right)B{B}^T{\left(B{B}^T\right)}^{†}$ satisfies

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{V}_1\left(t\right)}{\mathrm{d}t}& ={\left(A^{T}A\right)}^{†}{A}^{T}A\frac{\mathrm{d}V\left(t\right)}{\mathrm{d}t}B{B}^T{\left(B{B}^T\right)}^{†}\hfill \\ \hfill & =\gamma {\left(A^{T}A\right)}^{†}{A}^{T}A\left[A^{T}A\left(A^T(D-AV\left(t\right)B)B^T\right)B{B}^T\right]B{B}^T{\left(B{B}^T\right)}^{†}.\hfill \end{array}$$

According to the basic properties of the Moore–Penrose inverse [18,19], it follows

$${\left(B{B}^T\right)}^{T}B{B}^T{\left(B{B}^T\right)}^{†}={\left(B{B}^T\right)}^T=B{B}^T,{\left(A^{T}A\right)}^{†}{A}^{T}A{\left(A^{T}A\right)}^T={\left(A^{T}A\right)}^T=A^{T}A$$

which further implies

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{V}_1\left(t\right)}{\mathrm{d}t}& =\gamma A^{T}A\left(A^T(D-AV\left(t\right)B)B^T\right)B{B}^T\hfill \\ & =\frac{\mathrm{d}V\left(t\right)}{\mathrm{d}t}.\hfill \end{array}$$

Consequently, $V_2\left(t\right)=V\left(t\right)-V_1\left(t\right)$ satisfies $\frac{\mathrm{d}{V}_2\left(t\right)}{\mathrm{d}t}=\frac{\mathrm{d}V\left(t\right)}{\mathrm{d}t}-\frac{\mathrm{d}{V}_1\left(t\right)}{\mathrm{d}t}=0$, which implies

$$\begin{array}{cc}\hfill V_2\left(t\right)& =V_2\left(0\right)\hfill \\ \hfill & =V\left(0\right)-V_1\left(0\right)\hfill \\ \hfill & =V\left(0\right)-{\left(A^{T}A\right)}^{†}{A}^{T}AV\left(0\right)B{B}^T{\left(B{B}^T\right)}^{†}\hfill \\ \hfill & =V\left(0\right)-A^{†}AV\left(0\right)B{B}^{†},t\ge 0.\hfill \end{array}$$

(17)

Furthermore, from Theorem 1, $A^{T}AV\left(t\right)B{B}^T\to A^{T}D{B}^T$ and $V_1\left(t\right)$ converges to

$$\begin{array}{ll}\hfill V_1\left(t\right)& ={\left(A^{T}A\right)}^{†}{A}^{T}AV\left(t\right)B{B}^T{\left(B{B}^T\right)}^{†}\to {\left(A^{T}A\right)}^{†}{A}^{T}D{B}^T{\left(B{B}^T\right)}^{†}\\ & =A^{†}D{B}^{†}\end{array}$$

as $t\to +\infty$. Therefore, $V\left(t\right)=V_1\left(t\right)+V_2\left(t\right)$ converges to the equilibrium state

$$\tilde{V}\left(t\right)=A^{†}D{B}^{†}+V_2\left(t\right)=A^{†}D{B}^{†}+V\left(0\right)-A^{†}AV\left(0\right)B{B}^{†}.$$

The proof is finished. □

4. Numerical experiments on GNN and GGNN dynamics

Numerical examples in this section are represented based on the Simulink implementation of GGNN formula in Figure 1. Three activation functions $f(·)$ are used in numerical experiments:

1. linear function

$$f_{lin}\left(x\right)=x;$$

(18)

2. power-sigmoid activation function

$$f_{ps}(x,\rho ,\varrho )=\left\{\begin{array}{ccc}\hfill x^{\rho }& \mathrm{if}& \left|x\right|\ge 1\hfill \\ \hfill \frac{1+e^{-\varrho }}{1-e^{-\varrho }}·\frac{1+e^{-\varrho x}}{1-e^{-\varrho x}}& \mathrm{if}& \left|x\right|<1\hfill \end{array}\right.$$

(19)

where $\varrho >2$ and $\rho \ge 3$ is odd integer.

3. smooth power-sigmoid function

$$f_{sps}(x,\rho ,\varrho )=\frac{1}{2}{x}^{\rho }+\frac{1+e^{-\varrho }}{1-e^{-\varrho }}·\frac{1+e^{-\varrho x}}{1-e^{-\varrho x}},$$

(20)

where $\varrho >2$ and $\rho \ge 3$ is odd integer.

The parameter $\gamma$, initial state $V\left(0\right)$ and parameters $\rho$ and $\varrho$ of the nonlinear activation functions (19) and (20), are entered directly in the model, while matrices A, B and D are defined from the workspace. It is assumed $\rho =\varrho =3$ in all examples. The ode15s differential equation solver is used in configuration parameters.

The blocks powersig, smoothpowersig and transpmult include the codes described in [15].

Example 4.1.

Let us consider the idempotent matrix A from [32,33]

$$A=\left[\begin{array}{cccc}1& 0& 1& 1\\ 0& 1& 1& 2\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]$$

with the theoretical Moore-Penrose inverse

$$V^{*}=A^{†}=\frac{1}{3}\left[\begin{array}{cccc}2& -1& 0& 0\\ -1& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right]$$

under the input parameters $\gamma ={10}^8$, $B=D=I_4$, $V\left(0\right)=O_4$, where $I_4$ and $O_4$ denote the identity and zero $4\times 4$ matrix, respectively. The Simulink implementation from Figure 1 exports the graphical results in Figure 2 and Figure 3 which display the behavior of $\left|\right|A^T(D-AV\left(t\right)B)B^T{\left|\right|}_F$ and $\left|\right|V\left(t\right)-V^{*}{\left|\right|}_F$ respectively. It is observable that the norms generated by the application of GGNN formula vanish faster to zero against corresponding norms in the GNN model. Graphs in presented figures strengthen the fast convergence of the GGNN dynamical system and the important role which can include the application of this specific model (9) to problems that require the computation of the Moore-Penrose inverse.

Example 4.2.

Let us consider the matrices

$$A=\left[\begin{array}{ccc}-8& 8& -4\\ 11& 4& -7\\ 1& -4& 3\\ 0& 12& -10\\ 6& 12& -12\end{array}\right],B=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right],D=\left[\begin{array}{ccc}-84& 2524& 304\\ -2252& -623& 2897\\ 484& -885& -701\\ -1894& 2278& 2652\\ -2778& 1524& 3750\end{array}\right].$$

Ranks of input matrices are equal to $r=\mathrm{rank}\left(A\right)=2$, $\mathrm{rank}\left(D\right)=2$ and $\mathrm{rank}\left(B\right)=3$. The linear GGNN formula GGNN$(A,B,D)$ (9) is applied to solve the matrix equation $AXB=D$, which gives in the case $V\left(0\right)=0$

$$X=A^{†}D{B}^{†}=\left[\begin{array}{cccc}-113.9846& -147.1385& 137.7385& 0\\ -74.4615& 136.1538& 107.8462& 0\\ 100.0462& -64.4154& -135.7846& 0\end{array}\right].$$

The gain parameter of the model is $\gamma ={10}^9$, $V\left(0\right)=0$ and the final time is $t=0.00001$.

Elementwise trajectories of the variable state matrix $V\left(t\right)$ with red lines are shown in Figure 4a–c, for linear, power-sigmoid and smooth power-sigmoid activation functions, respectively. It is observable the convergence of elementwise trajectories to the black dashed lines of the theoretical solution X. The trajectories in figures indicate a usual convergence behaviour, so the system is globally asymptotically stable. The norm of the error matrix $E_G$ of both model for linear and non-linear activation function are shown on Figure 5a–c. The nonlinear activation function shows superiority in the convergence speed, comparing with the linear activation function. On each graph, Frobenius norm of the error from GGNN formula vanish faster to zero than GNN model which strengthens the fact that the proposed dynamical system (9) includes accelerated convergence property against (3).

Example 4.3.

Let us explore behavior of GGNN$(A{A}^T,A^{T}A,A)$ for computing the Moore-Penrose inverse. We consider matrix:

$$A=\left[\begin{array}{ccc}9& 3& -3\\ -1& 1& 0\\ 4& 7& 2\\ 2& 4& -4\\ 13& 5& 8\end{array}\right].$$

The error matrix $E\left(t\right)=A^T(I-AV)$ initiates the GNN$(A^{T}A,I,A^T)$ dynamics for computing $A^{†}$. Corresponding error matrix for GGNN$({\left(A^{T}A\right)}^2,I,A^{T}A{A}^T)$ is

$$\begin{array}{cc}\hfill E_G\left(t\right)& =A^{T}A\left(A^T-A^{T}AV\right)=A^{T}A{A}^T\left(I-AV\right).\hfill \end{array}$$

Rank of the input matrix is equal to $r=\mathrm{rank}\left(A\right)=3$. The gain parameter of the model is $\gamma =100$, initial state is $V\left(0\right)=0$ with stop time $t=0.00001$. The Moore–Penrose inverse of A is given by

$$A^{†}\approx \left[\begin{array}{ccccc}0.0775& -0.0235& -0.0538& -0.0069& 0.0390\\ -0.0445& 0.0376& 0.1310& 0.0655& -0.0167\\ -0.0826& 0.0077& 0.0252& -0.0577& 0.0589\end{array}\right]$$

The error matrix $E\left(t\right)=D-AV\left(t\right)B$ of GNN and GGNN model for both linear and non-linear activation functions are shown on Figure 6a–c, and $E_G$ of both model for linear and non-linear activation function are shown on Figure 7a–c. Like in previously example, we can conclude that GGNN converges faster compared to the GNN model.

Example 4.4.

Consider the matrices

$$A=\left[\begin{array}{ccccc}15& -352& -45& -238& 42\\ -5& 14& 8& 132& -65\\ 235& -65& 44& 350& -73\end{array}\right],D=\left[\begin{array}{ccc}-4& 4& 16\\ 3& 1& -9\\ 1& -7& 2\\ 2& 2& -4\\ 4& 1& -5\end{array}\right],A_1=DA,$$

that dissatisfy $\mathrm{rank}\left(A_1\right)=\mathrm{rank}\left(D\right)=3$. Now, we apply GGNN formula to solve the matrix equation $A_{1}X=D$. Notice, if we look at GGNN$(A,B,D)$ in this example, we actually operate with $A=A_1$ and $B=I_3$. So, we consider GNN$(A_1,I_3,D)$. The error matrix for the corresponding GGNN model is

$$\begin{array}{cc}\hfill E_G\left(t\right)& =A_1^T\left(D-A_{1}V{I}_3\right)I_3^T=A_1^T\left(D-A_{1}V\right).\hfill \end{array}$$

The gain parameter of the model is $\gamma ={10}^9$ and the final time is $t=0.00001$. The zero initial state $V\left(0\right)=0$ generates the best approximate solution $X=A_1^{†}D={\left(DA\right)}^{†}D$ of the matrix equation $A_{1}X=D$, given by

$$X=A_1^{†}D=\left[\begin{array}{ccc}-0.000742151758732906& -0.00913940509710217& 0.00337604725603635\\ -0.00281859443979587& -0.00383668701022646& -0.000268348265455338\\ -0.000377755170367094& -0.00134463893719598& 0.000458599151069315\\ -0.000177073414005706& 0.00403153859974667& 0.000743238210989296\\ -0.000956081467310272& -0.00748631706136413& 0.00124829423173208\end{array}\right].$$

The Frobenius norm of the error matrix $D-AV\left(t\right)B$ in GNN and GGNN models for both linear and non-linear activation functions are shown on Figure 8a–c, and error matrix $E_G$ in both models for linear and non-linear activation function are shown on Figure 9a–c. It is observable that GGNN converges faster than GNN.

Example 4.5.

Table 2 shows the results obtained during experiments we conducted with nonsquared matrices, where $m\times n$ is dimension of matrix. Table 1 is the input data in the Simulink model which are used to perform experiments that generate results in Table 2. Best cases in Table 2 are marked in bold text.

Numerical results arranged in Table 2 are divided into two parts by a horizontal line. The upper part corresponds to test matrices of dimensions $\le 10$ while the lower part corresponds to the dimensions $m,n\ge 10$. Considering first two columns, it is observable from the upper part that GGNN generates smaller values ${\left|\right|E\left(t\right)\left|\right|}_F$ compared to GGNN. Values ${\left|\right|E\left(t\right)\left|\right|}_F$ in the lower part generated by GNN and GGNN are equal. Considering the third and fourth columns, it is observable from the upper part that GGNN generates smaller values $\left|\right|E_G\left(t\right){\left|\right|}_F$ compared to GGNN. On the other hand, values $\left|\right|E_G\left(t\right){\left|\right|}_F$ in the lower part, generated by GGNN, are smaller than corresponding values generated by GNN. Last two columns show that the GGNN requires smaller CPU time compared to GNN. General conclusion is that GGNN model is more efficient on rank deficient test matrices of larger order $m,n\ge 10$.

4.1. Application of GGNN to electrical networks

It is really interesting to apply the novel GGNN formula (9) for the calculation and study of different parameters related to electrical networks. For this reason the following circuit in Figure 10 from [34] is considered.

Our goal is to estimate the currents ${\mathrm{I}}_1,{\mathrm{I}}_2,{\mathrm{I}}_3$ in amperes $\left(\mathrm{A}\right)$ while electrical potential $\mathrm{E}$ is measured in volts $\left(\mathrm{V}\right)$ and resistors with resistance $\mathrm{R}$ is measured in ohms $\left(\Omega \right)$. Applying the current law in the points $\mathrm{A},\mathrm{B}$, relationships ${\mathrm{I}}_1={\mathrm{I}}_2+{\mathrm{I}}_3$ and $8{\mathrm{I}}_1+4{\mathrm{I}}_2=20$ obtained respectively and from the voltage law the relationship $20{\mathrm{I}}_3-4{\mathrm{I}}_2=16$ which lead to the following system $A\mathrm{I}=D$ in matrix form

$$\left[\begin{array}{ccc}1& -1& -1\\ 2& 1& 0\\ 0& 1& 5\end{array}\right].\left[\begin{array}{c}{\mathrm{I}}_1\\ {\mathrm{I}}_2\\ {\mathrm{I}}_3\end{array}\right]=\left[\begin{array}{c}0\\ 5\\ -4\end{array}\right].$$

with theoretical solution ${\mathrm{I}}_1=2,{\mathrm{I}}_2={\mathrm{I}}_3=1$. For the parameters $\gamma ={10}^6$ and the zero initial condition simulink implementation from Figure 1 extracts Figure 11 and Figure 12.

It is observable from Figure 11 that GGNN$(A,1,D)$ initiates a faster convergence than GNN$(A,1,D)$ formula for the same parameter $\gamma$ as error $\left|\right|V\left(t\right)-{\mathrm{I}\left|\right|}_F$ vanishes faster to zero. Figure 12 presents state trajectories of ${\mathrm{I}}_1,{\mathrm{I}}_2,{\mathrm{I}}_3$ obtained by the exact solution with the trajectories resulting from the GGNN formula for $\gamma ={10}^6$ and $V\left(0\right)=0$. These observations indicate that the proposed GGNN formula for solving general linear matrix equations is usable in solving electrical networks, which is an interesting engineering problem.

5. Mixed GGNN-GZNN model for solving matrix equations

Let us define gradient error matrix of the matrix equation $AX=B$ by

$$E_{A,I,B}\left(t\right)=A^T\left(AV\left(t\right)-B\right).$$

The GZNN design (10) corresponding to the error matrix $E_{A,I,B}$, marked with GZNN$(A,I,B)$, is of the form:

$${\dot{E}}_{A,I,B}\left(t\right)=-{\gamma }_1\mathcal{F}\left(A^T\left(AV\left(t\right)-B\right)\right).$$

(21)

Now, the scalar-valued norm-based error function corresponding to $E_{A,I,B}\left(t\right)$ is given by

$$\epsilon \left(t\right)=\epsilon \left(V\left(t\right)\right)=\frac{1}{2}\left|\right|E_{A,I,B}\left(t\right){\left|\right|}_F=\frac{\left|\right|A^T\left(AV\left(t\right)-B\right){\left|\right|}_F}{2}.$$

The following dynamic state equation can be derived using the GGNN$(A,I,B)$ design formula derived from (9):

$$\dot{V}\left(t\right)=-{\gamma }_2{A}^{T}A\mathcal{F}\left(A^T\left(AV\left(t\right)-B\right)\right).$$

(22)

Further, it follows

$${\dot{E}}_{A,I,B}\left(t\right)=A^{T}A\dot{V}\left(t\right)=-{\gamma }_2{A}^{T}A{A}^{T}A\mathcal{F}\left(A^T\left(AV\left(t\right)-B\right)\right).$$

(23)

Next step is to define new hybrid model based on the summation of (21) and (23) in the case ${\gamma }_1={\gamma }_2=\frac{1}{2}\gamma$, as follows:

$$\begin{array}{cc}\hfill {\dot{E}}_{A,I,B}\left(t\right)& =A^{T}A\dot{V}\left(t\right)\hfill \\ \hfill & =-\gamma \left({\left(A^{T}A\right)}^2+I\right)\mathcal{F}\left(A^T\left(AV\left(t\right)-B\right)\right)\hfill \\ \hfill & =-\gamma \left({\left(A^{T}A\right)}^2+I\right)\mathcal{F}\left(E_{A,I,B}\left(t\right)\right).\hfill \end{array}$$

(24)

The model (24) is derived as a combination of the model GGNN$(A,I,B)$ and the model GZNN$(A,I,B)$. Hence, it is equally justified to use the term Hybrid GGNN (HGGNN shortly) or Hybrid GZNN (HGZNN shortly) model. But, the model (24) is implicit, so that it is not a kind of GGNN dynamics. On the other hand, it is designed for time-invariant matrices, which is not in accordance with the common nature of GZNN models, because usually GZNN is used for the time-varying case. A formal comparison of (24) and GZNN$(A,I,B)$ reveals that both these methods possess identical left hand sides and the right hand side of (24) can be derived multiplying the right hand side of GZNN$(A,I,B)$ by the term ${\left(A^{T}A\right)}^2+I$.

Formally, (24) is closer to the GZNN dynamics, so, we will denote the model (2.4) by HGZNN$(A,I,B)$, considering that this model is not the exact GZNN neural dynamics and it is applicable in time-invariant case. This is the case of constant coefficient matrices A, I, B. Figure 13 represents the Simulink implementation of HGZNN$(A,I,B)$ dynamics (24).

Now, we will take into account solving the matrix equation $XC=D$. The error matrix for this equation is defined by

$$E_{I,C,D}\left(t\right)=\left(V\left(t\right)C-D\right)C^T.$$

The GZNN design (10) corresponding to the error matrix $E_{I,C,D}$, marked with GZNN$(I,C,D)$, is of the form:

$${\dot{E}}_{I,C,D}\left(t\right)=\dot{V}C{C}^T=-{\gamma }_1\mathcal{F}\left(\left(V\left(t\right)C-D\right)C^T\right).$$

(25)

On the other hand, the GGNN design formula (9) produces the following dynamic state equation:

$$\dot{V}\left(t\right)=-{\gamma }_2\mathcal{F}\left((V\left(t\right)C-D)C^T\right)C{C}^T,V\left(0\right)=V_0.$$

(26)

The GGNN model (26) is denoted by GGNN$(I,C,D)$. It implies

$${\dot{E}}_{I,C,D}\left(t\right)=\dot{V}\left(t\right)C{C}^T=-{\gamma }_2\mathcal{F}\left((V\left(t\right)C-D)C^T\right)C{C}^{T}C{C}^T.$$

(27)

A new hybrid model based on the summation of (25) and (27) in the case $\gamma =2{\gamma }_1=2{\gamma }_2$ can be proposed as follows

$$\begin{array}{cc}\hfill {\dot{E}}_{I,C,D}\left(t\right)& =\dot{V}\left(t\right)C{C}^T\hfill \\ \hfill & =-\gamma \mathcal{F}\left((V\left(t\right)C-D)C^T\right)\left(I+{\left(C{C}^T\right)}^2\right)\hfill \\ \hfill & =-\gamma \mathcal{F}\left(E_{I,C,D}\left(t\right)\right)\left(I+{\left(C{C}^T\right)}^2\right).\hfill \end{array}$$

(28)

The model (28) will be denoted by HGZNN$(I,C,D)$. This is the case of constant coefficient matrices I, C, D.

For the purposes of the proof of the following results, we will denote by $ECR\left(\mathcal{M}\right)$ the exponential convergence rate of the model $\mathcal{M}$. With ${\lambda }_{min}\left(K\right)$ and ${\lambda }_{max}\left(K\right)$, we denote the smallest and largest eigenvalue of a matrix K, respectively. In the continuation of the work we use three types of activation functions $\mathcal{F}$: linear, power-sigmoid and smooth power-sigmoid.

Next theorem determines the equilibrium state of HGZNN$(A,I,B)$ and defines its global exponential convergence.

Theorem 3.

Let $A\in {\mathbb{R}}^{k\times n},B\in {\mathbb{R}}^{k\times m}$ be given and satisfy $A{A}^{†}B=B$ and $V\left(t\right)\in {\mathbb{R}}^{n\times m}$ be the state matrix of (24), where $\mathcal{F}$ is defined by $f_{lin}$, $f_{ps}$ or $f_{sps}$.

a) Then $V\left(t\right)$ achieves global convergence and satisfies $AV\left(t\right)\to B$ when $t\to +\infty$, starting from any initial state $X\left(0\right)\in {\mathbb{R}}^{n\times m}$. The state matrix $V\left(t\right)\in {\mathbb{R}}^{n\times m}$ of HGZNN $(A,I,B)$ is stable in the sense of Lyapunov.

b) The exponential convergence rate of the  HGZNN $(A,I,B)$ model (24) in the linear case is equal to

$$ECR\left(\mathrm{HGZNN}(A,I,B)\right)=\gamma \left(1+{\lambda }_{min}\left({\left(A^{T}A\right)}^2\right)\right).$$

(29)

c) The activation state variables matrix $V\left(t\right)$ of the model HGZNN $(A,I,B)$ is convergent when $t\to +\infty$ with the equilibrium state matrix

$$V\left(t\right)\to {\tilde{V}}_{V\left(0\right)}=A^{†}B+(I-A^{†}A)V\left(0\right).$$

(30)

Proof. 

a) With the assumption $A{A}^{†}B=B$ we have solvability of the matrix equation $AX=B$.

We can define the Lyapunov function as

$$\begin{array}{cc}\hfill \mathcal{L}\left(t\right)& =\frac{1}{2}\left|\right|E_{A,I,B}\left(t\right){\left|\right|}_F^2\hfill \\ \hfill & =\frac{1}{2}\mathrm{Tr}\left({\left(E_{A,I,B}\left(t\right)\right)}^T{E}_{A,I,B}\left(t\right)\right).\hfill \end{array}$$

Hence, from (24) and $\mathrm{d}\mathrm{Tr}\left(V^{T}V\right)=2\mathrm{Tr}\left(V^T\mathrm{d}V\right)$, it holds that

$$\begin{array}{cc}\hfill \dot{\mathcal{L}}\left(t\right)& =\frac{1}{2}\mathrm{d}\mathrm{Tr}\left({\left(E_{A,I,B}\left(t\right)\right)}^T{E}_{A,I,B}\left(t\right)\right)\hfill \\ \hfill & =\frac{1}{2}2\mathrm{Tr}\left({\left(E_{A,I,B}\left(t\right)\right)}^T\mathrm{d}{E}_{A,I,B}\left(t\right)\right)\hfill \\ \hfill & =\mathrm{Tr}\left({\left(E_{A,I,B}\left(t\right)\right)}^T{\dot{E}}_{A,I,B}\left(t\right)\right)\hfill \\ \hfill & =\mathrm{Tr}\left({\left(E_{A,I,B}\left(t\right)\right)}^T\left(-\gamma \left({\left(A^{T}A\right)}^2+I\right)\mathcal{F}\left(E_{A,I,B}\left(t\right)\right)\right)\right)\hfill \\ \hfill & =-\gamma \mathrm{Tr}\left(\left({\left(A^{T}A\right)}^2+I\right)\mathcal{F}\left(E_{A,I,B}\left(t\right)\right){\left(E_{A,I,B}\left(t\right)\right)}^T\right).\hfill \end{array}$$

In the linear case it follows

$$\dot{\mathcal{L}}\left(t\right)=-\gamma \mathrm{Tr}\left(\left({\left(A^{T}A\right)}^2+I\right)E_{A,I,B}\left(t\right){\left(E_{A,I,B}\left(t\right)\right)}^T\right).$$

We also consider next inequality [35], which is valid for a real symmetric matrix K and a real symmetric positive-semidefinite matrix L of the same size:

$${\lambda }_{min}\left(K\right)Tr\left(L\right)\le \mathrm{Tr}\left(KL\right)\le {\lambda }_{max}\left(K\right)\mathrm{Tr}\left(L\right).$$

Now, it can be chosen: $K={\left(A^{T}A\right)}^2+I$ and $L=E_{A,I,B}\left(t\right){\left(E_{A,I,B}\left(t\right)\right)}^T$. Let ${\lambda }_{min}\left({\left(A^{T}A\right)}^2\right)\ge 0$ be the minimal eigenvalue of ${\left(A^{T}A\right)}^2$. Then $1+{\lambda }_{min}\left({\left(A^{T}A\right)}^2\right)\ge 1$ is the minimal nonzero eigenvalue of ${\left(A^{T}A\right)}^2+I$, which implies

$$\dot{\mathcal{L}}\left(t\right)\le -\gamma \left(1+{\lambda }_{min}\left({\left(A^{T}A\right)}^2\right)\right)\mathrm{Tr}\left(E_{A,I,B}\left(t\right){\left(E_{A,I,B}\left(t\right)\right)}^T\right).$$

(31)

From (31), it can be concluded

$$\dot{\mathcal{L}}\left(t\right)=\left\{\begin{array}{cc}\hfill <0& \mathrm{if}{E}_{A,I,B}\left(t\right)\ne 0\\ \hfill =0& \mathrm{if}{E}_{A,I,B}\left(t\right)=0.\end{array}\right.$$

(32)

According to (32), the Lyapunov stability theory confirms that $E_{A,I,B}\left(t\right)=AV\left(t\right)-B=0$ is a globally asymptotically stable equilibrium point of the HGZNN$(A,I,B)$ model (24). So, $E_{A,I,B}\left(t\right)$ converges to the zero matrix, i.e. $AV\left(t\right)\to B$ from any initial state $X\left(0\right)$.

b) From a) it follows that

$$\begin{array}{cc}\hfill \dot{\mathcal{L}}& \le -\gamma \left(1+{\lambda }_{min}\left({\left(A^{T}A\right)}^2\right)\right)\mathrm{Tr}\left({\left(E_{A,I,B}\left(t\right)\right)}^T{E}_{A,I,B}\left(t\right)\right)\hfill \\ \hfill & =-\gamma \left(1+{\lambda }_{min}\left({\left(A^{T}A\right)}^2\right)\right)\left|\right|E_{A,I,B}\left(t\right){\left|\right|}_F^2\hfill \\ \hfill & =-2\gamma \left(1+{\lambda }_{min}\left({\left(A^{T}A\right)}^2\right)\right)\mathcal{L}\left(t\right).\hfill \end{array}$$

This implies

$$\begin{array}{cc}\hfill \mathcal{L}& \le \mathcal{L}\left(0\right)e^{-2\gamma \left(1+{\lambda }_{min}\left({\left(A^{T}A\right)}^2\right)\right)t}⟺\hfill \\ \hfill \left|\right|E_{A,I,B}\left(t\right){\left|\right|}_F^2& \le \left|\right|E_{A,I,B}\left(0\right){\left|\right|}_F^2{\mathrm{e}}^{-2\gamma \left(1+{\lambda }_{min}\left({\left(A^{T}A\right)}^2\right)\right)}⟺\hfill \\ \hfill \left|\right|E_{A,I,B}\left(t\right){\left|\right|}_F& \le \left|\right|E_{A,I,B}\left(0\right){\left|\right|}_F{\mathrm{e}}^{-\gamma \left(1+{\lambda }_{min}\left({\left(A^{T}A\right)}^2\right)\right)},\hfill \end{array}$$

which confirms the convergence rate (29) of HGZNN$(A,I,B)$.

c) This part of the proof can be verified by following an analogous result from [17]. □

Theorem 4.

Let $C\in {\mathbb{R}}^{m\times l},D\in {\mathbb{R}}^{n\times l}$ be given and satisfy $D{C}^{†}C=D$ and $V\left(t\right)\in {\mathbb{R}}^{n\times m}$ be the state matrix of (28), where $\mathcal{F}$ is defined by $f_{lin}$, $f_{ps}$ or $f_{sps}$.

a) Then $V\left(t\right)$ achieves global convergence $V\left(t\right)C\to D$ when $t\to +\infty$, starting from any initial state $V\left(0\right)\in {\mathbb{R}}^{n\times m}$. The state matrix $V\left(t\right)\in {\mathbb{R}}^{n\times m}$ of HGZNN $(I,C,D)$ is stable in the sense of Lyapunov.

b) The exponential convergence rate of the  HGZNN $(I,C,D)$ model (28) in the linear case is equal to

$$ECR\left(\mathrm{HGZNN}(I,C,D)\right)=\gamma \left(1+{\lambda }_{min}\left({\left(C{C}^T\right)}^2\right)\right).$$

(33)

c) The activation state variables matrix $V\left(t\right)$ of the model HGZNN $(I,C,D)$ is convergent when $t\to +\infty$ with the equilibrium state matrix

$$V\left(t\right)\to {\tilde{V}}_{V\left(0\right)}=D{C}^{†}+V\left(0\right)(I-C{C}^{†}).$$

(34)

Proof. 

a) With the assumption $D{C}^{†}C=D$ we have solvability of the matrix equation $XC=D$.

Lets define the Lyapunov function with

$$\begin{array}{cc}\hfill \mathcal{L}\left(t\right)& =\frac{1}{2}\left|\right|E_{I,C,D}\left(t\right){\left|\right|}_F^2\hfill \\ \hfill & =\frac{1}{2}\mathrm{Tr}\left({\left(E_{I,C,D}\left(t\right)\right)}^T{E}_{I,C,D}\left(t\right)\right).\hfill \end{array}$$

Hence, from (28) and $\mathrm{d}\mathrm{Tr}\left(X^{T}X\right)=2\mathrm{Tr}\left(X^T\mathrm{d}X\right)$, it holds that

$$\begin{array}{cc}\hfill \dot{\mathcal{L}}\left(t\right)& =\frac{1}{2}\mathrm{d}\mathrm{Tr}\left({\left(E_{I,C,D}\left(t\right)\right)}^T{E}_{I,C,D}\left(t\right)\right)\hfill \\ \hfill & =\mathrm{Tr}\left({\left(E_{I,C,D}\left(t\right)\right)}^T{\dot{E}}_{I,C,D}\left(t\right)\right)\hfill \\ \hfill & =\mathrm{Tr}\left({\left(E_{I,C,D}\left(t\right)\right)}^T\left(-\gamma \left({\left(C{C}^T\right)}^2+I\right)\mathcal{F}\left(E_{I,C,D}\left(t\right)\right)\right)\right)\hfill \\ \hfill & =-\gamma \mathrm{Tr}\left(\left({\left(C{C}^T\right)}^2+I\right)\mathcal{F}\left(E_{I,C,D}\left(t\right)\right){\left(E_{I,C,D}\left(t\right)\right)}^T\right).\hfill \end{array}$$

According to similar results from [36], one can verify the following inequality

$$\dot{\mathcal{L}}\left(t\right)=-\gamma \mathrm{Tr}\left(\left({\left(C{C}^T\right)}^2+I\right)E_{I,C,D}\left(t\right){\left(E_{I,C,D}\left(t\right)\right)}^T\right).$$

We also consider next inequality [35], which is valid for a real symmetric matrix K and a real symmetric positive-semidefinite matrix L of the same size:

$${\lambda }_{min}\left(K\right)\mathrm{Tr}\left(L\right)\le \mathrm{Tr}\left(KL\right)\le {\lambda }_{max}\left(K\right)\mathrm{Tr}\left(L\right).$$

Now, it can be chosen: $K={\left(C{C}^T\right)}^2+I$ and $L=E_{I,C,D}\left(t\right){\left(E_{I,C,D}\left(t\right)\right)}^T$.

Let ${\lambda }_{min}\left({\left(C{C}^T\right)}^2\right)\ge 0$ be the minimal eigenvalue of ${\left(C{C}^T\right)}^2$. Then $1+{\lambda }_{min}\left({\left(C{C}^T\right)}^2\right)\ge 1$ is the minimal nonzero eigenvalue of ${\left(C{C}^T\right)}^2+I$. This implies

$$\dot{\mathcal{L}}\left(t\right)=-\gamma \left(1+{\lambda }_{min}\left({\left(C{C}^T\right)}^2\right)\right)\mathrm{Tr}\left(E_{I,C,D}\left(t\right){\left(E_{I,C,D}\left(t\right)\right)}^T\right).$$

(35)

From (35), it can be concluded

$$\dot{\mathcal{L}}\left(t\right)=\left\{\begin{array}{cc}\hfill <0& \mathrm{if}{E}_{I,C,D}\left(t\right)\ne 0\\ \hfill =0& \mathrm{if}{E}_{I,C,D}\left(t\right)=0.\end{array}\right.$$

(36)

According to (36), the Lyapunov stability theory confirms that $E_{I,C,D}\left(t\right)=V\left(t\right)C-D=0$ is a globally asymptotically stable equilibrium point of the HGZNN$(A,I,B)$ model (28). So, $E_{I,C,D}\left(t\right)$ converges to the zero matrix, i.e. $V\left(t\right)C\to D$ from any initial state $V\left(0\right)$.

b) From a) it follows

$$\begin{array}{cc}\hfill \dot{\mathcal{L}}& \le -\gamma \left(1+{\lambda }_{min}\left({\left(C{C}^T\right)}^2\right)\right)Tr\left({\left(E_{I,C,D}\left(t\right)\right)}^T{E}_{I,C,D}\left(t\right)\right)\hfill \\ \hfill & =-\gamma \left(1+{\lambda }_{min}\left({\left(C{C}^T\right)}^2\right)\right)\left|\right|E_{I,C,D}\left(t\right){\left|\right|}_F^2\hfill \\ \hfill & =-2\gamma \left(1+{\lambda }_{min}\left({\left(C{C}^T\right)}^2\right)\right)\mathcal{L}\left(t\right).\hfill \end{array}$$

This implies

$$\begin{array}{cc}\hfill \mathcal{L}& \le \mathcal{L}\left(0\right)e^{-2\gamma \left(1+{\lambda }_{min}\left({\left(C{C}^T\right)}^2\right)\right)t}⟺\hfill \\ \hfill \left|\right|E_{I,C,D}\left(t\right){\left|\right|}_F^2& \le \left|\right|E_{I,C,D}\left(0\right){\left|\right|}_F^2{e}^{-2\gamma \left(1+{\lambda }_{min}\left({\left(C{C}^T\right)}^2\right)\right)}⟺\hfill \\ \hfill \left|\right|E_{I,C,D}\left(t\right){\left|\right|}_F& \le \left|\right|E_{I,C,D}\left(0\right){\left|\right|}_F{e}^{-\gamma \left(1+{\lambda }_{min}\left({\left(C{C}^T\right)}^2\right)\right)},\hfill \end{array}$$

which confirms that the convergence rate of HGZNN$(I,C,D)$ is

$$ECR\left(\mathrm{HGZNN}(I,C,D)\right)=\gamma \left(1+{\lambda }_{min}\left({\left(C{C}^T\right)}^2\right)\right).$$

c) This part of the proof can be verified by following an analogous result from [17]. □

Corollary 5.1.

a) Let the matrices $A\in {\mathbb{R}}^{k\times n},B\in {\mathbb{R}}^{k\times m}$ be given and satisfy $A{A}^{†}B=B$ and $V\left(t\right)\in {\mathbb{R}}^{n\times m}$ be the state matrix of (24), with an arbitrary nonlinear activation $\mathcal{F}$. Then $ECR\left(\mathrm{GZNN}\right(A,I,B\left)\right)=\gamma$.

b) Let the matrices $C\in {\mathbb{R}}^{m\times l},D\in {\mathbb{R}}^{n\times l}$ be given and satisfy $D{C}^{†}C=D$ and $V\left(t\right)\in {\mathbb{R}}^{n\times m}$ be the state matrix of (28) with an arbitrary nonlinear activation $\mathcal{F}$. Then $ECR\left(\mathrm{GZNN}\right(I,C,D\left)\right)=\gamma$.

5.1. Regularized HGZNN model for solving matrix equations

From Theorem 3 and Corollary 5.1 (a), it follows

$$\frac{ECR\left(\mathrm{HGZNN}\right(A,I,B\left)\right)}{ECR\left(\mathrm{GZNN}\right(A,I,B\left)\right)}=1+{\lambda }_{min}\left({\left(A^{T}A\right)}^2\right)\ge 1.$$

(37)

Similarly, according to Theorem 4 and Corollary 5.1 (b), it can be concluded that

$$\frac{ECR\left(\mathrm{HGZNN}\right(I,C,D\left)\right)}{ECR\left(\mathrm{GZNN}\right(I,C,D\left)\right)}=1+{\lambda }_{min}\left({\left(C{C}^T\right)}^2\right)\ge 1.$$

(38)

Convergence of HGZNN$(A,I,B)$ (resp. HGZNN$(I,C,D)$) can be improved in the case ${\lambda }_{min}\left({\left(A^{T}A\right)}^2\right)>0$ (resp. ${\lambda }_{min}\left({\left(C{C}^T\right)}^2\right)>0$). There exist two possible situations when the acceleration terms $A^{T}A$ and $C{C}^T$ improve the convergence. The first case assumes invertibility of A (resp. C), and the second case assumes left invertibility of A (resp. right invertibility of C). Still, in some situations the matrices A and C could be rank deficient. Hence, in the case when A and C are square and singular, it is useful to use the invertible matrices $A+\lambda I$ and $C+\lambda I$, $\lambda >0$ instead of A and C, and consider the models HGZNN$(A+\lambda I,I,B)$ and HGZNN$(I,C+\lambda I,D)$. In below are presented the convergence results considering nonsingularity of $A+\lambda I$ and $C+\lambda I$.

Corollary 5.2.

Let $A\in {\mathbb{R}}^{n\times n}$, $B\in {\mathbb{R}}^{n\times m}$ be given and $V\left(t\right)\in {\mathbb{R}}^{n\times m}$ be the state matrix of (24), where $\mathcal{F}$ is defined by $f_{lin}$, $f_{ps}$ or $f_{sps}$. Let $\lambda >0$ be a selected real number. Then the following statements are valid:

a) The state matrix $V\left(t\right)\in {\mathbb{R}}_r^{n\times m}$ of the model HGZNN $(A+\lambda I,I,B)$ converges globally to

$${\tilde{V}}_{V\left(0\right)}={(A+\lambda I)}^{-1}B,$$

when $t\to +\infty$, starting from any initial state $X\left(0\right)\in {\mathbb{R}}^{n\times m}$ and the solution is stable in the sense of Lyapunov.

b) The minimal exponential convergence rate of HGZNN $(A+\lambda I,I,B)$ in the case $\mathcal{F}=I$ is equal to

$$ECR\left(\mathrm{HGZNN}(A+\lambda I,I,B)\right)=\gamma \left(1+{\lambda }_{min}\left({\left({\left(A+\lambda I\right)}^T\left(A+\lambda I\right)\right)}^2\right)\right).$$

c) Let ${\tilde{V}}_{V\left(0\right)}$ be the limiting value of $V\left(t\right)$ when $t\to +\infty$. Then

$$\underset{\lambda \to 0}{lim}{\tilde{V}}_{V\left(0\right)}=\underset{\lambda \to 0}{lim}{\left(A+\lambda I\right)}^{-1}B.$$

(39)

Proof. 

Since $A+\lambda I$ is invertible, it follows $V={\left(A+\lambda I\right)}^{-1}B$.

From (30) and invertibility of $A+\lambda I$ we can get validity of a). In this case, it follows

$$\begin{array}{cc}\hfill {\tilde{V}}_{V\left(0\right)}& ={(A+\lambda I)}^{-1}B+(I-{(A+\lambda I)}^{-1}(A+\lambda I))V\left(0\right)\hfill \\ \hfill & ={(A+\lambda I)}^{-1}B+(I-I)V\left(0\right)\hfill \\ \hfill & ={(A+\lambda I)}^{-1}B.\hfill \end{array}$$

Part b) is proved analogously as in Theorem 3. Last part c) follows from a). □

Corollary 5.3.

Let $C\in {\mathbb{R}}^{m\times m}$, $D\in {\mathbb{R}}^{n\times m}$ be given and $V\left(t\right)\in {\mathbb{R}}^{n\times m}$ be the state matrix of (28), where $\mathcal{F}=I,\mathcal{F}={\mathcal{F}}_{ps}$ or $\mathcal{F}={\mathcal{F}}_{sps}$. Let $\lambda >0$ be a selected real number. Then the following statements are valid:

a) The state matrix $V\left(t\right)\in {\mathbb{R}}_r^{n\times m}$ of HGZNN $(I,C+\lambda I,D)$ converges globally to

$${\tilde{V}}_{V\left(0\right)}=D{(C+\lambda I)}^{-1},$$

when $t\to +\infty$, starting from any initial state $X\left(0\right)\in {\mathbb{R}}^{n\times m}$ and the solution is stable in the sense of Lyapunov.

b) The minimal exponential convergence rate of HGZNN $(I,C+\lambda I,D)$ in the case $\mathcal{F}=I$ is equal to

$$ECR\left(\mathrm{HGZNN}(I,C+\lambda I,D)\right)=\gamma \left(1+{\lambda }_{min}\left({\left(\left(C+\lambda I\right){\left(C+\lambda I\right)}^T\right)}^2\right)\right).$$

c) Let ${\tilde{V}}_{V\left(0\right)}$ be the limiting value of $V\left(t\right)$ when $t\to +\infty$. Then

$$\underset{\lambda \to 0}{lim}{\tilde{V}}_{V\left(0\right)}=\underset{\lambda \to 0}{lim}D{\left(C+\lambda I\right)}^{-1}.$$

(40)

Proof. 

It can be proved analogously to Corollary 5.2. □

Remark 1.

The notation $A_1=A+\lambda I$ and $C_1=C+\lambda I$ will be used. Main observations about the convergence properties of HGZNN $(A,I,B)$ and HGZNN $(I,C,D)$ are highlighted as follows.

 $1.$ Hybrid neural dynamics HGZNN $(A,I,B)$ (resp. HGZNN $(I,C,D)$) converge faster than GZNN $(A,I,B)$ (resp. GZNN $(I,C,D)$). The accelerated convergence rate is equal to $1+{\lambda }_{min}\left({\left(A^{T}A\right)}^2\right)>1$ (resp. $1+{\lambda }_{min}\left({\left(C{C}^T\right)}^2\right)>1$).

 $2.$ Regularized hybrid dynamics HGZNN $(A_1,I,B)$ and HGZNN $(I,C_1,D)$ are applicable even in the case if A and C are singular matrices.

 $3.$ HGZNN $(A_1,I,B)$(resp. HGZNN $(I,C_1,D)$) always faster converge than GZNN $(A_1,I,B)$ (resp. GZNN $(I,C_1,D)$). The accelerated convergence rate is $1+{\lambda }_{min}\left({\left(A_1^T{A}_1\right)}^2\right)>1$ (resp. $1+{\lambda }_{min}\left({\left(C_1{C}_1^T\right)}^2\right)>1$).

6. Numerical examples on hybrid models model

In this section the numerical examples are represented based on the Simulink implementation of HGZNN formula. The previously mentioned three types of activation functions $f(·)$ in (18), (19) and (20) will be used in the following examples. The parameters $\gamma$, initial state $V\left(0\right)$ and parameters $\rho$ and $\varrho$ of the nonlinear activation functions (19) and (20), are entered directly in the model, while matrices A, B, C and D are defined from the workspace. We assume that $\rho =\varrho =3$ in all examples. The ordinary differential equation solver in configuration paremeters is the ode15s. The blocks powersig, smoothpowersig and transpmult include the codes described in [16].

We present numerical examples in which we compare Frobenius norms $\left|\right|E_G{\left|\right|}_F$ and $\left|\right|A^{-1}B-V\left(t\right){\left|\right|}_F$ which are generated by HGZNN, GZNN and GGNN.

Example 6.1.

Consider the matrix

$$A=\left[\begin{array}{ccccc}0.49& 0.276& 0.498& 0.751& 0.959\\ 0.446& 0.68& 0.96& 0.255& 0.547\\ 0.646& 0.655& 0.34& 0.506& 0.139\\ 0.71& 0.163& 0.585& 0.699& 0.149\\ 0.755& 0.119& 0.224& 0.891& 0.258\end{array}\right].$$

In this example we compare HGZNN$(A,I,I)$ model with GZNN$(A,I,I)$ and GGNN$(A,I,I)$ considering all three types of activation functions. The gain parameter of the model is $\gamma ={10}^6$, initial state $V\left(0\right)=0$ and the final time is $t=0.00001$.

Elementwise trajectories of the state variable with red lines are shown in Figure 14a–c, for linear, power-sigmoid and smooth power-sigmoid activation functions, respectively, and is observable the converge to the black dashed lines of the theoretical solution X. Trajectories indicate a usual convergence behaviour, so the system is globally asymptotically stable. Error matrix $E_G$ of HGZNN, GZNN and GGNN model for both linear and non-linear activation functions are shown on Figure 15a–c, and error matrix $A^{-1}B-V\left(t\right)$ of both model for linear and non-linear activation function are shown on Figure 16a–c. On each graph, Frobenius norm of error from HGZNN formula vanish faster to zero than GZNN and GGNN model.

Example 6.2.

Consider the matrices

$$\begin{array}{cc}\hfill A& =\left[\begin{array}{cccccc}0.0818& 0.0973& 0.0083& 0.0060& 0.0292& 0.0372\\ 0.0818& 0.0649& 0.0133& 0.0399& 0.0432& 0.0198\\ 0.0722& 0.0800& 0.0173& 0.0527& 0.0015& 0.0490\\ 0.0150& 0.0454& 0.0391& 0.0417& 0.0984& 0.0339\\ 0.0660& 0.0432& 0.0831& 0.0657& 0.0167& 0.0952\\ 0.0519& 0.0825& 0.0803& 0.0628& 0.0106& 0.0920\end{array}\right],\hfill \\ \hfill B& =\left[\begin{array}{ccccccc}0.1649& 0.1813& 0.0851& 0.1197& 0.0138& 0.1437& 0.1558\\ 0.1965& 0.1759& 0.0625& 0.0942& 0.0639& 0.1937& 0.0847\\ 0.1460& 0.1636& 0.0323& 0.1392& 0.1062& 0.1063& 0.0182\\ 0.0688& 0.0521& 0.0358& 0.1400& 0.1309& 0.0650& 0.0533\\ 0.1168& 0.1189& 0.0846& 0.1277& 0.0815& 0.0211& 0.0307\\ 0.0216& 0.0045& 0.0188& 0.0067& 0.1640& 0.1222& 0.0562\end{array}\right].\hfill \end{array}$$

In this example, we compare HGZNN$(A,I,B)$ model with GZNN$(A,I,B)$ and GGNN$(A,I,B)$ considering all three types of activation functions. The gain parameter of the model is $\gamma =1000$, initial state $V\left(0\right)=0$ and the final time is $t=0.01$.

Elementwise trajectories of the state variable with red lines are shown in Figure 17a–c, for linear, power-sigmoid and smooth power-sigmoid activation functions, respectively, and is observable the converge to the black dashed lines of the theoretical solution X. We can see that trajectories indicate a usual convergence behaviour, so the system is globally asymptotically stable. The error matrix $E_G$ of HGZNN, GZNN and GGNN model for both linear and non-linear activation functions are shown on Figure 18a–c, and the residual matrix $A^{-1}B-X\left(t\right)$ of both models for linear and non-linear activation function are shown on Figure 19a–c. On each graph, for both error cases, the Frobenius norm of error HGZNN formula is similar to the Frobenius norm of error of GZNN model, and they both converges faster to zero than GGNN model.

7. Conclusions

We show that the error function lying in the basement of GNN and ZNN dynamical evolutions can be defined using the gradient of the Frobenius norm of the traditional error function $E\left(t\right)$. The result of this intention is usage of an original error function $E_G\left(t\right)$ for the basis of GNN dynamics which results in the proposed GGNN model. The results related to the GNN model (called GNN$(A,B,D)$) for solving the general matrix equation $AXB=D$ are extended on GGNN model (called GGNN$(A,B,D)$) in both theoretical and computational directions. In theoretical sense, the convergence of defined GGNN model is considered. It is shown that the neural state matrix $V\left(t\right)$ of the GGNN$(A,B,D)$ model asymptotically converges to the solution of the matrix equation $AXB=D$, for an arbitrary initial state matrix $V\left(0\right)$, and coincides with the general solution of the linear matrix equation. A number of applications of GNN(A, B, D) are considered. All applications are globally convergent. Several particular appearances of the general matrix equation are observed and applied in computing various classes of generalized inverses. Illustrative numerical examples and simulation results are obtained using Matlab Simulink implementation and presented to demonstrate validity of the derived theoretical results. The influence of various nonlinear activations on the GNN models is considered in both the theoretical and the computational direction. In the presented examples it can be concluded that GGNN model is faster and has smaller error compared to the GNN model.

Further research can be oriented to definition of finite-time convergent GGNN or GZNN models as well as definition of noise-tolerant GGNN or GZNN design.