### Abstract

Neural network (NN) approaches have been widely applied for modeling and identification of nonlinear multiple-input multiple-output (MIMO) systems. This paper proposes a stochastic analysis of a class of these NN algorithms. The class of MIMO systems considered in this paper is composed of a set of single-input nonlinearities followed by a linear combiner. The NN model consists of a set of single-input memoryless NN blocks followed by a linear combiner. A gradient descent algorithm is used for the learning process. Here we give analytical expressions for the mean squared error (MSE), explore the stationary points of the algorithm, evaluate the misadjustment error due to weight fluctuations, and derive recursions for the mean weight transient behavior during the learning process. The paper shows that in the case of independent inputs, the adaptive linear combiner identifies the linear combining matrix of the MIMO system (to within a scaling diagonal matrix) and that each NN block identifies the corresponding unknown nonlinearity to within a scale factor. The paper also investigates the particular case of linear identification of the nonlinear MIMO system. It is shown in this case that, for independent inputs, the adaptive linear combiner identifies a scaled version of the unknown linear combining matrix. The paper is supported with computer simulations which confirm the theoretical results.

##### Keywords:

Nonlinear system identification; Neural networks; Gradient descent; Statistical analysis### Introduction

Neural network [1] approaches have been extensively used in the past few years for nonlinear MIMO system modeling, identification and control where they have shown very good performances compared to classical techniques [2-6].

If these NN approaches are to be used in real systems, it is important for the algorithm designer and the user to understand their learning behavior and performance capabilities. Several authors have analyzed NN algorithms during the last two decades which considerably helped the neural network community to better understand the mechanisms of neural networks [1,7-15]. For example, the authors in [13] have studied a simple structure consisting of two inputs and a single neuron. The authors in [8] studied a memoryless single-input single-output (SISO) system identification model for the single neuron case. In [9] the authors proposed a stochastic analysis of gradient adaptive identification of nonlinear Wiener systems composed of a linear filter followed with a Zero-memory nonlinearity. The model was composed of a linear adaptive filter followed by an adaptive parameterized version of the nonlinearity. This study has been later generalized [16] for the analysis of stochastic gradient tracking of time-varying polynomial Wiener systems. In [12] the author analyzed NN identification of nonlinear SISO Wiener systems with memory for the case where the adaptive nonlinearity is a memoryless NN with an arbitrary number of neurons. The case of a nonlinear SISO Wiener-Hammerstein system (i.e., an adaptive filter followed by an adaptive Zero-memory NN followed by an adaptive filter) has been analyzed in [11].

This paper deals with a typical class of nonlinear MIMO systems (Figure
1) which is composed of *M* inputs, *M* memoryless nonlinearities, a linear combiner, and *L* outputs. This corresponds, for example, to MIMO channels used in wireless terrestrial
communications
[17-22], satellite communications
[23,24], amplifier modeling
[25], control of nonlinear MIMO systems
[6], etc. Recently, a neural network approach has been proposed to adaptively identify
the overall input–output transfer function of this class of MIMO systems and to characterize
each component of the system (i.e., the memoryless nonlinearities and the linear combiner)
[4]. The proposed NN model is composed of a set of memoryless NN blocks followed by an
adaptive linear combiner. Each part of the adaptive system aims at identifying the
corresponding part in the unknown MIMO system. The algorithm has been successfully
applied to system modeling, channel tracking, and fault detection.

**Figure 1.** **Nonlinear MIMO system.**

The purpose of this paper is to provide a stochastic analysis of NN modeling of this class of MIMO systems. The paper provides a general methodology that may be used to solve other problems in stochastic NN learning analysis. The methodology consists of splitting the study into simple structures, before studying the complete structure. Here, as a first step we start by analyzing a simple linear adaptive MIMO scheme (consisting of an adaptive matrix) that identifies the nonlinear MIMO system (i.e., problem of linear identification of a nonlinear MIMO system). Then we analyze a nonlinear adaptive system in which the nonlinearities are assumed to be known and frozen during the learning process, only the linear combiner is made adaptive. Finally, the complete adaptive scheme is analyzed taking into account the insights given by the analysis of the simpler structures. In our analytical approach, we derive the general formulas and recursions, which we apply to a case study that we believe is illustrative to the reader.

The paper is organized as follows. The problem statement is given in Section 2. The study of the simple structures is detailed in Section 3. Section 4 presents the analysis for the complete structure. Simulation results and illustrations are given in Section 5. Finally, conclusions and future work are given in Section 6.

### Problem statement

#### Nonlinear MIMO system

The class of nonlinear MIMO systems discussed in this paper is presented in Figure
1. Each input *x*_{i} (*n*) (*i =* 1,…,*M*) is nonlinearly transformed by a memoryless nonlinearity *g*_{i} (.). The outputs of these nonlinearities are then linearly combined by an *L* × *M* matrix *H =* [*h*_{ji}] (assumed in this paper to be constant). Matrix *H* is defined by the unknown system to be identified. For example, in wireless MIMO
communication systems, *M* is the propagation matrix representing the channel between *M* transmitting antennas and *L* receiving antennas.

The *j*^{th} output of the MIMO system is expressed as:

where *N*_{j} is a white Gaussian noise with variance σ_{0}^{2}. Let

The system input–output relationship can be expressed in a matrix form as:

#### Neural Network identification structure and algorithm

The neural network (Figure
2) is composed of *M* blocks. Each block *k* has a scalar input *x*_{k} (*n*) (*k =* 1,…,*M*), N neurons and a scalar output. The output of the k^{th} block is expressed as:

**Figure 2.** **NN identification structure.**

Where *f* is the NN activation function. *a*_{ki}, *c*_{ki}, *b*_{ki} are, respectively, the input weight, bias term, and output weight of the *i*^{th} neuron in the *k*^{th} block. The output *NN*_{k} of the *k*^{th} block is connected to the *j*^{th} output of the system through weight *w*_{jk}. The system *j*^{th} output is then expressed as:

Weights *w*_{jk} will be represented by an LxM matrix: *W =* [*w*_{jk}]. Let

Equations (4) can then be expressed in a matrix form as:

For the learning process, the NN parameters are updated so that to minimize the sum of the squared errors between the unknown system outputs and the corresponding outputs of the model (Figure 3):

**Figure 3.** **Adaptive learning scheme.**

The gradient descent recursions for weight adaptation are:

where *μ* is a small positive constant and

#### Case study

After the derivation of the general formulas, it is important that we apply them to
special cases in order to get closed-form expressions of the different recursions
that can be illustrated to the reader. We have chosen here a case study that we think
is good to illustrate our results. In this case study, the inputs *x*_{i} (*n*) will be assumed uncorrelated Zero-mean Gaussian variables with variance
*erf* function. The unknown nonlinear transfer functions are taken from a family of nonlinear
functions of the form
*α*_{i} and *β*_{i} are positive constants. These nonlinear functions are reasonable models for amplitude
conversions of nonlinear high power amplifiers (HPA) used in digital communications
[12,25,26]. Note that other nonlinear functions may be considered, however, explicit closed-form
solutions of the different derivations may not be possible.

### Study of simplified structures: Linear adaptation

Before analyzing the full structure, we will analyze the following simplified schemes which will help us understand the complete structure:

1. The adaptive system is composed of an adaptive linear combiner *W* (Section 3.1).

2. The adaptive system is composed of *W* and scaled versions of the unknown nonlinearities (Section 3.2).

#### Linear adaptive system

This section studies the linear adaptive system that tries to model the nonlinear MIMO system (Figure 4):

**Figure 4.** **Linear adaptive system.**

#### Mean weight behavior and Wiener solution

Since matrix W is linear, it will not be able to identify the nonlinear blocks. However,
we will see that it is able to identify matrix *H* to within a diagonal scaling matrix if the inputs are Zero-mean and independent.

The gradient descent update of matrix *W* is expressed as:

Averaging both sides of (11) and using the standard LMS assumption of small *μ*[10], we obtain:

Where

By setting the updating gradient term to Zero, it can be shown that this equation has a single stationary point (Wiener solution [10]) which is expressed by:

Following Equation (12), the mean weights can be expressed as a function of the initial condition as:

If *μ* is sufficiently small, the first term converges to 0 and the second term converges
to

Hence, the mean weights converge to the Wiener solution:

It can be easily shown that the stability condition on *μ* is
[10]:

where λ_{max} is the largest eigenvalue of the covariance matrix *R*_{XX}.

Note that for Zero-mean independent inputs, *U* is a diagonal matrix:

In this case, the linear adaptation allows the identification of matrix *W* to within a scaling matrix, which depends on the nonlinearities and the input signals.
As expected, the scaling matrix reduces to the identity matrix if *g*_{k} (*x*_{k}) = *x*_{k}.

- Application to the case study

For the particular nonlinear functions given in the case study (see Section 2.3), it is easy to show:

The mean weight transient recursions are expressed as:

Matrix *U* reduces to the following diagonal matrix:

#### Transient MSE and Wiener MSE

The transient MSE is determined by:

where:

Using the independence of noise and weights at time *n*, we get:

The total MSE is therefore expressed as:

Wiener MSE:

The Wiener MSE,
*W*_{0} = *HU*. It can be easily shown that:

It is clear from this equation that if the unknown functions are linear, then the
Wiener MSE reduces to the noise power. The MSE is always larger than ζ_{0} because of the misadjustment error introduced by the weight fluctuations.

Now we can write the MSE as a function of the Wiener MSE:

Let the instantaneous deviation of the matrix weights with respect to the Wiener solution be denoted by:

We have:

This expression is similar to that of the well-known LMS algorithm [10], and can be evaluated as the sum of the minimum error and excess error (or misadjustment) as:

where

The misadjustment is expressed as:

At the convergence, we have:

Derivation of the misadjustment:

From Equation (11) it is easy to show that the weight fluctuations follow the recursion:

Taking the mean of this equation and applying the orthogonality principle between the input vector and the Wiener error, we get:

Thus, as expected, if *μ* is sufficiently small *E*(*
*

*V*

*n*)) converges to 0.

Similarly, for each vector *V*_{j} we can obtain the following recursion:

The evaluation of the covariance matrix of the weight fluctuations is obtained by
multiplying both sides of Equation (34) by

These expectations are derived in Appendix III, which yields:

Taking into account that E(V_{j}(∞)) = 0, *K*_{VjVj} can be obtained by solving the following equation:

This expression holds for any input signal. It can be simplified if

It is now easy to determine the total misadjustment:

Note that, as expected, in the case of linear functions Δ(∞) reduces to:

The additional terms are due to the nonlinearities and they should be calculated specifically for each nonlinearity.

- Application to the case study

The MSE is expressed as:

The Wiener MSE is expressed in this case as:

#### Adaptive W, the nonlinearities are frozen and known with scale factors

In this section, matrix *W* is adaptive, the nonlinearities are frozen and known with scale factors (Figure
5).

**Figure 5.** **Simplified scheme, the nonlinearities are frozen and assumed known with a scaling
factor.**

#### Mean weight behavior and stationary points

The gradient descent update of matrix W is expressed as:

where

Averaging both sides of (43) and using the standard LMS assumption of small *μ*, we obtain:

These recursions have a single stationary point (Wiener solution) which is:

Following Equation (44), the mean weight behavior can be expressed as function of the initial condition as:

Hence, if *μ* is sufficiently small, it can be shown that the mean weights converge to the Wiener
solution:

The stability condition on *μ* is:

Where *λ*_{max} is the largest eigenvalue of the covariance matrix Ω^{2}*R*_{g (X)g(X)}.

Thus, if each nonlinear function *g*_{k} (.) is known with a scaling factor *η*_{k}, then weights *h*_{jk} will be identified by *w*_{jk} (to the inverse of the scaling factor).

#### MSE

We have:

where:

Using the independence of noise and weights at time *n*, we obtain:

The MSE is therefore expressed as:

Wiener MSE:

The Wiener MSE can be easily expressed as:

Therefore the Wiener MSE is equal to the noise floor: There are no terms due to the
nonlinearities. This is expected since the nonlinearities are known with a scaling
matrix Ω (we have seen that the scaling matrix is canceled by *W*_{0} since *W*_{0} = HΩ^{-1}).

Let Z (*n*) = Ω *g*(*X*(*n*)), we can then express the MSE as a function of ζ_{0}, the weight fluctuation vector *V*(*n*) = *W*(*n*)–*W*_{0}, and Z(*n*):

Similarly to Equation (29), the misadjustment is expressed as:

The steady state MSE is then expressed as:

Derivation of the misadjustment:

It is easy to show that the weight fluctuations follow the recursion:

Taking the mean of this equation and applying the orthogonality principle between the input vector and the Wiener error, we obtain:

Thus, as expected, if *μ* is sufficiently small, E(V(*n*)) converges to 0.

For each vector *V*_{j} we have similar recursions:

The evaluation of the covariance matrix of the weight fluctuations is obtained by
multiplying both sides of Equation (57) by

In a similar way as in Appendix III, *K*_{VjVj} (∞) can be obtained by solving the following equation:

This expression can not be further simplified because *R*_{ZZ} is not necessarily of the form

Therefore, tr(*R*_{ZZ}*K*_{VjVj} (∞)) should be calculated for each nonlinearity and for each Ω.

It is interesting to study the case where nonlinearities are known with the same scaling
factor, i.e., Ω = *ηI*. In this case, and if the input vectors are independent and the outputs of the nonlinearities
are Zero-mean and of equal variance

As expected, the total misadjustment reduces to:

Here the value of the misadjustment is similar to that of linear identification of a linear system (LMS algorithm). This is expected since in this case there are no errors due to the approximation of the nonlinearities.

#### Case study

For the case study, it is easy to show that

### Study of the full structure

This section deals with the full structure (Figures
2,
3). All the NN and matrix *W* weights are updated.

#### Mean weight transient behavior

We take the following notations for the weights:

The update of matrix W is expressed as:

Averaging both sides of (63) and using the standard LMS assumption of small *μ*, we obtain:

Where

These matrices are time-dependent since they depend on the NN block weights which are updated through time.

Using the scalar notation we have:

With these notations we have:

For the NN block weights we have:

These equations hold for any nonlinearity. In the following, we will calculate them explicitly for the case study described in Section 2.3

- Application to the case study:

Since the inputs are independent and Zero-mean, we have *K*_{i} (*a*_{km},*b*_{km})=0, *k*≠*i* , and (see Appendix I)

In the other hand we have:

Inserting these expressions in equations (66)-(69), we obtain:

The explicit expressions of the different derivatives are detailed in Appendix II.

#### Stationary points

We obtain the stationary points by setting to 0 the expectations of the updating gradient terms in (64) and(4.5-7).

For W, we obtain:

For *c*_{ki} we obtain the equations:

For *a*_{ki} we obtain the equations:

For *b*_{ki} we obtain:

The above equations are nonlinear in the NN variables. They can be solved numerically, but they are very difficult to solve analytically.

Convergence of the algorithm to the stationary points:

It is always interesting to show whether an algorithm is capable of converging to
its stationary points or not. In our case it is difficult to establish this, since
the updating equations of the weights are nonlinear, except for *W*.

In the case where the NN weights are frozen we can establish the convergence condition
for *W*.

In this case we have:

The covariance matrices are fixed, since in this case the NN weights are frozen.

*E*(*W*(*n*)) can be expressed as a function of the initial condition as:

If *μ* is sufficiently small, the steady state solution to (81) is:

Hence, the mean weights converge to the stationary point *W*_{0}, and the stability condition on *μ* is:

Where *λ*_{max} is the largest eigenvalue of the correlation matrix *R*_{NN(X)NN(X)}.

Application to the case study:

For the case study it can be shown that *U* reduces to a diagonal matrix:

where:

This indicates that weights *w*_{jk} are scaled versions of the unknown weights *h*_{jk}, the scale factor *γ*_{k} is the same for all the weights connecting the *k*^{th} NN block to the outputs and it depends only on block *k* weights. If the error is sufficiently small, the *k*^{th} block NN will approximate the *k*^{th} nonlinearity to the inverse of the scale factor.

### MSE expression

The transient MSE is determined by:

where:

Which can be expressed as:

Let G_{il} = E(*g*_{i}(*x*_{i} (*n*))*gl*(*x*_{i} (*n*))). Using the notations of Section 4.1, we have:

Application to the case study:

It can be easily shown that:

The *1*^{st} term of
^{nd} term is the signal power of the *j*^{th} MIMO output, the *3*^{rd} term is the sum of the individual contributions of the neurons weighed by *W* and *H* weights, the 4^{th} term represents the sum of the coupling terms between neurons inside the same block
weighed by W. Note that since the inputs are Zero-mean and independent, there are
no coupling terms between neurons in different blocks (as in Eq. (89)).

The total MSE is then expressed as:

Case of frozen NN weights:

It is interesting to see the behavior of the MSE in the case where the NN weights are frozen.

In this case we have:

Here the minimum MSE depends on the noise floor and on the NN approximation error
of the nonlinearities. It is clear from this equation and from Section 3.2 that, if
the NN blocks ideally identify the nonlinearities (to within scale factors), then
*ζ*_{0} reduces to the noise floor.

The MSE can be written as a function of *ζ*_{0} as:

The steady state MSE is in this case:

The misadjustment can be derived similarly as in Sections 3.1 and 3.2. We obtain a
similar equation as (53), by replacing *R*_{ZZ} by *R*_{NN(X)NN(X)}. The equation can not be simplified further.

It is interesting to notice, however, that if the NN blocks perfectly identify the nonlinearities and if the conditions above equation (60) are fulfilled, then:

### Simulation examples

In this section we present some simulation results which are applied to the case study
described in Section 2.3. In these simulations, we have considered a 2 × 2 MIMO system
(i.e., *M = L = 2*). For the parameterized nonlinearities we have chosen *α*_{1}=*α*_{2}=1, *β*_{1}=1, *β*_{2}=2. Unless otherwise specified, the inputs are uncorrelated Zero-mean white Gaussian
processes with σ_{xi}*=* 1. In the simulations, the unknown combining matrix was fixed and was taken as
*H* can be seen as the propagation matrix between 2 transmitting antennas and 2 receiving
antennas.

#### Linear adaptation

For the linear adaptation case (Section 3.1), the adaptive system is composed of a
2x2 matrix W. For the noise we have taken σ_{0} = 0.001. The mean weight recursions and the MSE transient behaviors (Figures
6,
7) have been estimated over 20 Monte Carlo (MC) simulations and compared to the theoretical
derivations (Equations (19) and (41)). This chosen number of MC simulations shows
excellent fit between the Theory and MC estimations which confirms the validity of
the different assumptions made. A larger number of MC simulations allows a better
smoothing of the curves, but the conclusions remain the same.

**Figure 6.** **Smoothed MSE vs. iteration number - Linear adaptation case.**

**Figure 7.** **Mean weight transient behavior - Linear adaptation case.**

Matrix W converges to a scaled version of H:

#### MSE surface for the full NN algorithm

We move now to the study of the full NN algorithm. In this simulation, we have taken
N = 3 neurons in each of the two NN blocks. The learning rate was fixed to *μ* =0.0045. Figure
8 shows the MSE surface (i.e., Eq. (90), with no time dependence) as a function of
*w*_{11} and *w*_{12} (the other parameters were fixed). It is clear that the MSE is quadratic in *w*_{11} and *w*_{21}. It presents a single global minimum (as shown in Equations (76 and 84)).

**Figure 8.** **MSE surface, ***MSE ***(***w*_{11}**, ***w*_{21}**).**

Figure
9 (resp. Figure
10) shows the MSE surface as a function of *a*_{11} and *c*_{11} (resp. *a*_{11} and *a*_{12}) (the other parameters were fixed). It can be noted the flat areas (plateau regions)
around the minima of the MSE surface. This explains the slow evolution of the NN weights
when the algorithm gets close to its convergence point.

**Figure 9.** **MSE surface, ***MSE ***(***a*_{11}**, ***c*_{21}**).**

**Figure 10.** **MSE surface, ***MSE ***(***a*_{11}**, ***a*_{12}**).**

The MSE evolution during the learning process (Equation (90)) has been compared to 20 MC estimations (Figure 11a): The theory shows very good fit with the simulation results. In the Figure, we notice that the MSE presents several phases (each phase is controlled by a time constant) which end by a plateau phase where the MSE decreases very slowly. This is a typical behavior of the backpropagation algorithm [1] which is fundamentally different from that of the linear adaptation scheme (Figure 6). Here the MSE error is much smaller. This is expected, since here the additional MSE error due to the nonlinearities (Eq. 92) is highly reduced because our NN blocks have correctly identified the unknown nonlinearities (Figure 12, 13). Here we are in a situation close to that of Section 3.2 (Equations 51-52).

**Figure 11.** **a: MSE during the learning process (theory and simulation), σ**_{0}**=0.001 and ***μ***=0.0045. b**: MSE for different values of noise variance, the learning rate was *μ*=0.0045 (theoretical results). **c** : MSE for different values of learning rate *μ* (the noise variance was set to σ_{0}=0.001), theoretical results.

**Figure 12.** **Identification of ***g*_{1 }**(***x***), by the 1**^{st }**NN block.**

**Figure 13.** **Identification of ***g*_{2 }**(***x***) by the 2**^{nd }**NN block.**

Figure
11b shows the MSE evolution during the learning process for different values of the
noise variance σ_{0}. It can be seen that, as the noise variance decreases, the MSE decreases. However,
below a certain value of σ_{0} (here σ_{0}=0.0005), the MSE curves are almost identical. This is because in this case, the weight
misadjustment error (for the linear part) and the nonlinear approximation error (of
the nonlinear memoryless part) are much higher than the error caused by the presence
of noise (see Eqs. 92-93).

Figure
11c investigates the influence of the learning rate *μ*. It can be seen that as *μ* increases (up to *μ*=0.002), the algorithm is faster and the MSE is lower at the end of the simulation
time. However, for *μ*>0.002, as *μ* increases, the algorithm is faster at the beginning of the learning process, but
the MSE is higher at the end of the simulation time. This is due to the misadjustment
error which is higher for higher *μ* (see, e.g., Eq. 95).

#### Mean weight transient behavior for the full NN algorithm

Here we keep the system described in Section 5.2. The mean weight recursions for the
linear combiner *W* and the two NN blocks are shown in Figure
14,
15,
16 for both theory and MC estimations. The theoretical and estimated curves are indistinguishable.
This confirms the validity of the different assumptions made in Sections 4.1 and 4.2.

**Figure 14.** **Mean weight behavior: Matrix W (theory and simulation).**

**Figure 15.** **Mean weight behavior: First NN block (theory and simulation).**

**Figure 16.** **Mean weight behavior: Second NN block (theory and simulation vs. iterations).**

Notice that, in Figure
14, *W* weights have a fast evolution at the beginning of the learning process (with values
approaching *H*×*U*(*n*) where U is a diagonal matrix). They then evolve slowly till the end of the learning
process. The slow evolution is justified by the plateau regions presented by the MSE
surface. At the end of this simulation, matrix U was close to a diagonal matrix:

Figure
12,
13 show that functions *g*_{1} (*x*) and *g*_{2} (*x*) have been correctly identified by the corresponding NN blocks (the NN functions
are normalized by the scaling factors γ_{1}=1/1.2702, γ_{1}=1/1.0946, respectively).

#### Impact of correlated inputs

In the simulations below we study the impact of correlated inputs. The input signal
vector is chosen here as a 2D Gaussian process with covariance matrix of the form
*μ*=0.0075. We have run several simulations for different values of the cross correlation.
Note that (*ρ*=0) corresponds to independent inputs, and (*ρ*=1) corresponds to the same input (i.e., *x*_{1}=*x*_{2}). The values of matrix U and the MSE after *n*=210^{5} iterations are shown in Table
1. It can be seen from Table
1 that, in practice, matrix *U* remains very close to a diagonal matrix, even for high correlation between inputs.
This indicates that the system is capable of correctly identifying the nonlinearities
even when the inputs are highly correlated. The identification performances for the
cases (*ρ*=0.6) and ( *ρ*=0.99) are illustrated in Figures
17,
18 and
19,
20, respectively. As expected, the MSE increases as the correlation between inputs increases
(Table
1). When *ρ*=1 (i.e., the two inputs are the same) the system is capable of correctly identifying
the overall MIMO input–output transfer function. However, in this case, it is not
capable of separating the nonlinearities (as U is not diagonal). The reason is that
in this case, the system is seen by the learning algorithm as a 1x2 SIMO system which
has several equivalent structures. Figure
21 shows an example of two equivalent structures. Therefore, the adaptive system is
structurally not able to separate the nonlinearities. It is worth to note that for
the case (*ρ*=0.999), the inputs look like noisy versions of each other (i.e., this is equivalent
to a 1x2 system identification problem with noisy inputs). Thus, the MSE for the case
(*ρ*=0.999 ) is larger than the MSE for the case (*ρ*=1).

**Table 1.** **Effect of correlated inputs**

**Figure 17.** **Identification of ***g*_{1 }**(***x***), case ***ρ***=0.6.**

**Figure 18.** **Identification of ***g*_{2 }**(***x***), case ***ρ***=0.6.**

**Figure 19.** **Identification of ***g*_{1 }**(***x***), case ***ρ***=0.99.**

**Figure 20.** **Identification of ***g*_{2 }**(***x***), case ***ρ***=0.99.**

**Figure 21.** **Example of two equivalent structures when the inputs are the same.**

### Conclusion and future work

The paper provides a statistical analysis of NN modeling and identification of a class of nonlinear MIMO systems. The study investigates the MSE error, mean weight behavior, stationary points, misadjustment error, and stability conditions. The unknown system is composed of a set of single-input memoryless nonlinearities followed by a combining matrix. The NN model is composed of a set of single-input memoryless NN blocks followed by an adaptive linear combiner. The paper is supported with simulation results which show good agreement between the theoretical recursions and MC simulations. Future work will focus on 3 research directions. The first will explore the theoretical findings in order to express the effect of the number of neurons on the transient and steady state behavior of the algorithm. The second research axis will investigate the case where matrix H is time-varying and/or with memory (this may have applications, for example, in adaptive control of nonlinear dynamical MIMO systems). Finally, we will study the algorithm behavior and performance for specific inputs (such as space-time coded signals used in wireless communications and their impact on the system capacity).

### Appendix I

1) **Calculation of ***F*_{kk} Let *x*_{1} and *x*_{2} be two zero-mean Gaussian variables such that

Thus, using the un-correlation criteria between *x*_{1} and *x*_{2} for *ρ*=0 , we have:

where

Note that in the biasless case (i.e. all the bias terms are set to 0) this expression reduces to:

The integral is then simple to calculate:

In the other hand, since

2) Calculation of *K*

The inside integral can be eliminated by integrating by parts on variable *x*.The integral is then evaluated by combining the terms in the exponentials and completing
the squares. This yields:

Again, a Taylor series expansion can be used to simplify this expression.Note that in the biasless case we have:

### Appendix II

The derivatives needed to compute the different recursions are expressed as follows: