Nonlinear System Identification Using Type-2 Fuzzy Recurrent Wavelet Neural Network

In this paper, the integration of Type-2 fuzzy set theory and recurrent wavelet neural network(WNN) is proposed to allow managing of non-uniform uncertainties for identifying non-linear dynamic system. The proposed Type-2 fuzzy WNN is inherently a recurrent multilayered network which constructed based on a set of Type-2 fuzzy rules and recurrent connections in the second layer of the FWNN. Each rule comprises a wavelet function in the consequent part. The structure has both advantages of recurrent and wavelet neural network which expand the basic ability of fuzzy neural network to deal with temporal problems. Both antecedent and consequent parameters update rules are derived based on the gradient descent method. The structure is applied in the identification of dynamic plants which is commonly used in the literature. Simulation result from the identification of a second-order non-linear plant confirms the better performance and effectiveness of the proposed structure.

Introduction
The use of conventional modeling approaches in the study of nonlinear dynamical system identification suffers from some deficiencies, including lack of precise, strongly nonlinear and time-varying behavior, formal knowledge about the system and a high degree of uncertainty. Under such conditions, modelfree approaches such as artificial neural networks, fuzzy logic systems and fusion of them which are universal approximators, have been developed to compensate for the effects of nonlinearities and system uncertainties without significant prior knowledge of system dynamics [1,2].
Type-2 fuzzy logic systems as an extension of its type-1 have better ability to deal with and model uncertainties and measurement noise [3]. Different optimization methods are applied to estimate the parameters of neural network and fuzzy neural systems from numerical data. The methods can be categories into two types: derivative-based (such as gradient descent, least square) and derivative-free optimization methods(such as particle swarm optimization and genetic algorithm).
Back propagation error method is mainly used for training feed-forward neural networks. Because of some deficiencies of BP algorithm like slow convergence rate, being able to trap to local optima and dependent to the initial condition, some differ-1 ent methods based on BP algorithm are developed such as Momentum, adaptive learning rate. These modifications accelerate the convergence of network training, but they have not the ability to avoid areas of local optima [4].
Wavelet function as the activation function can enhance the advantages of neural networks for faster learning ability and wavelet decomposition for identification purposes. In literature, wavelet neural network (WNN) and its integration with fuzzy logic (FWNN) to determine an optimal definition of the premise and consequent part of fuzzy rules, are applied in the identification and control of nonlinear dynamical systems [5][6][7][8].
To enhance the approximation accuracy of the network and capture the dynamic behavior of the system, different kinds of feedback loops in the topology of the network are embedded. The combination of RNN and fuzzy logic (RFNN) are applied to various practical application [9][10][11], and different methods for their parameters tuning are discussed.
In this study, a type-2 fuzzy recurrent wavelet neural network(T2FRWNN) is presented which integrates FWNN and recurrent neural network to benefits the advantages of both. The proposed dynamic identification structure is equipped the networks to process temporal information and learn sequences with using local internal feedback loops by feeding the lower and upper of type-2 membership function back to itself. Besides, network response in handling the uncertainties and imprecision can be achieved by utilization of type-2 FNN and wavelet activation function. For the parameter adaptation of the T2FRWNN, gradient descent is used.
The main body of this paper is organized as follows. In Section 2, the structure of T2FRWNN is introduced. In Section 3, the parameters update rules based on gradient descent are derived. In Section 4, simulation results are shown and discussed.
Finally, concluding remarks are given in Section 5.

Network structure
The T2FRWNN structure benefits from recurrent type-2 membership function in the premise part and crisp wavelet numbers for the consequent part. In other words, the structure is the recurrent version of the so-called A2-C0 fuzzy system [12] which implements a recurrent wavelet fuzzy model. Each fuzzy if-then rule base with I input variables considered here is as follows: where K is the number of membership function and f r is the output of rth rule (r = {1, 2..., N }).Ã ik is the kth type-2 fuzzy membership function related to ith input variable. ρ is the weight coefficient between the input and the hidden layer. Consequent part of above rule involved a wavelet function of input variables. Wavelets are defined by a family of functions a and b (a > 0, b ϵ R) as dilation and translation respectively. In this manuscript, among several families of wavelets, Mexican Hat is considered as mother wavelet function. where where, z ri = xi−bri ari . The structure of T2FRWNN proposed in this paper includes seven layers. In the first layer, the number of nodes is equal to the number of input signals. These 2 nodes are used for distributing input signals. In the second layer, each node corresponds to one linguistic term. For each input signal entering into the system, type-2 membership functions are used. The type-2 MFs have uncertain standard deviation and fixed center and the membership degree µ ik (x i ) and µ ik (x i ) are calculated according to (4) and (5).
where c ik is center of type-2 MF, σ ik and σ ik are the upper and lower standard deviation of the kth type-2 MF of ith input. Moreover, ξ ik is the recurrent parameter defined as (6) which store the past information of the network. It is to be noted that the feed back weights of nodes in layer 2 are interval values.
where θ ik and θ ik are considered as feed back weights of the nodes in this layer. Nodes in the third layer represent one fuzzy rule and perform a fuzzy meet operation on inputs from layer 2 using an algebraic product operation to obtain upper and lower firing strength w r and w r which are calculated as follows: Layer 4 determines the normalized values of the lower and the upper firing strength corresponding to each node in layer 3: Layer 5 is consequent layer. Nodes in this layer compute the product of normalized firing strengthw r ,w r and wavelet function of input variables. Layer 6 consists of two summation blocks, one for summation of upper and the other for summation of lower outputs of the previous layer. Layer 7 calculates the output of the network using (10) where q is the design parameter which enables to adjust the lower or the upper portions, depending on the level of certainty of the system [13]. Consequently, the adaptation laws for the parameters and the proof of the stability of the learning process are given using a time-varying q.

Parameter Update Rules For Learning
In order to obtain the learning algorithm for the parameters gradient descent method with adaptive learning rate is applied. The use of adaptive learning rate give an assurance of convergence and speeds up the learning of the network. Moreover, a momentum is used to accelerate the learning process. The parameters which should be updated are the parameters of premise and consequent part of fuzzy rules, including c ik , σ ik , σ ik , θ ik and θ ik as premise part of kth type-2 MF of ith input and a r , b r , and ρ r as the consequent part of the rth rule. The error function E(k) of T2FRWNN is defined by: where K is the number of input signals of the network. The parameters mentioned above are adjusted as follows: and for the coefficient q: where α is the learning rate, λ is the momentum, i is the number of input signals of the network (input neurons), r is the number of rules (r = {1, ..., N }) and k is the number of MFs (k = {1, ..., m}). The values of derivatives in (12) can be calculated by the following formulas (15)- (17): the coefficient q is updated as follows: The derivatives in (12) are determined by the following formulas (20)-(24): Figure 1: where,

Simulation Example
In this section, one simulation studies of nonlinear system identifications are considered regarding the proposed model and learning algorithm to evaluate the performance of proposed T2FRWNN. Fig.1 depicts the structure of the identification system. The inputs of the model are delayed values from the nonlinear plant output and delayed input signals from the plant. Here, the prob-lem is to find such values of parameters of proposed T2FRWNN by using them in the system for all input values the difference between plant output y(t) and the T2FRWNN output y N (t) will be minimum. As an example, identification of a second-order non-linear plant that has been used in [14] is considered. The process is described by the following difference equation: where y(k−1), y(k−2), y(k−3) are one-, two-and threestep delayed outputs of the plant, respectively, u(k) and u(k − 1) are current and one-step delayed inputs of the plant. In order to have a comparison, the identification of the same plant by using different models are considered in [14] and [15]. For the simulation studies, the input signal in this manuscript has the following expression:

Conclusion
This paper has proposed an T2FRWNN structure for identification of dynamic plants. According to the design aspect of the proposed structure, the Type-2 Tagaki-Sugeno fuzzy logic is integrated with wavelet neural network which makes a fuzzy input space into different wavelet-based subspaces. Besides, the recurrent part of the structure gives the abilities to attract dynamics and store temporary information. In the training aspect, the parameter update rules of the structure are derived based on the gradient descent algorithm. In the simulation example of system identification, the proposed T2FRWNN has a better performance in comparison with other models, despite smaller number of parameters.


extfieldsofstudy:
Computer Science
Engineering

text:
Non-Coherent Sensor Fusion via Entropy Regularized Optimal Mass Transport

This work presents a method for information fusion in source localization applications. The method utilizes the concept of optimal mass transport in order to construct estimates of the spatial spectrum using a convex barycenter formulation. We introduce an entropy regularization term to the convex objective, which allows for low-complexity iterations of the solution algorithm and thus makes the proposed method applicable also to higher-dimensional problems. We illustrate the proposed method's inherent robustness to misalignment and miscalibration of the sensor arrays using numerical examples of localization in two dimensions.

INTRODUCTION
Source localization using an array of sensors is a well established and still highly active field of research in signal processing, with applications in areas such as radar, sonar, and audio processing [1][2][3][4]. The localization problem is often formulated as a spatial spectral estimation problem, where signal sources are identified with peaks of power in the spectral representation. In contrast to the closely related problem of direction of arrival (DoA) estimation, where one aims to find the directions of incoming planar waves from far field sources, the problem of localization also includes the estimation of the range (i.e., distance to the sources) and thus needs to utilize that the incoming wavefronts are non-planar. Classical methods, originally used mainly for DoA estimation, such as beamforming techniques [5] and subspace methods like MU-SIC [6] and ESPRIT [7], may be readily extended to the localization problem, with more recent efforts including works exploiting spectral sparsity [8][9][10]. Many of these methods rely on the assumption of well-calibrated arrays, i.e., known sensor positions, in order to produce reliable estimates. However, This work was supported in part by the Swedish Research Council, Carl Trygger's foundation, and the Royal Physiographic Society in Lund. this assumption may in some applications be violated due to, e.g., calibration errors, perturbations of the sensor localizations, or synchronization errors. Therefore, considerable effort has been directed towards formulating estimators robust to such calibration errors (see, e.g., [11][12][13][14][15]). In some of the approaches to address this problems, techniques have been developed that divide the available sensors into subarrays, each assumed to be perfectly calibrated, although relative errors between the different subarrays are allowed for [16]. The DoA estimation and/or localization is then performed by minimizing measures of discrepancy between the measurements, i.e., the estimated joint covariance matrix of the complete set of sensors, and the assumed array geometry [15].
In this work, we instead consider localization by fusing information from separate sensor arrays, where we only assume to have access to measurements in the form of covariance matrices for the individual arrays, so-called noncoherent processing (see, e.g., [17][18][19][20]), i.e., we assume no information about cross-correlation between different arrays. Such a scenario may arise in distributed sensor networks, wherein sensor nodes transmit information to a central processing unit that performs the estimation [21,22]. In such settings, it is common to have restrictions on the communication bandwidth, which requires transmitting only aggregated data [21], e.g., in the form of covariance matrices. Also, as the separate sensor arrays may not be jointly synchronized, cross-correlations cannot be estimated from the raw data.
Extending on the work presented in [23], which proposes a robust method for DoA sensor fusion using optimal mass transport (OMT), we here generalize this approach to the source localization problem. One of the main ideas is to use OMT as a notion of distance between spectra (see, e.g., [24][25][26]), and in particular, the proposed method performs sensor fusion by computing the corresponding barycenter solutions based on incomplete measurements. In addition, to make the scheme computationally efficient, we introduce an entropy regularization term in the problem formulation and show how to solve the resulting optimization problem iteratively, similar to the popular Sinkhorn iteration scheme [27] (see also [28]). The robustness of the proposed method is demonstrated in numerical examples.

Near field localization
Consider an array consisting of p sensors located at y k ∈ R d , for k = 1, 2, . . . , p, and assume that point-like sources, located in a space S ⊂ R d , emit spherical waves impinging on the sensors. Then, the spatial distribution of the energy of the source signals may be representated by a non-negative function, or measure, Φ on S, referred to as the power spectrum. We will denote the space of such functions by M + (S). Letting λ denote the wavelength of the impinging waves, the covariance matrix of the sensor measurements is R = Γ(Φ), where the linear operator Γ : M(S) → C p×p is given by (see, e.g., [29]) where a(x) ∈ C p is the array manifold vector [30], i.e., with dm(x) denoting the Lebesgue measure. In other words, a positive semi-definite matrix R is a valid covariance for a sensor array if and only if R = Γ(Φ), for some Φ ∈ M + (S). From this, it may be noted that performing source localization based on sensor measurements corresponds to inferring a power spectrum Φ ∈ M + (S) from observations of the array covariance matrix R.
We will here consider a scenario wherein several sensor arrays measure the superposition of the impinging waves, i.e., when covariance matrices R j , for j = 1, . . . , J, are available, corresponding to J separate sensor arrays. For the case when all arrays are calibrated, i.e., when the position of the sensors are known perfectly, source localization thus corresponds to finding a single power spectrum Φ that is consistent with the observed matrices, i.e., Γ j (Φ) = R j , for j = 1, . . . , J, where Γ j is the linear operator corresponding to array j. However, for the case of miscalibration, the (assumed) operators Γ j contain errors, which may result in there being no single Φ ∈ M + (S) consistent with all covariance matrices. We will here adress this problem by utilizing OMT in order to induce robustness to the sensor fusion problem.

Optimal mass transport
Let Φ 0 and Φ 1 be two non-negative mass distributions on the space S. A transport plan M ∈ M + (S 2 ) for such a pair is a non-negative distribution on S 2 , with S 2 = S × S, where M (x 0 , x 1 ) describes the amount of mass transported from location x 0 ∈ S to location x 1 ∈ S. By defining a cost function c : S 2 → R describing the cost of moving mass, we may define the total cost associated with a transport plan as The OMT problem is then to find the feasible transport plan from Φ 0 to Φ 1 incurring minimal transportation cost (see, e.g., [31]), i.e., to solve where the constraints ensure that the transport plan is valid, i.e., matches the marginals Φ 0 and Φ 1 . The suitability of utilizing OMT costs as measures of distance has previously been considered for imposing metric structure on the space of power spectra [24], as well as to induce distances and smooth interpolants for covariance matrices [25]. In this work, the optimal transport cost T (Φ 0 , Φ 1 ) will be used for interpolating between mass distributions Φ 0 and Φ 1 . It is worth noting that the OMT cost quantifies location errors based on the distances on the underlying space, in this case S ⊂ R d , thereby allowing for offsets due to, e.g., calibration errors. This is in contrast to standard metrics (e.g., L 1 , L 2 ) which only compare functions point by point without taking the geometry of the underlying space, i.e., relative location of the points, into account.

PROPOSED METHOD
Consider the localization scenario with the source signals impinging on a set of J arrays, with array manifold vectors a j (x), generating the covariance matrices R j , j = 1, . . . , J.
Given the array geometries and covariance matrix estimates of each individual array, we aim to estimate the spatial spectrum Φ. That is, we assume that each sensor array is calibrated within itself, but there might be global calibration errors resulting in inaccurate knowledge of each sensor array's position and orientation. For the far field case, i.e., DoA estimation, an OMT formulation for fusing information from different sensor arrays has been proposed in [23]. There, the estimate was formed as the spectrum closest, in the OMT sense, to a set of J spectra, each consistent with a corresponding covariance matrix. Correspondingly, in the localization setting the spatial spectrum Φ may be found as the solution to the convex barycenter problem where Γ j denotes the linear operator (1) of array j. In order to find practical solutions to the barycenter problem in (5), the space S may be discretized as to be represented by n points x k ∈ S. Accordingly, the marginal spatial spectra are represented by the vectors Φ j ∈ R n + , and the cost function and transport plans are represented by matrices C = [c k, ] ∈ R n×n and M = [m k, ] ∈ R n×n + , respectively, where c k, = c(x k , x ), and m k, denotes the mass transported from x k to x . Although such a discretization yields finite-dimensional linear programs, finding an optimal transport plan can be computationally demanding due to the large number of variables.
For discretizations of the standard OMT problem (4), this was in [27] addressed by regularizing the problem by introducing an entropy term to the objective function, yielding the strictly convex program where D(M ) = n k, =1 (m k, log(m k, ) − m k, + 1) is an entropy term, > 0 a small constant, and 1 ∈ R n denotes a vector of ones. In [27], it was shown that this problem may be solved in a computationally efficient way, using so-called Sinkhorn iterations.

Efficient sensor fusion
In order to arrive at a computationally efficient scheme for computing spectral barycenters, we propose to regularize the problem in (5)  subject to A j Φ j = ∆ j + r j for j = 1, . . . , J.
Here, r j = vec(R j ) ∈ C p 2 j , for j = 1, . . . , J, are the vectorizations of the covariances matrices, and A j ∈ C p 2 j ×n are matrix representations of the discretized operators Γ j . Also, the vectors ∆ j ∈ C p 2 j allow for noisy covariance matrix estimates, penalized by the standard 2 -norm, as controlled by the positive penalty parameter γ. The addition of the entropy term, D(M ), renders the problem strictly convex and allows for efficiently computing the barycenter Φ, using the observation in [28], that Sinkhorn iterations are equivalent to a block coordinate ascent in a Lagrange dual of the entropy regularized transport problem. This then allows for increasing the dimensionality of the sensor fusion problem, i.e., making it feasible to include also range estimation or considering higher spatial dimensions. The extended Sinkhorn iterations, constituting a globally convergent solution algorithm for (7), are detailed in the following proposition.
Proposition 1 A block coordinate ascent scheme for the dual of the barycenter problem (7) is to iteratively perform 1. For j = 1, . . . , J, let λ j be the solution to Algorithm 1 Sinkhorn-Newton method for the optimal mass transport barycenter problem with partial information. Require: Initial guess λ j , v j , and let u j = exp(A j λ j / ) for j = 1, . . . , J while Sinkhorn not converged do for j=1:J do while Newton not converged do

Update
v j = y./(K T u j (λ j )) for j = 1, . . . , J, Here, K = exp(−C/ ), and exp(·), ./, and denotes the elementwise exponential function, elementwise division, and elementwise product, respectively. The barycenter solution is then given by Φ = y, where y is the limit point of the iterations. Note that we here convert all complex quantities to an equivalent real-valued representation such that the vectorized covariances, r j , are mapped according to where Re(·) and Im(·) denote the real and imaginary parts, respectively, and where ∆ j and the discretized operators A j are converted as to be consistent with this.
A proof of the proposition may be found in Appendix A. It may be noted that for the noise-free case with complete observations of the underlying spectra, i.e., ∆ j = 0 and A j = I, as well as J = 2, one obtains the standard Sinkhorn iterations.  Fig. 1. Spectral estimate as given by the proposed estimator in (7). The alignment error is 6.7 degrees.
The roots of (8) are determined by the means of a Newton method, using the Jacobian The complete estimation method for computing the spectral barycenter is summarized in Algorithm 1.

NUMERICAL RESULTS
In this section, we illustrate the proposed method's behavior in sensor fusion scenarios for localization in 2D, demonstrating the robustness to array alignment errors that the OMT criterium induces. Specifically, we consider two uncorrelated sources impinging on two asynchronous sensor arrays; one ellipsoidal shaped array consisting of 8 sensors, and one linear array consisting of 7 sensors. To simulate misalignment between the two arrays, we introduce an unknown rotation to the ellipsoidal array. The scenario, as shown in Figure 1, is then used in a simulation study. Varying the misalignment, i.e., the rotation angle, between 0 and 10 degrees, we for each considered angle perform 100 Monte Carlo simulations, were the locations of the signal sources are randomized uniformly on the square [−0.5, 0.5] × [−0.5, 0.5]. We let the source signals be uncorrelated circularly symmetric Gaussian white noises and add an uncorrelated Gaussian noise to the sensors. The source signals have variance 100 and the variance of the sensor noise is 1. The wavelength of the impinging waves is twice that of the smallest sensor spacing in the linear array.
Estimating the covariance matrices of the separate sensor arrays using the sample covariance matrix and 500 signal snapshots, we apply the proposed estimator in (7) to the covariance matrix estimates and record the distance from the location estimates to the actual target positions. Throughout, we use the parameters γ = 0.01 and = 0.005, and the cost function , for grid points x k , x ∈ R 2 . As comparison, we also obtain estimates using the noncoherent MUSIC and MVDR estimators, as described in [18], as well as the least-squares (LS) estimator from [32], and the non-coherent SPICE estimator from [19]. Figures 1 and 2 show the estimates obtained in one simulation using the proposed and MUSIC estimators, respectively. The misalignment angle is 6.7 degrees. As can be seen, for the proposed estimator, the misalignment only causes a small perturbation on the estimated source locations, demonstrating the ability of the OMT formulation to perform smooth interpolation on the underlying domain. In contrast, the MUSIC (pseudo-) spectral estimate contains spurious peaks, as well as larger devi-ations from the true source locations. The robustness of the proposed estimator is further demonstrated in Figure 3, displaying the results from the simulations. From this, it may be noted that the proposed estimator deviates slightly more from the ground truth for the non-perturbed case as compared to the other methods. However, the robustness to misalignment is significantly higher; the rate of error increase is considerably lower than for the other methods.
A. PROOF OF PROPOSITION Proof 1 In the following, we denote the set {M 1 , . . . , M J } by M, and similarly for ∆, λ and µ. The Lagrangian of (7) with the dual variables λ and µ is Minimizing the Lagrangian with respect to M j gives the expression for the mass transport matrix solutions, with K = exp(−C/ ), u j = exp(A T j λ j / ) and v j = exp(µ j / ), for j = 1, . . . , J. Further, minimizing (9) with respect to the variables Φ and ∆, one gets a dual problem to where L(Φ * , M * , ∆ * , λ, µ) j = (λ T j r j − e (A T j λj ) T / e −C/ e µj / − 1 4γ λ T j λ j + n 2 ).
A block coordinate ascent method is to maximize (11) iteratively with respect to λ and µ. Maximization with respect to λ j requires r j − A j diag(e −C/ e µj / )e A T j λj / − 1 2γ λ j = 0, which gives the first step (8). To maximize (11) with respect to µ j , consider the corresponding Lagrangian with dual variable y, given by Maximization with respect to µ j requires −diag(e −C T / e A T j λj / )e µj / + y = 0 ⇔ v j = y./(K T u j ).
Further, due to the constraint in (11), it follows This completes the second step of the method. Finally, with (10), we get the expression for the barycenter marginal


Passive Copy-Move Tamper Detection Methods for Digital Images 5883

5882 Published By: Blue Eyes Intelligence Engineering & Sciences Publication Retrieval Number: C4730098319/19©BEIESP DOI:10.35940/ijrte. C4730.098319 Abstract: Copy-move tamper discovery in digital image is a significant research application in forensic investigation. Developing an imperative and reliable means to discover copy-move tampering in order to guarantee the authenticity of digital images is a solitary research area in image processing. This category of tampering is performed to hide some unwanted information or duplicate certain region of image. This article introduces a generic algorithmic skeleton to copy-move tamper discovery Experiments are carried out on three states of art copy-move tamper detection methods with respect to three different data sets. Empirical results indicates that PCA based method is better than Frequency transformation and Zernike moment based method in terms of detection accuracy in all three data sets. This article also identified the key challenges and opportunities in copy-move tamper detection.

I. INTRODUCTION
Copy-move is a familiar category of image exploitation which is performed along with different kind of geometric transformations and image distortions to match the irregularities among genuine and tampered region of an image [1]. In general, tampered patterns will combine the background details of an image, thus the manual detection of this tampering is a critical task. Hence, development of competent counterfeit discovery scheme becomes a demanding task due to its higher influence. There exist different methods in literature to discover this forgery but accuracy, false positive rate, type and size of image features, robustness against different geometric operations (translation, rotation, scaling) and image distortions (Noise addition, blurring, brightness change, color reduction and JPEG compression) get varied [2] [3]. Organization of this article includes section II which introduces materials and methods in support of tamper finding. Section III depicts the empirical outcome on existing tamper detection methods. Section IV depicts the discussion and Challenges in Copy-Move tamper recognition. Conclusion of the manuscript is presented in section V,

A. Algorithmic framework for Copy-move tamper detection
Discovering copy-move tampering requires profound analysis of local patterns or regions in an image. One preliminary approach is to split the image of dimension M X N into different chunks (b X b) and comparing these chunks for similarity estimation. Dimension of the block will not exceed dimension of assumed tampered region. Each block of pixels are analyzed by means of four methods such as frequency transformation based method, dimensionality reduction method, moment and key point technique. Since the tampered image segments are in same image there exists a correlation between these two image segments. This kind of correlation can be detected with different feature matching and extensive search algorithms [4] [5]. Figure 1 exemplifies this type of manipulation and figure 2 illustrate the whole algorithmic framework for this category.

) Exhaustive search
The elementary way to perceive copy-move tampering is to employ exhaustive search. The image has been shifted in to clockwise direction and shifted image has been compared with Original image as shown in figure 3. If P ij is a pixel value of a gray scale image with dimension M X N at location (i, j), the subsequent disparity is calculated with exhaustive search, (1) Where k=0, 1,……M-1, l=0,1,….N-1 designed for all i and j. Comparing P ij by means of its clockwise shift [k,l] is similar to comparing P ij by its recurring shift [k ' , l ' ] where k ' =M-k and l ' =N-l. Hence, it is sufficient to check merely clockwise shift [k,l] with 1≤k≤M/2, 1≤l≤N/2 . The technique is straightforward and successful for tiny scale images. Comparing segments with its cyclic shift leads high computational cost and not suitable for medium and large scale images [6].  For images with M rows and N columns the correlation coefficients between image segments are computed using the formula, is the auto correlation function, f(x, y) indicates the significance of pixels at position (x,y). Variable i and j represent pixel lags. During the initial stage suspected region of the image is extracted and whole image is divided to number of sub-regions of size same as the extracted region. Autocorrelation function is performed among every sub-region of the image. Sub region having highest correlation is identified as tampered region. This method is a viable option for small scale images and not suitable for large scale images since, it has large computational overhead and often not succeeded to perceive the forgery [7]. Both Exhaustive and auto correlation methods are traditional methods having its own advantages and implementation issues.

III. RESULTS
Most of the methods reported in the literature encompass challenges like accuracy in tamper detection and false detection. In this section, three well known state of art methods such as frequency transformation based, PCA based and Zernike moment based methods are successfully implemented using MATLAB and experiments were carried out with three different bench mark data sets namely Kodak image dataset [8], CoMoFoD image data base [9] and Muhammad et al. dataset [10]. Detail descriptions of the above three data sets are presented in the table 1.

A. Data Set 1) KODAK image data set
The first bench mark KODAK image data set restrains 24 uncompressed true color images. All images are portable Network Graphic (PNG) format of size 768 X 512 pixels. Images regions are arbitrarily tampered with adobe Photoshop image editing tool. 2) The CoMoFoD image data base CoMoFoD data set comprise 260 authentic and 260 tampered images, out of those 200 images are under small category of size 512 X 512 pixels and 60 images are under big group of dimension 3000 X 2000 pixels. Different transformations are applied to tampered images and classified into five groups as follows, 1. Translation -The copied region is shifted into fresh position not including further transformations. 2. Rotationtampered segment is spin with different angles also pasted in another location. 3. Scalingtampered segment size is enlarged or diminished also pasted in another position. 4. Distortiontampered province is deformed by various image distortions and interprets in to another position. 5. Combination -Copied region is deformed by two or more image transformations before it is moved to new location.

Different image distortions
Various image distortions performed on original and tampered images and their parameters are exposed in table 2.
False negative FN designates the amount of tampered images classified as authentic. F coalesce precision and recall value while it is close [11] P * R F = 2 * P + R

B. Frequency transformation based copy move tamper detection method
In this method, suspected input image is alienated into overlapping blocks and Discrete Cosine Transformation (DCT) coefficients are extracted from each overlapping blocks. Extracted coefficients are quantized with user specific parameter value Q for robust representation of image data and discover similar pattern in an image. Parameter value Q is indicate quality factor, the Q value is low false matches are high and the Q value is high false matches are low. Quantized Where (i1, i2) and (j1, j2) indicates position matching blocks. The algorithm finds normalized shift vector whose incidence surpass user specific threshold (T). Value of T is related to size of smallest fragment. Detection result of this method is shown in figure 4 [6]. It is optimal in the sense of energy compaction and well known method for reducing dimensionality of data. PCA renovate higher dimensional data into meaningful lower dimensional representation by means of Eigen vector decomposition of co-variance (or correlation) matrix. If data is projected on to eigenvector corresponding first largest Eigen value (Sometimes known as first principle component) then, it will reduce the data in to one dimension. If the data is projected on to Eigen vector corresponding to first and second largest Eigen value then, it will reduce the data into two dimensions. After reducing the dimension from d to k (k<d) the new set of coordinates are the linear combination of the original variables. First principle component has largest variance, second principle component has next largest variance and orthogonal to first one. This process will be continued until the dimension vector is reduced to k. When the dimension vector reached the maximum value all other variables are neglected without loss of information. Any square, symmetric, non-singular matrix can be transformed to a diagonal matrix using Eigen vectors. PCA based technique is competent and stout practice to perceive copy move tampering. In this technique, the suspected gray scale image is alienated to small preset size blocks of B pixels and PCA is applied to all blocks to diminish the image dimension representation. Tampered regions are detected by matching similar PCA coefficients. This method is used to identifies patterns in large set of data but infeasible for large size image blocks. Detection result of this method against simple and multiple copy-move tampering is shown in figure 5 [12].

D. Zernike moment based tamper detection
Moments and its invariant properties are widely used in image registration, reconstruction and its related fields. Different types of moments are introduced in the literature. In this method, Seung et al. used Zernike moments to discover most common nature of counterfeit such as copy-move; it shows better result than other moments. In order to perceive the manipulation the suspected image f with dimension M X N is alienated into overlapping chunks of dimension Zernike moments are extorted for every block B ij , where i and j designate the initial points. Extracted Zernike moments are included in the two dimensional matrix and lexicographical sorting is applied. Analogous rows are closest to each other therefore Euclidean distance between two similar rows are calculated to match the tampered image. This method effectively detect simple, multiple tampering and rotated manipulation. Detection results are shown in figure  6 [13].

IV. DISCUSSION
Experimental results demonstrate that the above three methods are efficient scheme to discover copy-move tampering. Performances of all three methods are evaluated under simple copy-move tampering at image level. Overall performance comparison of these methods using different bench mark data base are shown in figure 7 to 9. The performance curve of PCA method out performs other two methods in all three dataset. Computational complexity is yet another problem which is associated with block matching. Lexicographical matrix sorting is identified as a key factor related with computational difficulty. Table 3 shows image representation and feature dimension of three methods. Out of three methods, Zernike method has reduced feature dimension but accuracy of the method is low over PCA method.   Developing a distinctive as well as reliable method to perceive copy move tampering with different angel of rotation and various scaling factors.  Detection of the forged region which is affected by illumination.  Implementing detection system with dynamic view point change.  Design of robust method to reduce false matches and improve the accuracy while perceive copy-move manipulation.
 Developing a method to distinguish the copy-move forgery in an image with dissimilar image deformation such as Blurring, brightness change, contrast adjustment, color reduction, noise addition and JPEG compression with different quality factor

V. CONCLUSION
This article focused the experimental analysis of three well known copy-move tamper detection techniques and its implementation issues. The general algorithmic framework and two traditional approaches in copy-move tamper detection are reviewed. From the existing methods, it is found that the reliable techniques need to be developed to identify the copy-move counterfeit when the image underwent some geometric transformation like translation, rotation and scaling. Apart from the transformation, the presence of forgery region and size of image will affect the performance of tamper detection method. Empirical result shows that a method having good accuracy rate may not have reduced features and vice versa. Hence, development of an optimized method for accurate tamper detection with reduced feature set will have higher influence in image based applications such as security and authentication.


extfieldsofstudy:
Computer Science
Mathematics
Biology
Engineering
Physics
Medicine

text:
Inter- and intra- patient ECG heartbeat classification for arrhythmia detection: a sequence to sequence deep learning approach

Electrocardiogram (ECG) signal is a common and powerful tool to study heart function and diagnose several abnormal arrhythmias. While there have been remarkable improvements in cardiac arrhythmia classification methods, they still cannot offer acceptable performance in detecting different heart conditions, especially when dealing with imbalanced datasets. In this paper, we propose a solution to address this limitation of current classification approaches by developing an automatic heartbeat classification method using deep convolutional neural networks and sequence to sequence models. We evaluated the proposed method on the MIT-BIH arrhythmia database, considering the intra-patient and inter-patient paradigms, and the AAMI EC57 standard. The evaluation results for both paradigms show that our method achieves the best performance in the literature (a positive predictive value of 96.46% and sensitivity of 100% for the category S, and a positive predictive value of 98.68% and sensitivity of 97.40% for the category F for the intra-patient scheme; a positive predictive value of 92.57% and sensitivity of 88.94% for the category S, and a positive predictive value of 99.50% and sensitivity of 99.94% for the category V for the inter-patient scheme.). The source code is available at https://github.com/SajadMo/ECG-Heartbeat-Classification-seq2seq-model.

I. INTRODUCTION
An Electrocardiogram (ECG) is a common non-invasive tool to record heart activities and detect different abnormalities in heart functionality. Classification of the arrhythmic heartbeats in the ECG signal can be a challenging and time-consuming task for a physician, therefore, such heartbeat hand-annotating is often prone to error. This calls for automatic heartbeat classification methods that are able to diagnose arrhythmic heartbeats in real-time with high accuracy. Several machine learning algorithms such as support vector machines (SVM), multilayer perceptron (MLP), reservoir computing with logistic regression (RC) and decision trees have been utilized for arrhythmia detection [1], [2], [3], [4], [5], [6], [7], [8]. These shallow machine learning methods for ECG processing usually follow three main steps, including 1) signal pre-processing, which includes noise removal methods, heartbeat segmentation, etc; 2) feature extraction; and 3) learning/classification. Even though these methods with hand-engineered features and applying noise removal techniques have achieved acceptable performances, deep learning approaches (i.e., automated feature extractions) have shown impressive results in various domains ranging from computer vision and reinforcement learning to natural language processing [9], [10], [11], [12], [13], [14], [15], [16] as well as more applicable outcomes in biomedical signal processing [17], [18], [19], [20].
One of the main limitations of the current heartbeat classification methods including shallow and deep machine learning methods is their poor performance when dealing with imbalanced datasets. In particular, they attain a low positive predictive value and sensitivity for the classes with lower sample size in the dataset. For instance, the majority of existing ECG analysis techniques achieve a low sensitivity in the MIT-BIH arrhythmia database for ventricular escape beat (S) and fusion of ventricular, and normal beat (F) classes. Furthermore, most previously reported works in the literature have been evaluated based on intrapatient paradigm rather than the inter-patient scheme which is an obviously more realistic scenario to prevent training and test the model using the samples from the same patients. Therefore, although some of these methods achieved good accuracies using the intra-patient scheme, their results are not reliable as their evaluation process was biased [21].
As mentioned above, the conventional arrhythmia classification systems can be generally divided into two categories of inter-patient and intra-patient paradigms in terms of their evaluation mechanism. In intra-patient paradigm, the training and evaluation datasets can include heartbeats from the same patients, while in inter-patient paradigm, a more realistic evaluation mechanism is used where the heartbeat sets for test and training come from different individuals. One of our aims in this paper is to evaluate the proposed method with both the paradigms.
Inspired by the aforementioned issues with the previous works, this paper proposes a novel and effective approach for automatic ECG-based heartbeat classification by leveraging a sequence to sequence deep learning method and an oversampling method named Synthetic Minority Over-sampling Technique (SMOTE) to address the aforementioned challenge with minority classes. The proposed model is evaluated using inter-patient and intra-patient paradigms where it achieves the best results compared to the existing works in the literature.
The rest of this paper is organized as follows. Section II introduces the database utilized in this study. Section III describes the proposed method. Section IV presents the experimental setup and shows the achieved results by the proposed method along with a performance comparison to the state-of-the-art algorithms. Finally, Section V concludes the paper.

II. DATASET
In this study, we used the PhysioNet MIT-BIH Arrhythmia database to evaluate the performance of our proposed method [22], [23]. The MIT-BIH dataset includes the ECG signals for 48 different subjects recorded at the sampling rate of 360Hz. Each record contains two ECG leads; ECG lead II and lead V1. Usually, the lead II is used to detect heartbeats in the literature. Similarly, here in all experiments, we have applied ECG lead II. This database is recommended by the American association of medical instrumentation (AAMI) [24], since it includes the five essential arrhythmia groups as described in Table I.
We considered two main paradigms of inter-patient and intra-patient to evaluate the proposed model. In the intra-patient paradigm, two sets of data samples (beats) are chosen randomly as training and test samples in which there may be the heartbeat samples of the same patient in the training and test sets. While, in the inter-patient paradigm, the training and test set are constructed from different patients, following the protocol proposed by de Chazal et al. [21]. In this method, the heartbeats from the MIT-BIH database (44 records based on AAMI) are divided into two sets of records: DS1 = {101, 101, 106, 108, 109,112,114,115,116,118,119,122,124,201,203,205,207,208,209 is no concern about including the heartbeats from the same patient in both training and test sets.

III. METHODOLOGY
In the following sections, we present a detailed description of our proposed novel model to automatically classify each heartbeat of a given ECG signal.

A. Pre-processing
The input of this method is a sequence of ECG beats. In order to extract the ECG beats from a given ECG signal, we follow a few simple steps: 1) Normalizing the given ECG signal to the range of between zero and one.
2) Finding the set of t waves regarding the ECG R-peaks of its corresponding annotation file in the MIT-BIH Arrhythmia database.  3) Splitting the continuous ECG signal to a sequence of heartbeats based on the extracted t waves and assigning a label to each heartbeat based on the annotation file. 4) Resizing each heartbeat to a predefined fixed length (280 samples). We would like to note that these pre-processing steps for beat extraction are very simple and do not involve any form of filtering or noise removal methods.

B. The architecture
The sequence to sequence models have shown very impassive results in neural machine translation applications, nearly similar to human-level performance [25]. The architecture of sequence to sequence networks is usually composed of two main parts of the recurrent neural network (RNN) encoder and decoder. In this study, we leverage an RNN sequence to sequence model along with a convolutional neural network (CNN) to perform a heartbeat detection task. Fig. 1 illustrates the proposed network architecture for automatic beat classification. The CNN consists of three consecutive one-dimensional convolutional layers. The first layer is composed of 32 1-D convolution filters with a kernel size of 2 × 1 and a stride 1, followed by a Rectified Linear Unit (ReLU) non-linearity. The second layer consists of 64 1-D convolution filters with a kernel size of 2 × 1 and a stride 1, again followed by an ReLU. Finally, the third layer is comprised of 128 1-D convolution filters with a kernel size of 2 × 1 and a stride 1, followed by a rectifier non-linearity. Each convolutional layer except the last layer is also followed by a max pooling layer of pooling region of size 2 × 1 with a stride 1. At each time-step of training/testing the model, a sequence (size of maxtime) of ECG beats is fed into the CNN in order for feature extraction. The last convolutional layer outputs the maxtime of F feature maps of size k × 1 (e.g, here, we reached 128 feature maps 3 × 1). In the end, each beat of the input sequence is associated with a vector c ∈ R d . Figure 2 depicts the detailed network.
The sequence to sequence model is designed based on the encoder-decoder abstract ideas. The encoder encodes the input sequence, while the decoder computes the category of each beat of the input sequence. The encoder is actually composed of long short-term memory (LSTM) units, which is also called the many to one LSTM. The LSTMs can capture the complex and long short-term context dependencies between the inputs and the targets [26]. This is due to the fact that they capture non-linear dependencies on entire observation when predicting a target. The (time) sequence of input feature vectors herein are fed to the LSTMs and then the last hidden state calculated by the LSTM is considered as the encoder representation and is used to initialize the fist hidden state of the decoder, as depicted in Fig. 1  We have utilized the bidirectional recurrent neural network (BiRNN) units in the network architecture instead of the standard LSTM (i.e., RNN). Standard RNNs are unidirectional, hence they are restricted to the use of the previous input state. To address this limitation, the BiRNN have been proposed [27], which can process data in both forward and backward directions. Thus, the current state has access to previous and future input information simultaneously. The BiRNN consists of a forward network and a backward network. The input sequence is fed in normal time order, t = 1, ..., T for the forward network, and in reverse time order, t = T, ..., 1 for the backward network. Finally, the weighted sum of the outputs of the two networks is computed as the output of the BiRNN. This mechanism can be formulated as follow: where ( − → h t , − → b ) are the hidden state and the bias of the forward network, and ( ← − h t , ← − b ) are the hidden state and the bias of the backward network. Also, x t and y t are the input and the output of the BiRNN, respectively. The decoder is used to generate the target sequence beat by beat. Similar to the encoder, the building block of the decoder is an LSTM, but a many-to-many LSTM. The decoder gets a new representation of the input sequence generated by the encoder to initialize its hidden state. It also takes the same given target shifted by one and started with a special feature vector < GO > as input. We should note that  Figure 2: A diagram of convolutional layers used in the proposed model.
the input (the shifted target) is just used during the training phase and is not applied for the testing phase. Then, a softmax is applied to the output of the LSTM to convert it to a vector of probabilities p ∈ R C , where C represents the number of categories (i.e., the heartbeat types) and each element of p indicates the probability of each class in the category.

A. Experimental setup
The performance of the proposed heartbeat detection method has been evaluated using the MIT-BIH arrhythmia database using both inter-patient and intra-patient evaluation methods. Heartbeat category distribution of extracted beats from the MIT-BIH arrhythmia is not uniform and the numbers of normal beats are much greater than other categories. Machine learning approaches usually have trouble learning when one class dominates the others. To cope with this, the dataset has been oversampled to nearly reach a balanced number of beats in each category (group). To this end, we used the synthetic minority over-sampling technique (SMOTE) which generates the synthetic data points by considering the similarities between existing minority samples [28].
Similar to several works in the literature, we evaluated the arrhythmia classifier on four cardiac cycles including N (normal and bundle branch block beats), S (supraventricular ectopic beats), V (ventricular ectopic beats), and F (fusion of N and V beats) (for more details refer to Table I). For the intra-patient scheme, we extracted totally 101,290 heartbeats, including 90,494 beats for the category N, 2,777 for the category S, 7,217 for the category V and 802 for the category F. To reduce the impact of class imbalance problem, the SMOTE technique was applied to upsample the categories with the small factions of heartbeats. We have trained the proposed model with 80% of the dataset and evaluated it with the remaining 20%. While for the inter-patient scheme, we had 45,796 beats of the category N, 941 beats of the category S, and 3,780 beats of the category V in DS1 (i.e., training set), and 44,196 beats for the category N, 1,836 beats for the category S, and 3,216 beats for the category V in DS2 (i.e., test set). Again, the SMOTE approach was utilized to compensate the number of heartbeats of the categories with the lower heartbeats.
The augmented data (i.e., synthetic data generated based on the original available data) should not be considered in the validation phase as the information of the validation set may exist in the training set. For instance, although the work reported in [18], followed this approach, we believe that the results are not quite accurate. It is worth mentioning that we generated new data samples by the SMOTE algorithm after splitting the heartbeats into the training and test datasets. Indeed, the oversampling was only performed on the training dataset for the intra-patient paradigm (and on DS1 for the inter-patient paradigm) and the information in the test dataset was not used to create synthetic data points. Therefore, the evaluation process is reliable and generalizable.
Four main measures usually were considered in the literature to evaluate the performance of heartbeat classification techniques including the sensitivity (SEN), positive predictive value (PPV), specificity (SPEC) (1 -False Positive Rate (FPR)) and accuracy (Acc) as defined in following: P P V = T P/(T P + F P ) (5) where TP (True Positive), TN (True Negative), FP (False Positive) and FN (False Negative) indicate the number of heartbeats correctly labeled, number of heartbeats correctly identified as not correspond to the heartbeats, number of heartbeats that incorrectly labeled, and number of heartbeats which were not identified as the heartbeats that they should have been, respectively. The network was trained for a maximum of 300 epochs and the initial LSTM hidden and cell states were set to zero. All network weights were updated by the RMSProp algorithm with mini batches of size 20 and a learning rate of α = 0.001.  [17] 93

B. Results and Discussion
The results are presented for two evaluation scenarios of intar-patient paradigm, in which the training and test sets were randomly chosen from all available patients' heartbeats, and the inter-patient paradigm, in which the training and test have been performed on the heartbeats of extracted from DS1 and DS2, respectively (i.e., no common individual in test and training sets). Table II presents a comparison of heartbeat classification results for the proposed method and the existing algorithms, considering intra-patient scheme. As it is clear, our sequence to sequence model outperforms all state-of-the-art algorithms significantly in terms of all evaluation metrics including overall accuracy, sensitivity, positive predictive value and specificity for the considered groups, N, S, V, F. In addition, it is worth noticing that the work done by Acharya et al. [18] used artificially augmented dataset to build the model, and their model evaluation was performed using augmented data, while we evaluated our model on real data samples without including any augmentations in the test set.
As confirmed by the results, our proposed method can provide a robust solution for class imbalance problem as one of the key challenges in dealing with medical data, which is due to the limited availability of abnormal classes compared to the normal classes in biomedical datasets. It is shown in Table II that our model achieves remarkable outcomes for the category  F with only 802 heartbeats and the category S with 2,777 heartbeats. We also validated our method using the more realistic evaluation method of inter-patient for heartbeat classifiers based on using DS1 set for training and DS2 set for testing [21], [32]. Table III presents the performance comparison between our proposed method and several state-of-the-art methods using MIT-BIH arrhythmia database where the inter-patient paradigm is considered. As it can be seen from the table, overall the proposed method has better performance for classifying all heartbeat categories compared to works listed in Table III. In spite of the low number of S heartbeats in the training set, our proposed method obtained significant evaluation results.
The proposed method is generic in nature and it is expected to achieve a promising performance in several biomedical applications dealing with class imbalance problem. In the proposed classification structure, first, the CNN extracts a set of meaningful features of the given ECG heartbeats (Note, the model was trained by added synthetic data generated by SMOTE algorithm to the available samples to compensate the number of small categories like F). Then, the encoder maps the features to new feature representations, capturing temporal patterns, and finally, the decoder takes the feature representations and produces the outputs (i.e., the labels for each heartbeat of the input sequence), considering complex context dependencies between the inputs and the targets.

V. CONCLUSION
In this study, we presented a novel and effective automatic heartbeat classification/annotation, considering intar-and interpatient schemes and validated its performance using the MIT-BIH arrhythmia database. The proposed method leverages the ability of deep convolutional neural network and encoder-decoder network in which we have used a bidirectional recurrent neural network as its building blocks. According to the results, the suggested method significantly outperforms the existing algorithms in the literature for both intar-patient paradigm and inter-patient paradigm. Furthermore, the proposed method can be applied to several biomedical applications such as sleep staging where there are strong dependencies between each stage and sufficient data are not available. Also, the proposed network with a low number of parameters (i.e., with a maximum size of 5.5MB) can be used with wearable devices.


The problem of neural networks communication

In spite of the successful application of artificial neural networks (ANNs) for the solution of multiple problems (forecasting, language translation, image classification, voice recognition etc.), ANNs are still autonomous entities incapable of communication or exchange of their knowledge. Meanwhile, the ability to communicate is critical for further development of methods of artificial intelligence. We propose and test several methods of communication and knowledge fusion of ANNs. These methods do not require the presence of the initial training data and use only the internal parameters of ANNs. We propose generative iterative and non-iterative methods of ANNs communication. Noniterative methods show the classification accuracy similar to that provided by an ensemble of ANNs. The accuracy of generative methods is similar to a network trained on the joint dataset.

Introduction
The ability to communicate and exchange their knowledge is crucial for artificial neural networks (ANNs) and is necessary for further development of artificial intelligence and self-learning systems. Modern ANNs, however, behave as completely autonomous units incapable of any information exchange. Additional difficulty comes from the fact that ANNs store the acquired knowledge in a distributed way amongst the millions of internal parameters (weights). Such an implicit representation of knowledge complicates its interpretation, interferes with its update, and obstructs its exchange with other neural networks. If joint use of knowledge is needed, there exists two current approaches: 1) training new neural network with a joint dataset; 2) using an ensemble of networks to produce joint prediction [1]. However, there could be situations when both of these approaches become unacceptable: previously used training data may become unavailable; storing and maintaining an ensemble of ANNs may become too resource demanding. Ideally, neural networks should be able to exchange their knowledge without any external help; the result of this exchange would be a single ANN containing the fused knowledge (in the following we denote this as fANN).
The problem of knowledge fusion of several ANNs may be divided into two separate subproblems: extraction of knowledge from one neural network and its incorporation into another ANN. The problem of knowledge extraction is essential for understanding of the functioning of neural networks. Without this understanding, ANNs cannot be used in such critical applications as IOP Publishing doi: 10.1088/1742-6596/1368/5/052033 2 transportation and health care. The problem of new knowledge incorporation into a neural network is closely related to the problem of catastrophic interference: when trained on new piece of data without an access to previous datasets, ANN quickly forgets all the previously learned knowledge [2]. Thus, the development of algorithms for communication and knowledge exchange between neural networks would simultaneously solve other important problems of machine learning. At the same time, one can try to find methods that do not require explicit extraction and embedding of knowledge, but instead operate with internal parameters of networks.
This paper is devoted to the development of methods of knowledge fusion of neural networks that do not lead to the forgetting of previously learned content and do not require previous training data. The need in the integration of knowledge from different ANNs without an access to previous training data may appear in various situations. One example is medical expert systems, which training data could be confidential, but joint knowledge can lead to better diagnosis. Another example is autonomous vehicles, which can benefit from the experience of peers. Previously, the problem of knowledge fusion was solved only for ANNs with similar weights [3] and for Support Vector Machines [4]. We developed both noniterative methods and generative iterative approaches for the solution of knowledge fusion problem for ANNs trained independently. Noniterative methods allow fast knowledge fusion; generative methods provide more accuracy.
Let us note that the problem of knowledge fusion should not be confused with the term "knowledge transfer" (or "transfer learning") [5]. According to transfer learning concept, the weights of deep network trained on one set of data can be used for training on another dataset from a similar area to speed up the training [5]. Thus, transfer learning is more concerned with the learning of new task rather than maintaining all the knowledge. Examples of transfer learning can be found in such areas as computer vision [6], language processing [7] etc. In contrast to knowledge transfer, knowledge fusion tries to integrate all the knowledge of initial networks into a single ANN. Another difference with knowledge transfer is the absence of training data, which makes the problem of knowledge fusion even more challenging.

Problem statement
In its simplest form, the problem of knowledge fusion can be stated as follows. Two ANNs (A and B) with a similar structure are trained on a classification task by some standard method (for example, by backpropagation). Each neural network is trained on its own dataset. After the training is complete, each neural network contains different knowledge. It is necessary to modify the parameters of one of the networks (for example, by changing the weights) in such a way that the resulting network would be able to reproduce the response of both of the initial networks as if it was trained on a joint dataset (figure 1). During fusion of knowledge, we assume that all the previous training data are no longer accessible. In experiments performed in this work, we use one-hot encoding in the output layers of ANNs. This encoding assumes activation of one output neuron corresponding to a specific class. Thus, the For the fusion of knowledge of deep networks, we additionally utilize the concept of transfer learning [5]. A typical architecture of deep networks consists of a convolutional part, which extracts the features from raw data, and a fully connected classifier on top that performs the separation of the features into classes. The process of training of deep networks consisted of several stages. First, the weights of some network trained on large image dataset (for example, ImageNet [8]) were used for initialization of the weights of the convolutional part. The large set of classes of ImageNet database makes the pretrained network an ideal candidate for transfer learning. At the second stage, the bottom convolutional part was supplemented by an untrained classifier, which was later trained on various classification tasks (for example, classification of animals or urban objects). At the third stage, the convolutional part of each network underwent fine tuning leading to the final slight modification of the weights. As a result, we obtained a set of deep networks with different set of weights and containing different knowledge.
Next section describes the methods for the merging the knowledge contained in several networks in more details.

Method of weights summation
The method of weights summation (WS) is based on a probabilistic process of summation of a large number of random variables (weights of ANNs) [9]. The weights of the network with a fused knowledge F can be obtained by direct summation of weights of networks and : In spite of its simplicity, this algorithm results in a rather high classification accuracy of the final fused network. We propose the following explanation of the efficiency of WS method. For some feature vector belonging to class of network , the presynaptic activity in fANN is determined by the following expression: In principle, the trained neural networks are robust under the perturbation of the weights. Thus, if the second term at the right-hand side of equation (1) is significantly smaller that the first one, the fused network will still be able to perform correct classification. To ensure the inequality , ≫ , for every node (where ( , ), are the weights connected to the node ), several conditions should be met. First, networks and should be well trained in a sense that they produce noticeable presynaptic activity in the presence of native features and low activity in the presence of foreign ones. Second, the mean of the weights' probability distribution within one layer should be zero 〈 〉 〈 〉 0. Indeed, the weights of one layer , , , and components of the feature vector , may be considered as random variables taken from some distribution. For network , feature vector contains independent random variables. The product of independent random variables ( ) can be rewritten as where is the number of components in the feature vector and average 〈 , 〉 is taken over presynaptic weights of corresponding layer. It follows from equation (2), that the term , will be close to zero if the expected value of weights distribution is zero 〈 〉 ≈ 0. The property 〈 〉 0 is valid for common training techniques [10] and is maintained in current research.
In our experiments, the WS algorithm was tested on shallow fully-connected feed-forward networks and on deep convolutional networks. The experiments with fully-connected networks were performed on a set of handwritten digits MNIST [11]. The whole dataset consists of 60000 training images of size 28 by 28 and 10000 test samples corresponding to 10 classes. We arbitrarily divided the whole dataset into two equal sets containing 5 classes each. The intensity of every input image was normalized to be in the range [0,1]. Each network and was trained on data corresponding only to the target classes of particular network. Thus, each network never saw the samples corresponding to the classes of another network and could not meaningfully classify corresponding images. The training algorithm was implemented in Keras 2.2 with the TensorFlow backend. Both networks and had an input layer with 784 nodes corresponding to the number of pixels in input images, one hidden layer with 800 nodes, and output layer with 5 nodes corresponding to the number of classes. All hidden layers had an additional bias term. The nodes in the hidden layer had ReLU activation functions. For the output layer, we used softmax activation; the cross-entropy was used as a loss function. Before training, all the weights were initialized randomly from a normal distribution with zero mean with a standard deviation 0.05. L2 regularization was used during training with L2 0.001. Adam iterative method was used for training with all the parameters set to the Keras' default values. In each experiment, a batch size of 200 was used. Early stopping was used while testing the accuracy on a validation set consisting of 12000 samples. We quantify the efficiency of classification of each network through accuracy (percent of correctly classified images). Simulations were repeated 10 times with random initialization of weights and with random subdivision of digits into classes to accumulate proper statistics.
The accuracy of the original networks A and B on partial datasets achieved 0.983±0.002. After the fusion of knowledge by weights summation, the accuracy of fANN turned out to be 0.864±0.035. This number is certainly smaller than the accuracy of a network trained on a joint dataset 0.968±0.002, but still much larger than a chance level.
We also compared the accuracy of fANN to the classification accuracy of an ensemble of networks A and B. The ensemble prediction was chosen on a base of maximum activity among all output neurons of all networks. The accuracy of fANN appeared to be just a little lower than the accuracy of an ensemble 0.900±0.021. However, one should remember that fANN replaces the ensemble of networks with just one neural net.
Experiments with deep networks were performed on separate classes of images of ImageNet project [8] (www.image-net.org). Considering the variety of the objects indexed in the project, one can construct the networks trained to classify either similar or dissimilar classes of objects. In particular, for current experiment, two "indigenous" neural networks were trained to classify 11 African and 11 Australian animals (table 1). The third "urban" neural network was trained to classify 11 objects pertinent to urban life (table 1). Each class had approximately 1500 samples.
As in experiments with MNIST dataset, it was assumed that the original trained networks A and B observed only the objects of their own classes and never saw the objects from the classes of other networks. The aim was to embed the knowledge contained in one neural network into another one in such a way that the resulting network can recognize the objects (e.g., animals or household item) from the previously unknown category even without ever "seeing" the objects from those categories before.
For the architecture of networks and , we chose VGG-16 model [12] included in Keras version 2.2 with its weights pre-trained on 1000 ImageNet classes. Those 1000 classes include large variety of objects, but, still, do not contain all the classes used in our experiments.
With the concept of knowledge transfer, the pretrained convolutional part of VGG-16 network was a good starting point for the training on previously unknown objects. In the pretrained VGG-16 network, we kept only its bottom convolutional part that we used as universal feature extractor. On top of this convolutional part, we placed a custom classifier consisting of one hidden layer with 128 nodes and one output layer with 11 neurons corresponding to the number of classes of the networks being fused. Dropout was used for the weights connecting hidden to output nodes of the classifier. ReLU activation function was used for the hidden layer, softmax activation was used for the output layer.
The input to VGG-16 network has to be a fixed-size 224 × 224 RGB image that was cropped out of a center of input raw images. Augmentation of dataset was not used. All the images were preprocessed by scaling the intensity of every image to the range [0,1]; 15% of the images were dedicated to validation and another 15% were used for the test. For the training, we used 64 pictures batch, and the training was performed with Adam optimizer on a computer with GPU support (Ubuntu 16.04 х64, Intel core i3, 6GB RAM, Nvidia GeForce GTX 1050 Ti 4GB).
The training procedure for the initial networks A and B was performed in several steps. First, the bottom convolutional part was frozen, and the top classifier was initialized with random weights and trained with small learning rate (10 -5 ) over 4 epochs. Then we froze the weights of the classifier and unfroze the top block (3 layers) of VGG-16 convolutional part. This block was trained over 2 epochs with the same small learning rate. Such a training procedure resulted in networks and different not only in their top fully connected parts, but also in their convolutional layers. As a result, the networks trained to classify African and Australian animals achieved 0.951±0.005 and 0.933±0.004 accuracy, correspondingly. The network trained to classify urban objects achieved 0.955±0.003 accuracy. The main source of classification errors was the mixing between classes with similar features. For example, the African network often misclassified chimpanzees and gorillas. The Australian network, in its turn, misclassified kangaroos and wallabies.
As the weights of the convolutional parts of both networks A and B originated from the same set of pretrained weights, during fusion by WS method, those weights were averaged according to the approach adopted by Smith & Gashler [3]. The weights of the classifiers on the top of networks were obtained from an independent set of weights and thus were added together according to the weights summation method. Experiments demonstrated that the classification accuracy of the network with merged knowledge depends on the similarity of the features of the objects, classified by original networks and . For example, if both networks and were trained to recognize the objects of the same category, then the network with fused knowledge would confuse the objects with similar sets of features. Figure 2 shows the result of knowledge fusion between two networks in the case when they were trained to recognize African and Australian animals. fANN is confident (classification accuracy is of the order of 90%) about animals with unique features, for example, zebras, but confuses similar animals from different continents (for example, emus and ostriches). An average classification accuracy of fANN is 76%. For comparison, the classification accuracy of the same set of classes by an ensemble of networks A and B is 81%.
In the case of very different set of features of the objects recognized by networks and (for example, in the case of living and inanimate objects), the classification accuracy of fANN increases. Figure 3 shows the result of knowledge fusion from two networks in a case when one of the networks was trained to recognize African animals and the other was trained to recognize urban objects. An average classification accuracy of fANN was 81%, which is better than in previous example. The classification accuracy by an ensemble of networks was 88%.

Method of elastic weight consolidation
The second noniterative method that we used for the fusion of knowledge of ANNs was based on an idea of selective modification of weights. If one changes the weights of networks A and B until these weights coincide with additional condition that loss functions ( ) and ( ) increase insignificantly The increase of total loss function F ( F ) ( F ) + ( F ) is minimal if we change only less important weights for current task. Following paper [13], we call this method Elastic Weight Consolidation (EWC).
To estimate the change of the loss function, one can decompose it into Taylor's series with respect to weights. Here is an example of such a decomposition for network : Here are the optimal weights of network trained on corresponding dataset; the summation in equation (3) is performed over all the weights. As the number of parameters reaches millions in modern neural networks, the computation and the storage of the whole matrix of second order partial derivatives (Hessian) 2 ( ) ⁄ is impossible. Following paper [13], we approximated the whole Hessian by its diagonal part: where , ≡ 2 2 ⁄ . Expressions similar to (3) and (4) should be written for network as well. Factors , , , in equation (4) express the importance of weights , , , . Usually, one needs to calculate these coefficients in advance during the training of neural network. It is easy to show that the weights of fANN F that minimize the loss function F can be calculated with the following formula In principle, equation (5) solves the problem of fusion knowledge of two ANNs. However, weights , of independently trained networks are not necessarily located in a vicinity of one another, and decomposition (3) may become insufficient. In attempt to bring closer the values of and , we performed an "alignment" of networks by rearranging hidden neurons. The contribution to the overall loss function F from a pairing of a node of network and a node of network is To find the optimal pairing of neurons from networks A and B, one needs to solve, so called, Assignment Problem [14]. The coefficients compose the cost matrix that serves as input for the Assignment Problem algorithm.
For experiments with EWC method, we used the networks with the same architecture as in WS experiments. Similarly, the experiments were performed on shallow networks with MNIST data and on deep convolutional networks with ImageNet data. In spite of heavier computations, EWC demonstrated performance very similar to the one of weights summation method. This tells us that noniterative methods of knowledge fusion are limited in achieving high classification accuracy. An alternative approach can be developed on the base of generative models.

Generative methods
Recently, neural networks were developed that can generate data of the same type and similar probabilistic distribution ( ) as the initial training data. One of such ANNs are generativeadversarial networks (GANs) [15]. Another type of generative neural networks is variational autoencoder (VAE) [16].
The idea behind the generative knowledge fusion is in the following. The knowledge of the trained classifier is its ability to predict the probability ( | ) of an object belonging to class . If another network somehow learns to reproduce this probability than it acquires the knowledge of the first network. Moreover, original dataset { } is not necessary; one just need to reproduce probability ( |̃) for any ̃ sampled from probability ( ). Generative models can provide both samples ̃ from probability ( ) and soft targets ̃ representing probabilities of these samples to belong to certain classes. For our experiments, we used variational autoencoders [16]. The architecture of VAE (figure 4a) was similar to the one used by van de Ven & Tolias [17] for prevention of catastrophic interference in sequential learning. Regular classifier was supplemented by generative part (figure 4a). The generative part includes standard VAE components [16]: sampler from latent distribution ( , ) and decoder (figure 4a). The loss function of the resulting VAE total consisted of three parts: categorical crossentropy classification error class , binary cross-entropy reconstruction loss recon , and latent variable regularization term KL [16,17]  To prevent catastrophic forgetting, synthetic data should be generated not only from knowledge donor but also from knowledge acceptor (figure 4b).
The experiments with generative approach were performed on MNIST dataset. The classifier of the same structure as the one described in section 3.1 was supplemented by generative part (figure 4a) consisting of two layers with 8 nodes each encoding latent distribution mean and variance , one hidden layer with 800 neurons (decoder) and ReLU activation function, and output layer with sigmoid activation function and 784 nodes corresponding to the dimension of the original images.
As in previous experiments, MNIST dataset was arbitrarily divided into two parts containing 5 classes each. Networks A and B were trained over 40 epochs with Adam optimizer to properly classify and generate images of digits. The classification accuracy of networks A and B achieved 0.983±0.004. After the networks A and B have being trained, 60000 artificial samples ̃ with corresponding soft labels ̃ were generated from the decoder by supplying it with random variable sampled from standard normal distribution. The examples of generated images are shown in figure 5. As usually happens with images generated by VAEs, the digits appear somewhat blurred. The generated datasets {̃,̃}, {̃,̃} were supplied to network B (figure 4b); network B was trained on this combined dataset over 40 epochs. The accuracy of the retrained network B tested on complete MNIST test set was 0.929±0.006. This is considerably higher than the accuracy of noniterative knowledge fusion, higher than the accuracy of an ensemble of networks, and it is getting close to the accuracy of a network trained on joint dataset (0.9698±0.0006). The shortcoming of generative method however is larger network, that should incorporate generative part and longer training time.

Conclusion
In this work, we demonstrated the possibility of knowledge fusion of several neural networks without an access to the original training data. The developed methods are applicable both for fully-connected and for convolutional networks. The methods were tested in classification tasks using publicly available datasets.
Noniterative methods (summation of weights and elastic weight consolidation) demonstrate the possibility of fast fusion of knowledge from several neural networks. The classification accuracy of the network with fused knowledge depends on the similarity of features of objects in the fused datasets: the classification accuracy increases when fusing knowledge from different areas. The classification accuracy of fANN obtained with WS and EWC methods is similar to one obtained from an ensemble of networks.
Generative iterative methods are slower than noniterative ones, however, they show higher accuracy in knowledge fusion. The accuracy of fANN obtained with generative methods is similar to the accuracy of the network trained on a joint set of training data. The quality of knowledge fusion is determined by the quality of reconstruction of original data by generative models and can be increased by improving VAE's generation quality [18] or by using different generative models [19]. In this paper, we presented generative methods only for shallow networks. In the future, we plan to extend the generative approaches for deep convolutional networks.
In this work, we considered examples of knowledge fusion between two ANNs. However, the proposed methods allow simple generalization for several networks and for networks with dissimilar architecture.
In addition to applications in machine learning, development of methods of interaction and information exchange between neural networks will be useful for other scientific areas such as neuro-