# Study on Parameter Optimization for Support Vector Regression in Solving the Inverse ECG Problem

Authors: Mingfeng Jiang; Shanshan Jiang; Lingyan Zhu; Yaming Wang; Wenqing Huang; Heng Zhang
Journal: Computational and Mathematical Methods in Medicine
 (2013)

---

## Abstract

The typical inverse ECG problem is to noninvasively reconstruct the transmembrane potentials (TMPs) from body surface potentials (BSPs). In the study, the inverse ECG problem can be treated as a regression problem with multi-inputs (body surface potentials) and multi-outputs (transmembrane potentials), which can be solved by the support vector regression (SVR) method. In order to obtain an effective SVR model with optimal regression accuracy and generalization performance, the hyperparameters of SVR must be set carefully. Three different optimization methods, that is, genetic algorithm (GA), differential evolution (DE) algorithm, and particle swarm optimization (PSO), are proposed to determine optimal hyperparameters of the SVR model. In this paper, we attempt to investigate which one is the most effective way in reconstructing the cardiac TMPs from BSPs, and a full comparison of their performances is also provided. The experimental results show that these three optimization methods are well performed in finding the proper parameters of SVR and can yield good generalization performance in solving the inverse ECG problem. Moreover, compared with DE and GA, PSO algorithm is more efficient in parameters optimization and performs better in solving the inverse ECG problem, leading to a more accurate reconstruction of the TMPs.


---

## Body

## 1. Introduction

The inverse ECG problem is to obtain myocardial transmembrane potential (TMPs) distribution from body surface potentials (BSPs), thus noninvasively imaging the electrophysiological information on the cardiac surface [1, 2]. Generally, approaches to solving this inverse ECG problem can be relied on potential-based model, including epicardial, endocardial, or transmembrane potentials, which is used to evaluate the potential values on the cardiac surface [3] at certain time instants. Moreover, the cardiac electrophysiological information is closely associated with the transmembrane potentials (TMPs) of the myocardial cells. Compared to body surface potentials (BSPs) recordings, TMPs can provide more detailed and complicated electrophysiological information. In this study, we focus on implementing the reconstruction of TMPs from BSPs.To study inverse ECG problems, various numerical methods have been proposed. In the last decades, regularization methods have been engaged for dealing with the inherent ill-posed property. The regularization techniques include truncated total least squares (TTLS) [4], GMRes [5], and the LSQR [6], which require an appropriate selection of regularization parameters so as to relax the ill-posedness of the inverse ECG problem and to produce a well-posed problem. However, the robustness and the quality of the inverse solution are not always guaranteed; despite that they can more or less deal with the geometry and measurement noises for the inverse ECG problems. In addition, during the solution procedure, the inverse ECG problem can be treated as a regression problem with multi-inputs (BSPs) and multi-outputs (TMPs). Therefore, an alternative method, support vector regression (SVR) method [7], was proposed to solve the inverse ECG problem. Compared with conventional regularization methods, SVR method can produce more accurate results in terms of reconstruction of the transmembrane potential distributions on epi- and endocardial surfaces [1, 8]. SVR is an extension version of support vector machine (SVM) to solve nonlinear regression estimation problems, which aims at minimizing an upper bound of the generalization error instead of the training error [9] by adhering to the principle of structural risk minimization (SRM).Although SVR is a powerful technique to solve the nonlinear regression problem, it has received less attention in relation to the inverse ECG problems, due to the fact that SVR algorithm is sensitive to users’ defined free parameters. The involved hyperparameters of SVR model consist of penalty parameterC, insensitive loss function parameter ε, and the parameter σ for kernel function. Inappropriate parameters in SVR can lead to overfitting or underfitting problems. How to properly set the hyperparameters is a major task, which has a significant impact on the optimal generalization performance and the SVR regression accuracy [10], especially when it comes to predicting the diagnosis of cardiac diseases, because even a slight improvement of prediction accuracy could have a significance impact on the patient’s diagnosis [11, 12].In early development stage of the algorithm, a grid search optimizing method [13] and cross-validation method [8] are employed to optimize the hyperparameters. However, these methods are computationally expensive and data intensive. Recently, a number of new algorithms have been proposed for the optimization of the SVR parameters [14, 15]. For example, Wang et al. [16] has proposed a hybrid load forecasting model combining differential evolution (DE) algorithm and support vector regression to deal with the problem of annual load forecasting. In the analysis of predicting tourism demand, the study [17] applies genetic algorithm (GA) to seek the SVR’s optimal parameters and then adopts the optimal parameters to construct the SVR models. The experimental results demonstrate that SVR outperforms the other two neural network techniques. Another study [18] tries a new technology, particle swarm optimization (PSO), to automatically determine the parameters of SVR, then applies the hybrid model (PSO-SVR) to grid resource prediction, and the experimental results indicate high predictive accuracy.In this paper, the above mentioned optimization algorithms (GA, DE, and PSO) are all adopted to  dynamically optimize the hyperparameters of SVR model in solving the inverse ECG problem. For convenience, the SVR model with GA parameter selection is referred to as GA-SVR method, and the other two are termed DE-SVR and PSO-SVR, respectively. In this paper, we attempt to investigate which one is the most effective in reconstructing the cardiac TMPs from BSPs, and a full comparison of the performance for solving the inverse ECG problem will be evaluated.


## 2. Theory and Methodology

### 2.1. Brief Overview of the SVR

In this section, the basic SVR concepts are concisely described; for detailed description, please see [19, 20]. Suppose a given training data of N elements {(xi,yi),i=1,2,…,N}, where xi denotes the ith element in n-dimensional space; that is, xi={x1i,…,xni}∈Rn,   and  yi∈R is the output value corresponding to xi. According to mathematical notation, the SVR algorithm builds the linear regression function as the following form:
(1)f(x,ω)=(ω·φ(x)+b),φ:Rn⟶F,ω∈F,
where ω and b are the slope and offset coefficients; φ(x) denotes the high-dimensional feature space, which is nonlinearly mapped from the input space x. And the previous regression problem is equivalent to minimizing the following convex optimization problem [see (2)]:
(2)min12∥ω∥2subjecttoyi-ωφ(xi)-b≤εωφ(xi)+b-yi≤ε.
In this equation, an implicit assumption is that a function f essentially approximates all pairs (xi,  yi) with ε precision, but sometimes this may not be the case. Therefore, by introducing two additional positive slack variables ξi and ξi*, the minimization is reformulated as the following constrained optimization problem [shown in (3)]:
(3)minR(ω,ξ,ξ*)=12∥ω∥2+C∑i=1N(ξi+ξi*)s.t.yi-(ω,φ(x))-b≤ε+ξi*(ω,φ(x))+b-yi≤ε+ξiξi,ξi*≥0,i=1,2,…,N,ε≥0,
where the parameter C is the regulator which is determined by the user, and it influences a tradeoff between an approximation error and the weights vector norm ||ω||; ξi and ξi* are slack variables that represent the distance from actual values to the corresponding boundary values of ε-tube.According to the strategy outlined by Schölkopf and Smola [10], by applying Lagrangian theory and the KKT condition, the constrained optimization problem can be further restated as the following equation:
(4)f(x,ω)=∑i=1N(αi-αi*)K(xi,x)+bs.t.∑i=1N(αi-αi*)=0.
Here αi and αi* are the Lagrange multipliers. The term K(xi,x) is defined as the kernel function. The nonlinear separable cases could be easily transformed to linear cases by mapping the original variable into a new feature space of high dimension using K(xi,x). The radial basis function (RBF) was applied in the study, which has the ability to universally approximate any distribution in the feature space. With an appropriate parameter, RBF usually provides a better prediction performance, so it is adopted in this study as shown in (5):
(5)K(xi,xj)=exp(-∥xi-xj∥22σ2),
wherexiandxiare input vector spaces and σ2 is the bandwidth of the kernel function.In the above equations, there exist three hyper-parameters to be determined in advance, that is, the penalty parameterC, insensitive parameter ε, and the related kernel function parameters σ2. They heavily affect the regression accuracy and computation complexity of SVR. The penalty parameter C controls the degree of punishing the samples whose errors go beyond the given value. The insensitive parameter ε controls the width of the ε-insensitive zone used to fit the training data. The value of ε can enhance the generalization capability; with the increase of ε, the number of support vectors will decrease, and the algorithmic computation complexity will also reduce. The bandwidth σ of the kernel function has a great influence on the performance of learning machine.In this study, three optimization methods, that is, GA, DE, and PSO, are presented to determine the optimal hyper-parameters of the SVR model. According to [15], the general range of C,  σ2, and ε has been given. In the trial operation, we narrowed it to avoid blindness in the optimization process. In this paper, the set of hyperparameter (C,σ2,ε) is initialized in the given range C∈[0,10000],  σ2∈[0,2], and ε∈[0,0.0001], where optimization methods (GA, DE, and PSO) are to seek the global optimal solutions.In order to calculate the error uniformly, we adopt the same fitness function which plays a critical role in measuring these algorithms performance. The fitness function is determined as follows:(6)minf=∑i=1nt(|ai-pi|/ai)nt*100%,
where nt is the number of training data samples, ai is the actual TMPs of train data, and pi is the predicted TMPs. The solution with a smaller fitness f of the training dataset has a better chance of surviving in the successive generations. The main tool, LIBSVM, was used for training and validating the SVR model [21].


### 2.2. Genetic Algorithm (GA) Optimization Method

The GA [22] is a biologically motivated optimization technique guided by the concepts of biological evolution and Darwin’s principle of survival of the fittest. It is a computer model of an evolution of a population of artificial individuals. In this study, for the specific optimizing problem of hyperparameters (C,σ2,ε), the process is defined as follows.Step 1 (Initialize parameters for GA).
The population size NP is equal to 20; the dimension of parameter vectorsD is 3. The termination criterion is set as follows: the number of iterations is set as 30, and the fitness tolerance value is set as 0.001. The probabilities of selection (Stochastic universal sampling), crossover (multipoint crossover) and mutation (mutation operator of the Breeder Genetic Algorithm) that were used herein are 0.9, 0.8, and 0.05, respectively.Step 2 (encode chromosomes).
According to the possible range of parametersC,σ2, and ε given before, the GA utilizes binary encoding method, and each parameter is encoded by 20 bits of binary number. Therefore, the search space is defined as the solution space in which each potential solution is encoded as a distinct chromosome.Step 3.
The parametersC,σ2, and ε of each individual are used to build the SVR model. With the BSPs of the training data, the cardiac TMPs are reconstructed. Then the performances of individuals in the generation are evaluated by the specific objective fitness function according to (6). The individual with the minimum fitness value will be selected, and then the chromosome of selected individual is preserved as the best result.Step 4.
The basic genetic search operators comprise selection, mutation, and crossover, which are applied orderly to obtain a new generation where the new individual(C,σ2,ε) with the best performance is retained.Step 5.
The new best fitness value will be compared to that of the best result, and then select the better one to update the best result.Step 6.
This process will not come to an end until the termination criterion is met, and the best chromosome is presented as the final solution; otherwise, go back to Step4.


### 2.3. Differential Evolution (DE) Optimization Method

Differential evolution [23] is a population-based and parallel direct search method which is used to approximate the global optimal solution. As is mentioned above, the optimization of hyperparameter (C,σ2,ε) can be transformed into solving the minimization of the fitness function Min: f(x),  x=xiG,  iG=1,2,…,  and NP is three-dimensional parameter vectors including (C,σ2,  and  ε). Subsequently, each generation evolves by employing the evolutionary operators involving mutation, crossover, and selection operations to produce a trail vector for each individual vector [24], and the detailed evolutionary strategies can be described as follows.Step 1 (initialize DE parameters).
The population size NP is set as 20 [25], and the dimension of parameter vectors D is 3. The termination criterion is set as follows: the number of iterations is set as 30, and the fitness tolerance value is set as 0.001. The mutation factor F is selected in [0.5, 1], and the crossover rate CR is selected in [0,1].Step 2 (initialize population).
SetG=0. Generate an NP*Dgeneration which consists of individuals (C,σ2,  and  ε) with uniform probability distribution random values within the control variable bounds of parameter space.Step 3.
With the parameters, train the SVR model, and forecast the training datasets. Then calculate the fitness target of all the individuals in the generation, and record the minimum values of the fitness function and set(C,σ2,  and ε) correspondingly as the xbest.Step 4 (mutation operation).
For each target vectorxi,G, an associated mutant vector is generated according to
(7)vi,G+1=xbest+F[(xr1,G-xr2,G)+(xr3,G-xr4,G)].
The random indexesr1,r2,r3, andr4have to be distinguished from each other and from the current trial index i.Step 5 (crossover operation).
In order to increase the diversity of the perturbed parameter vectors, crossover is introduced by the following formula:(8)uji,G+1={vji,G+1if(randb(j)≤CR)orj=rnbr(i)xji,G+1if(randb(j)>CR)andj≠rnbr(i),
where ui,G+1=(u1i,G+1,u2i,G+1,u3i,G+1), randb(j)∈[0,1] is the jth evaluation of a uniform random number generator, and rnbr(i) is a randomly chosen integer in the range in [1, 3] to ensure that ui,G+1 gets at least one element from vi,G+1.Step 6 (selection operation).
After the fitness values of the new generation being calculated, then selecting the new best individual asubest,G+1, the selection operation is performed. To decide whether or not it could replace the xbest, the fitness value of f(ubest,G+1) is compared to that of f(xbest). The operation is given as follows:
(9)xbest={ubest,G+1f(ubest,G+1)≤f(xbest)xbestf(ubest,G+1)>f(xbest).Step 7.
Termination condition checking: if max iterationsG is met or the fitness value of f(xbest) is reached within the fitness tolerance value, return the recorded global optimal parameter (C,σ2, and ε); otherwise, go to Step 4.


### 2.4. Particle Swarm Optimization (PSO) Optimization Method

The Particle swarm optimization (PSO) is a bioinspired stochastic optimization technique which was basically developed through simulating social behaviour of birds and insects that keep living by maintaining swarm actions. A particle is considered as a bird in a swarm consisting of a number of birds, and all particles fly through the searching space by following the current optimum particle to find the final optimum solution of the optimization problem [26].The detailed experimental procedure for selecting the best hyperparameters (C,σ2, and ε) is listed as follows.Step 1.
Particle initialization and PSO parameters setting [27]: set the PSO parameters including the scope of C,σ2, and ε; the number of particles is 20; particle dimension is 3. The termination criterion is set as follows: the number of iterations is set as 30 and the fitness tolerance value as 0.001.Step 2.
Randomly generate primal population of random particles (C,σ2, and ε) and velocities inside the searching space.Step 3.
Construct SVR model with the parameters in each particle, and perform the prediction with the training data. Then evaluate the regression accuracy based on the defined fitness function. Then maintain the record of the best performance with the minimum error as the global best solution.Step 4.
Particle manipulations: suppose thatXi(k) and Vi(k) are the particle and velocity at step k, and they are updated into their new values Xi(k+1) and Vi(k+1) at step k+1 by applying the two equations below:
(10)Vi(k+1)=ωVi(k)+c1r1(pi(k)-Xi(k))+c2r2(pg(k)-Xi(k)),Xi(k+1)=Xi(k)+Vi(k+1)Δt,
where Vi∈[Vimin,Vimax],  Δt is the unit time step value, and ω is the inertia weight. pi(k) and pg(k) are the best coordinates of particle number i and the whole swarm until step k, respectively. Positive constants c1 and c2 are learning rates while r1 and r2 are random numbers between 0 and 1.Step 5.
Train the SVR model using the particles (C,σ2, and ε) in the new generation and forecast the training datasets. After evaluating the fitness function, pick the new best particle as the local best.Step 6.
Compare the local best and global best, and choose the better one to update the global best.Step 7.
Stop condition checking: if a termination criterion (the number of iterations) is met or the fitness value of the global best is reached within the extent predefined, return the recorded global best as the best parameter (C,σ2, and ε); otherwise, go to Step 4.


### 2.5. Simulation Protocol and Data Set

The SVR model is tested with our previously developed realistic heart-torso model [11, 12]. In this simulation protocol, a normal ventricular excitation is illustrated as an example to calculate the data set for the SVR model. The considered ventricular excitation period from the first breakthrough to the end is 357 ms, and the time step is 1 ms, and thus 358 BSPs φB and TMPs φm temporal data sets are calculated. Sixty data sets (φB and φm) are chosen at times of 3 ms, 9 ms, 15 ms, …, 357 ms after the first ventricular breakthrough which will be used as testing samples to assess the SVR models prediction capability and robustness. The rest 298 in 358 data sets are employed as the training samples for the training and optimal parameter selection procedures. Moreover, during each ventricular excitation period, 412 potentials on the BSPs and 478 potentials on the TMPs are captured. That is to say, the matrix of BSPs is divided into the training data (298 × 412) and testing data (60 × 412), and the matrix of TMPs is also divided into the training data (298 × 478) and testing data (60 × 478). Each SVR regression model is built by using the mapping relations between BSPs training data (298 × 412) and one TMPs training point (298 × 1). Based on all the training data, 478 different regression models were built by using GA-SVR, DE-SVR, and PSO-SVR, respectively. Using the corresponding SVR model, we can reconstruct the corresponding TMPs (60 × 1) from the BSPs testing data (60 × 412). For all testing samples, we can reconstruct the TMPs from the all 478 regression SVR models.


## 2.1. Brief Overview of the SVR

In this section, the basic SVR concepts are concisely described; for detailed description, please see [19, 20]. Suppose a given training data of N elements {(xi,yi),i=1,2,…,N}, where xi denotes the ith element in n-dimensional space; that is, xi={x1i,…,xni}∈Rn,   and  yi∈R is the output value corresponding to xi. According to mathematical notation, the SVR algorithm builds the linear regression function as the following form:
(1)f(x,ω)=(ω·φ(x)+b),φ:Rn⟶F,ω∈F,
where ω and b are the slope and offset coefficients; φ(x) denotes the high-dimensional feature space, which is nonlinearly mapped from the input space x. And the previous regression problem is equivalent to minimizing the following convex optimization problem [see (2)]:
(2)min12∥ω∥2subjecttoyi-ωφ(xi)-b≤εωφ(xi)+b-yi≤ε.
In this equation, an implicit assumption is that a function f essentially approximates all pairs (xi,  yi) with ε precision, but sometimes this may not be the case. Therefore, by introducing two additional positive slack variables ξi and ξi*, the minimization is reformulated as the following constrained optimization problem [shown in (3)]:
(3)minR(ω,ξ,ξ*)=12∥ω∥2+C∑i=1N(ξi+ξi*)s.t.yi-(ω,φ(x))-b≤ε+ξi*(ω,φ(x))+b-yi≤ε+ξiξi,ξi*≥0,i=1,2,…,N,ε≥0,
where the parameter C is the regulator which is determined by the user, and it influences a tradeoff between an approximation error and the weights vector norm ||ω||; ξi and ξi* are slack variables that represent the distance from actual values to the corresponding boundary values of ε-tube.According to the strategy outlined by Schölkopf and Smola [10], by applying Lagrangian theory and the KKT condition, the constrained optimization problem can be further restated as the following equation:
(4)f(x,ω)=∑i=1N(αi-αi*)K(xi,x)+bs.t.∑i=1N(αi-αi*)=0.
Here αi and αi* are the Lagrange multipliers. The term K(xi,x) is defined as the kernel function. The nonlinear separable cases could be easily transformed to linear cases by mapping the original variable into a new feature space of high dimension using K(xi,x). The radial basis function (RBF) was applied in the study, which has the ability to universally approximate any distribution in the feature space. With an appropriate parameter, RBF usually provides a better prediction performance, so it is adopted in this study as shown in (5):
(5)K(xi,xj)=exp(-∥xi-xj∥22σ2),
wherexiandxiare input vector spaces and σ2 is the bandwidth of the kernel function.In the above equations, there exist three hyper-parameters to be determined in advance, that is, the penalty parameterC, insensitive parameter ε, and the related kernel function parameters σ2. They heavily affect the regression accuracy and computation complexity of SVR. The penalty parameter C controls the degree of punishing the samples whose errors go beyond the given value. The insensitive parameter ε controls the width of the ε-insensitive zone used to fit the training data. The value of ε can enhance the generalization capability; with the increase of ε, the number of support vectors will decrease, and the algorithmic computation complexity will also reduce. The bandwidth σ of the kernel function has a great influence on the performance of learning machine.In this study, three optimization methods, that is, GA, DE, and PSO, are presented to determine the optimal hyper-parameters of the SVR model. According to [15], the general range of C,  σ2, and ε has been given. In the trial operation, we narrowed it to avoid blindness in the optimization process. In this paper, the set of hyperparameter (C,σ2,ε) is initialized in the given range C∈[0,10000],  σ2∈[0,2], and ε∈[0,0.0001], where optimization methods (GA, DE, and PSO) are to seek the global optimal solutions.In order to calculate the error uniformly, we adopt the same fitness function which plays a critical role in measuring these algorithms performance. The fitness function is determined as follows:(6)minf=∑i=1nt(|ai-pi|/ai)nt*100%,
where nt is the number of training data samples, ai is the actual TMPs of train data, and pi is the predicted TMPs. The solution with a smaller fitness f of the training dataset has a better chance of surviving in the successive generations. The main tool, LIBSVM, was used for training and validating the SVR model [21].


## 2.2. Genetic Algorithm (GA) Optimization Method

The GA [22] is a biologically motivated optimization technique guided by the concepts of biological evolution and Darwin’s principle of survival of the fittest. It is a computer model of an evolution of a population of artificial individuals. In this study, for the specific optimizing problem of hyperparameters (C,σ2,ε), the process is defined as follows.Step 1 (Initialize parameters for GA).
The population size NP is equal to 20; the dimension of parameter vectorsD is 3. The termination criterion is set as follows: the number of iterations is set as 30, and the fitness tolerance value is set as 0.001. The probabilities of selection (Stochastic universal sampling), crossover (multipoint crossover) and mutation (mutation operator of the Breeder Genetic Algorithm) that were used herein are 0.9, 0.8, and 0.05, respectively.Step 2 (encode chromosomes).
According to the possible range of parametersC,σ2, and ε given before, the GA utilizes binary encoding method, and each parameter is encoded by 20 bits of binary number. Therefore, the search space is defined as the solution space in which each potential solution is encoded as a distinct chromosome.Step 3.
The parametersC,σ2, and ε of each individual are used to build the SVR model. With the BSPs of the training data, the cardiac TMPs are reconstructed. Then the performances of individuals in the generation are evaluated by the specific objective fitness function according to (6). The individual with the minimum fitness value will be selected, and then the chromosome of selected individual is preserved as the best result.Step 4.
The basic genetic search operators comprise selection, mutation, and crossover, which are applied orderly to obtain a new generation where the new individual(C,σ2,ε) with the best performance is retained.Step 5.
The new best fitness value will be compared to that of the best result, and then select the better one to update the best result.Step 6.
This process will not come to an end until the termination criterion is met, and the best chromosome is presented as the final solution; otherwise, go back to Step4.


## 2.3. Differential Evolution (DE) Optimization Method

Differential evolution [23] is a population-based and parallel direct search method which is used to approximate the global optimal solution. As is mentioned above, the optimization of hyperparameter (C,σ2,ε) can be transformed into solving the minimization of the fitness function Min: f(x),  x=xiG,  iG=1,2,…,  and NP is three-dimensional parameter vectors including (C,σ2,  and  ε). Subsequently, each generation evolves by employing the evolutionary operators involving mutation, crossover, and selection operations to produce a trail vector for each individual vector [24], and the detailed evolutionary strategies can be described as follows.Step 1 (initialize DE parameters).
The population size NP is set as 20 [25], and the dimension of parameter vectors D is 3. The termination criterion is set as follows: the number of iterations is set as 30, and the fitness tolerance value is set as 0.001. The mutation factor F is selected in [0.5, 1], and the crossover rate CR is selected in [0,1].Step 2 (initialize population).
SetG=0. Generate an NP*Dgeneration which consists of individuals (C,σ2,  and  ε) with uniform probability distribution random values within the control variable bounds of parameter space.Step 3.
With the parameters, train the SVR model, and forecast the training datasets. Then calculate the fitness target of all the individuals in the generation, and record the minimum values of the fitness function and set(C,σ2,  and ε) correspondingly as the xbest.Step 4 (mutation operation).
For each target vectorxi,G, an associated mutant vector is generated according to
(7)vi,G+1=xbest+F[(xr1,G-xr2,G)+(xr3,G-xr4,G)].
The random indexesr1,r2,r3, andr4have to be distinguished from each other and from the current trial index i.Step 5 (crossover operation).
In order to increase the diversity of the perturbed parameter vectors, crossover is introduced by the following formula:(8)uji,G+1={vji,G+1if(randb(j)≤CR)orj=rnbr(i)xji,G+1if(randb(j)>CR)andj≠rnbr(i),
where ui,G+1=(u1i,G+1,u2i,G+1,u3i,G+1), randb(j)∈[0,1] is the jth evaluation of a uniform random number generator, and rnbr(i) is a randomly chosen integer in the range in [1, 3] to ensure that ui,G+1 gets at least one element from vi,G+1.Step 6 (selection operation).
After the fitness values of the new generation being calculated, then selecting the new best individual asubest,G+1, the selection operation is performed. To decide whether or not it could replace the xbest, the fitness value of f(ubest,G+1) is compared to that of f(xbest). The operation is given as follows:
(9)xbest={ubest,G+1f(ubest,G+1)≤f(xbest)xbestf(ubest,G+1)>f(xbest).Step 7.
Termination condition checking: if max iterationsG is met or the fitness value of f(xbest) is reached within the fitness tolerance value, return the recorded global optimal parameter (C,σ2, and ε); otherwise, go to Step 4.


## 2.4. Particle Swarm Optimization (PSO) Optimization Method

The Particle swarm optimization (PSO) is a bioinspired stochastic optimization technique which was basically developed through simulating social behaviour of birds and insects that keep living by maintaining swarm actions. A particle is considered as a bird in a swarm consisting of a number of birds, and all particles fly through the searching space by following the current optimum particle to find the final optimum solution of the optimization problem [26].The detailed experimental procedure for selecting the best hyperparameters (C,σ2, and ε) is listed as follows.Step 1.
Particle initialization and PSO parameters setting [27]: set the PSO parameters including the scope of C,σ2, and ε; the number of particles is 20; particle dimension is 3. The termination criterion is set as follows: the number of iterations is set as 30 and the fitness tolerance value as 0.001.Step 2.
Randomly generate primal population of random particles (C,σ2, and ε) and velocities inside the searching space.Step 3.
Construct SVR model with the parameters in each particle, and perform the prediction with the training data. Then evaluate the regression accuracy based on the defined fitness function. Then maintain the record of the best performance with the minimum error as the global best solution.Step 4.
Particle manipulations: suppose thatXi(k) and Vi(k) are the particle and velocity at step k, and they are updated into their new values Xi(k+1) and Vi(k+1) at step k+1 by applying the two equations below:
(10)Vi(k+1)=ωVi(k)+c1r1(pi(k)-Xi(k))+c2r2(pg(k)-Xi(k)),Xi(k+1)=Xi(k)+Vi(k+1)Δt,
where Vi∈[Vimin,Vimax],  Δt is the unit time step value, and ω is the inertia weight. pi(k) and pg(k) are the best coordinates of particle number i and the whole swarm until step k, respectively. Positive constants c1 and c2 are learning rates while r1 and r2 are random numbers between 0 and 1.Step 5.
Train the SVR model using the particles (C,σ2, and ε) in the new generation and forecast the training datasets. After evaluating the fitness function, pick the new best particle as the local best.Step 6.
Compare the local best and global best, and choose the better one to update the global best.Step 7.
Stop condition checking: if a termination criterion (the number of iterations) is met or the fitness value of the global best is reached within the extent predefined, return the recorded global best as the best parameter (C,σ2, and ε); otherwise, go to Step 4.


## 2.5. Simulation Protocol and Data Set

The SVR model is tested with our previously developed realistic heart-torso model [11, 12]. In this simulation protocol, a normal ventricular excitation is illustrated as an example to calculate the data set for the SVR model. The considered ventricular excitation period from the first breakthrough to the end is 357 ms, and the time step is 1 ms, and thus 358 BSPs φB and TMPs φm temporal data sets are calculated. Sixty data sets (φB and φm) are chosen at times of 3 ms, 9 ms, 15 ms, …, 357 ms after the first ventricular breakthrough which will be used as testing samples to assess the SVR models prediction capability and robustness. The rest 298 in 358 data sets are employed as the training samples for the training and optimal parameter selection procedures. Moreover, during each ventricular excitation period, 412 potentials on the BSPs and 478 potentials on the TMPs are captured. That is to say, the matrix of BSPs is divided into the training data (298 × 412) and testing data (60 × 412), and the matrix of TMPs is also divided into the training data (298 × 478) and testing data (60 × 478). Each SVR regression model is built by using the mapping relations between BSPs training data (298 × 412) and one TMPs training point (298 × 1). Based on all the training data, 478 different regression models were built by using GA-SVR, DE-SVR, and PSO-SVR, respectively. Using the corresponding SVR model, we can reconstruct the corresponding TMPs (60 × 1) from the BSPs testing data (60 × 412). For all testing samples, we can reconstruct the TMPs from the all 478 regression SVR models.


## 3. The Proposed System Framework

### 3.1. Overview of the SVR Parameters Optimization Framework

As shown in Figure1, in the first stage, the original input data is preprocessed by scaling the training data and feature extraction. In the second stage, the hyperparameters (C, σ2, and ε) of SVR method are set carefully by using the GA, DE and PSO optimization methods. Moreover, GA-SVR, DE-SVR, and PSO-SVR are applied to construct the hybrid SVR model, respectively. Finally, according to the effective hybrid SVR models, we can reconstruct the transmembrane potentials (TMPs) from body surface potentials (BSPs) by using these testing samples.Figure 1
The framework of proposed parameters optimization method for SVR method in solving the inverse ECG problem.


### 3.2. Preprocessing the Data Set

#### 3.2.1. Scaling the Training Set

During preprocessing stage, each input variable is scaled in order to avoid the potential value from spanning a great numerical range and to prevent numeric difficulties. Generally, all input values including 358 BSPs and TMPs temporal data sets are linearly scaled to the range (0,1) by the following formula [see (11)]:
(11)φt′=φt-φtminφtmax-φtmin,
where φt is the original potential value of each time t, φt′ is the scaled potential value, φtmax is the maximum potential value of each time t, and φtmin is the minimum potential value of each time t.


### 3.3. Feature Extraction by Using KPCA

The kernel principal component analysis (KPCA) is a nonlinear feature extraction method which is one type of a nonlinear PCA developed by generalizing the kernel method into PCA. Specifically, the KPCA initially maps the original inputs into a high-dimensional feature spaceF using the kernel method and then calculates the PCA in the high-dimensional feature space F. The detailed theory of KPCA can be consulted in the paper [28].Feature extraction is a significant task for preprocessing the original input data in developing an effective SVR model. As is investigated in our previous study [12], the SVR with feature extraction by using KPCA has superior performances to that of using PCA or the single SVR method in reconstructing the TMPs. In this paper, the kernel principal component analysis (KPCA) was proposed to implement the feature extraction. And the dimension of BSPs dataset is reduced from 412 to 200 properly, which can reduce the dimensions of the inputs and improve the generalization performances of the SVR method.


## 3.1. Overview of the SVR Parameters Optimization Framework

As shown in Figure1, in the first stage, the original input data is preprocessed by scaling the training data and feature extraction. In the second stage, the hyperparameters (C, σ2, and ε) of SVR method are set carefully by using the GA, DE and PSO optimization methods. Moreover, GA-SVR, DE-SVR, and PSO-SVR are applied to construct the hybrid SVR model, respectively. Finally, according to the effective hybrid SVR models, we can reconstruct the transmembrane potentials (TMPs) from body surface potentials (BSPs) by using these testing samples.Figure 1
The framework of proposed parameters optimization method for SVR method in solving the inverse ECG problem.


## 3.2. Preprocessing the Data Set

### 3.2.1. Scaling the Training Set

During preprocessing stage, each input variable is scaled in order to avoid the potential value from spanning a great numerical range and to prevent numeric difficulties. Generally, all input values including 358 BSPs and TMPs temporal data sets are linearly scaled to the range (0,1) by the following formula [see (11)]:
(11)φt′=φt-φtminφtmax-φtmin,
where φt is the original potential value of each time t, φt′ is the scaled potential value, φtmax is the maximum potential value of each time t, and φtmin is the minimum potential value of each time t.


## 3.2.1. Scaling the Training Set

During preprocessing stage, each input variable is scaled in order to avoid the potential value from spanning a great numerical range and to prevent numeric difficulties. Generally, all input values including 358 BSPs and TMPs temporal data sets are linearly scaled to the range (0,1) by the following formula [see (11)]:
(11)φt′=φt-φtminφtmax-φtmin,
where φt is the original potential value of each time t, φt′ is the scaled potential value, φtmax is the maximum potential value of each time t, and φtmin is the minimum potential value of each time t.


## 3.3. Feature Extraction by Using KPCA

The kernel principal component analysis (KPCA) is a nonlinear feature extraction method which is one type of a nonlinear PCA developed by generalizing the kernel method into PCA. Specifically, the KPCA initially maps the original inputs into a high-dimensional feature spaceF using the kernel method and then calculates the PCA in the high-dimensional feature space F. The detailed theory of KPCA can be consulted in the paper [28].Feature extraction is a significant task for preprocessing the original input data in developing an effective SVR model. As is investigated in our previous study [12], the SVR with feature extraction by using KPCA has superior performances to that of using PCA or the single SVR method in reconstructing the TMPs. In this paper, the kernel principal component analysis (KPCA) was proposed to implement the feature extraction. And the dimension of BSPs dataset is reduced from 412 to 200 properly, which can reduce the dimensions of the inputs and improve the generalization performances of the SVR method.


## 4. Results and Discussion

### 4.1. The Result of Parameters Selection

In this study, the GA, DE, and PSO algorithms were proposed to seek the corresponding optimal hyperparameter (C,σ2,  and ε) of the 478 SVR models. Table 1 displays the average value and the standard deviations of three parameters in 478 SVR models by using GA, DE, and PSO optimization algorithms. From it, we can find that the gaps among these optimal parameters are obvious by using three different optimization methods. The standard deviations of ε are neglected because they are very small.Table 1
The average and the standard deviations of optimal hyperparameters (C, σ2, and ε) by using GA, DE, and PSO optimization algorithms.
Parameter
Method
GA
DE
PSO
C
1
447.61
±
594.32
7624.14
±
913.67
5.3
1
±
3.02
σ
2
0
.
64
±
0.55
1
.
59
±
1.38
1
.
2
2
±
1.02
ε
1.4
8
e
-
00
5
1.80
e
-
00
5
1.98
e
-
005


### 4.2. Comparison of Reconstruction Accuracy

This study compares the reconstruction performance based on the proposed three intelligent optimization algorithms for selecting the hyperparameters of SVR. By using GA-SVR, DE-SVR, and PSO-SVR hybrid models, two types of experiments for reconstructing the TMPs have been conducted. The first experiment reconstructs the TMPs for one representative source point (the 50th point) on the epicardium over 60 testing times, as is depicted in Figure2. Since the gaps among these three methods are so slight that we have to add a statistical analysis to validate the results. The mean values and standard deviations of these 60 testing errors between the simulated TMPs and the reconstructed TMPs by using GA-SVR, DE-SVR, and PSO-SVR methods are investigated, as shown in the Table 2. From Figure 2 and Table 2, the PSO-SVR outperforms the GA-SVR and DE-SVR methods in reconstructing the TMPs for one representative source point over all the testing times. Furthermore, the DE-SVR is superior to the GA-SVR method in reconstructing the TMPs.Table 2
The mean values and standard deviations of these 60 testing errors between the simulated TMPs and the reconstructed TMPs by using GA-SVR, DE-SVR, and PSO-SVR methods at one representative source point (the 50th point on the epicardium).
Indexes
Method
GA-SVR
DE-SVR
PSO-SVR
Mean value
2.96
2.69
1.55
Standard deviation
1.92
1.67
1.03Figure 2
The temporal course of the TMPs for one representative source point on epicardium (50th) over the 60 testing times. The reconstruction TMPs with GA-SVR, DE-SVR, and PSO-SVR methods are all compared with original TMPs.The second experiment reconstructs the TMPs on all the heart surface points at one time instant. These inverse ECG solutions are shown in Figure3, in which two sequential testing time points (27 and 51 ms after the first ventricular breakthrough) are presented to illustrate the performances of the GA-SVR, DE-SVR, and PSO-SVR methods. It can be seen that, among the three parameters optimization methods, the hybrid PSO-SVR performs better than GA-SVR and DE-SVR methods because its reconstructed solutions are much closer to the simulated TMPs distributions.The TMP distribution on the ventricular surface at two sequential testing time points. (a) shows the comparison of reconstructed TMP distributions when the time point is 27 ms and (b) shows the comparison of reconstructed TMP distributions when the time point is 51 ms. Moreover, the upper row shows the TMPs from an anterior view and the lower from a posterior view. In each row, the first figure is the original TMPs, and the rest three are the reconstructed TMPs by using the GA-SVR, DE-SVR, and PSO-SVR methods.
(a)
(b)Based on the simulation information of TMPs, we can evaluate the accuracy of the reconstructed TMPs of these three different intelligent optimization algorithms at the testing time by mean square error (MSE), the relative error (RE), and the correlation coefficient (CC):(12)MSE=1n∑i=1n(φc-φe)2,RE=∥φc-φe∥φe,CC=∑i=1n[(φc)i-φ¯c][(φe)i-φ¯e]∥φc-φ¯e∥∥φe-φ¯e∥.
Here,φeis the simulated TMP distribution at time t, and φc is the reconstructed TMPs. The quantities φ¯e and φ¯c are the mean value ofφeandφcover the whole ventricular surface nodes at time t. n is the number of nodes on the ventricular surface.The MSE, RE, and CC of the three sets of reconstructed TMPs over the 60 testing samples can be found in Figure4. In contrast to GA-SVR and DE-SVR, the PSO-SVR method can yield rather better results with lower RE, MSE, and a higher CC. To validate the robustness of the three methods in reconstructing the TMPs, the quantities m-MSE, m-RE and m-CC, that is, the mean value of MSE, RE, and CC at different testing time instants, are presented. Moreover, std-MSE, std-RE and std-CC, that is, the standard deviation of MSE, RE, and CC over the whole testing samples, are also investigated. As shown in Figure 5, the mean value of the MSE, RE, and CC over the 60 testing samples obtained by GA-SVR, DESVR, and PSO-SVR methods are presented, and the standard deviations of those 60 MSEs, REs and CCs are also provided.The performances of the reconstructed TMPs over 60 sampling times by using GA-SVR, DE-SVR, and PSO-SVR, respectively. (a) The mean square error (MSE) of the reconstructed TMPs; (b) The relative errors (RE) of the reconstructed TMPs; (c) The correlation coefficient (CC) of the reconstructed TMPs.
(a)
(b)
(c)The m-MSE, m-RE, m-CC, std-MSE, std-RE, and std-CC of 60 testing sampling times by using GA-SVR, DE-SVR, and PSO-SVR, respectively, in terms of mean ± standard deviation. (a) m-MSE and std-MSE, (b) m-RE and std-RE, and (c) m-CC and std-CC.
(a)
(b)
(c)


### 4.3. Comparison of Reconstruction Efficiency

In solving the inverse ECG problem, seeking the optimal parameters in the settled range is time consuming. Therefore, how to successfully excogitate a relatively timesaving optimal method matters significantly in the practical application, particularly in the diagnosis of cardiac disease. Here, we figure out the mean time of the three optimization methods in selecting the optimal hyperparameters of SVR. In order to forecast the 478 TMPs over 60 testing samples on an epi- and endocardial surface, 478 sets of parameters (C,σ2, and ε) have to be selected in each optimization method and with the testing data, 478 models are built accordingly. In this study, the mean time for selecting the optimal parameters by adopting GA, DE, and PSO is 208.27 s, 342.53 s, and 130.16 s, respectively. It can be seen that PSO can lead to a convergence more quickly and lessen more time than GA and DE. Therefore, according to the comparison of their efficiencies, PSO-SVR method significantly outperforms the GA-SVR and DE-SVR in the reconstruction of the TMP, which is the most efficient one among these three optimization methods.


## 4.1. The Result of Parameters Selection

In this study, the GA, DE, and PSO algorithms were proposed to seek the corresponding optimal hyperparameter (C,σ2,  and ε) of the 478 SVR models. Table 1 displays the average value and the standard deviations of three parameters in 478 SVR models by using GA, DE, and PSO optimization algorithms. From it, we can find that the gaps among these optimal parameters are obvious by using three different optimization methods. The standard deviations of ε are neglected because they are very small.Table 1
The average and the standard deviations of optimal hyperparameters (C, σ2, and ε) by using GA, DE, and PSO optimization algorithms.
Parameter
Method
GA
DE
PSO
C
1
447.61
±
594.32
7624.14
±
913.67
5.3
1
±
3.02
σ
2
0
.
64
±
0.55
1
.
59
±
1.38
1
.
2
2
±
1.02
ε
1.4
8
e
-
00
5
1.80
e
-
00
5
1.98
e
-
005


## 4.2. Comparison of Reconstruction Accuracy

This study compares the reconstruction performance based on the proposed three intelligent optimization algorithms for selecting the hyperparameters of SVR. By using GA-SVR, DE-SVR, and PSO-SVR hybrid models, two types of experiments for reconstructing the TMPs have been conducted. The first experiment reconstructs the TMPs for one representative source point (the 50th point) on the epicardium over 60 testing times, as is depicted in Figure2. Since the gaps among these three methods are so slight that we have to add a statistical analysis to validate the results. The mean values and standard deviations of these 60 testing errors between the simulated TMPs and the reconstructed TMPs by using GA-SVR, DE-SVR, and PSO-SVR methods are investigated, as shown in the Table 2. From Figure 2 and Table 2, the PSO-SVR outperforms the GA-SVR and DE-SVR methods in reconstructing the TMPs for one representative source point over all the testing times. Furthermore, the DE-SVR is superior to the GA-SVR method in reconstructing the TMPs.Table 2
The mean values and standard deviations of these 60 testing errors between the simulated TMPs and the reconstructed TMPs by using GA-SVR, DE-SVR, and PSO-SVR methods at one representative source point (the 50th point on the epicardium).
Indexes
Method
GA-SVR
DE-SVR
PSO-SVR
Mean value
2.96
2.69
1.55
Standard deviation
1.92
1.67
1.03Figure 2
The temporal course of the TMPs for one representative source point on epicardium (50th) over the 60 testing times. The reconstruction TMPs with GA-SVR, DE-SVR, and PSO-SVR methods are all compared with original TMPs.The second experiment reconstructs the TMPs on all the heart surface points at one time instant. These inverse ECG solutions are shown in Figure3, in which two sequential testing time points (27 and 51 ms after the first ventricular breakthrough) are presented to illustrate the performances of the GA-SVR, DE-SVR, and PSO-SVR methods. It can be seen that, among the three parameters optimization methods, the hybrid PSO-SVR performs better than GA-SVR and DE-SVR methods because its reconstructed solutions are much closer to the simulated TMPs distributions.The TMP distribution on the ventricular surface at two sequential testing time points. (a) shows the comparison of reconstructed TMP distributions when the time point is 27 ms and (b) shows the comparison of reconstructed TMP distributions when the time point is 51 ms. Moreover, the upper row shows the TMPs from an anterior view and the lower from a posterior view. In each row, the first figure is the original TMPs, and the rest three are the reconstructed TMPs by using the GA-SVR, DE-SVR, and PSO-SVR methods.
(a)
(b)Based on the simulation information of TMPs, we can evaluate the accuracy of the reconstructed TMPs of these three different intelligent optimization algorithms at the testing time by mean square error (MSE), the relative error (RE), and the correlation coefficient (CC):(12)MSE=1n∑i=1n(φc-φe)2,RE=∥φc-φe∥φe,CC=∑i=1n[(φc)i-φ¯c][(φe)i-φ¯e]∥φc-φ¯e∥∥φe-φ¯e∥.
Here,φeis the simulated TMP distribution at time t, and φc is the reconstructed TMPs. The quantities φ¯e and φ¯c are the mean value ofφeandφcover the whole ventricular surface nodes at time t. n is the number of nodes on the ventricular surface.The MSE, RE, and CC of the three sets of reconstructed TMPs over the 60 testing samples can be found in Figure4. In contrast to GA-SVR and DE-SVR, the PSO-SVR method can yield rather better results with lower RE, MSE, and a higher CC. To validate the robustness of the three methods in reconstructing the TMPs, the quantities m-MSE, m-RE and m-CC, that is, the mean value of MSE, RE, and CC at different testing time instants, are presented. Moreover, std-MSE, std-RE and std-CC, that is, the standard deviation of MSE, RE, and CC over the whole testing samples, are also investigated. As shown in Figure 5, the mean value of the MSE, RE, and CC over the 60 testing samples obtained by GA-SVR, DESVR, and PSO-SVR methods are presented, and the standard deviations of those 60 MSEs, REs and CCs are also provided.The performances of the reconstructed TMPs over 60 sampling times by using GA-SVR, DE-SVR, and PSO-SVR, respectively. (a) The mean square error (MSE) of the reconstructed TMPs; (b) The relative errors (RE) of the reconstructed TMPs; (c) The correlation coefficient (CC) of the reconstructed TMPs.
(a)
(b)
(c)The m-MSE, m-RE, m-CC, std-MSE, std-RE, and std-CC of 60 testing sampling times by using GA-SVR, DE-SVR, and PSO-SVR, respectively, in terms of mean ± standard deviation. (a) m-MSE and std-MSE, (b) m-RE and std-RE, and (c) m-CC and std-CC.
(a)
(b)
(c)


## 4.3. Comparison of Reconstruction Efficiency

In solving the inverse ECG problem, seeking the optimal parameters in the settled range is time consuming. Therefore, how to successfully excogitate a relatively timesaving optimal method matters significantly in the practical application, particularly in the diagnosis of cardiac disease. Here, we figure out the mean time of the three optimization methods in selecting the optimal hyperparameters of SVR. In order to forecast the 478 TMPs over 60 testing samples on an epi- and endocardial surface, 478 sets of parameters (C,σ2, and ε) have to be selected in each optimization method and with the testing data, 478 models are built accordingly. In this study, the mean time for selecting the optimal parameters by adopting GA, DE, and PSO is 208.27 s, 342.53 s, and 130.16 s, respectively. It can be seen that PSO can lead to a convergence more quickly and lessen more time than GA and DE. Therefore, according to the comparison of their efficiencies, PSO-SVR method significantly outperforms the GA-SVR and DE-SVR in the reconstruction of the TMP, which is the most efficient one among these three optimization methods.


## 5. Conclusion

In the study of the inverse ECG problem, SVR method is a powerful technique to solve the nonlinear regression problem and can serve as a promising tool for performing the inverse reconstruction of the TMPs. This study introduces three optimization methods (GA, DE, and PSO) to determine the hyperparameters of the SVR model and utilizes these models to reconstruct the TMPs from the remote BSPs. Feature extraction using the KPCA is also adopted to preprocess the original input data. The experimental results demonstrate that PSO-SVR is a relatively effective method in terms of accuracy and efficiency. By comparison, PSO-SVR method is superior to GA-SVR and DE-SVR methods, whose reconstructed TMPs are close to the simulated TMPs. Besides, PSO-SVR is much more efficient in determining the optimal parameters and in building the predicting model. According to these results, the PSO-SVR method can serve as a promising tool to solve the nonlinear regression problem of the inverse ECG problem.



---

*Source: 158056-2013-07-25.xml*