A Hybrid Model of Maximum Margin Clustering Method and Support Vector Regression for Noninvasive Electrocardiographic Imaging

Abstract

Noninvasive electrocardiographic imaging, such as the reconstruction of myocardial transmembrane potentials (TMPs) distribution, can provide more detailed and complicated electrophysiological information than the body surface potentials (BSPs). However, the noninvasive reconstruction of the TMPs from BSPs is a typical inverse problem. In this study, this inverse ECG problem is treated as a regression problem with multi-inputs (BSPs) and multioutputs (TMPs), which will be solved by the Maximum Margin Clustering- (MMC-) Support Vector Regression (SVR) method. First, the MMC approach is adopted to cluster the training samples (a series of time instant BSPs), and the individual SVR model for each cluster is then constructed. For each testing sample, we find its matched cluster and then use the corresponding SVR model to reconstruct the TMPs. Using testing samples, it is found that the reconstructed TMPs results with the MMC-SVR method are more accurate than those of the single SVR method. In addition to the improved accuracy in solving the inverse ECG problem, the MMC-SVR method divides the training samples into clusters of small sample sizes, which can enhance the computation efficiency of training the SVR model.

1. Introduction

The technique of noninvasive imaging of the heart’s electrical activity from the body surface potentials (BSPs) constitutes one form of the inverse problem of ECG [1, 2]. Approaches to solving the inverse ECG problem have been usually based on either an activation-based model or a potential-based model, which includes epicardial, endocardial, or transmembrane potentials. Activation-based models are used to investigate the arrival time of the propagation wavefront within the myocardium [3, 4]. The potential-based models are used to evaluate the potential values on the cardiac surface [5–7] or within the myocardium [8] at certain time instants. In this study, we explore a new solution for ECG inverse problem using the potential-based approach.

Due to its inherent ill-posed property, the inverse ECG problem is usually solved by “regularization” techniques. In the last decades, numerous regularization methods have been proposed to solve this ill-posed problem, including truncated total least squares (TTLS) [9], GMRes [10], and the LSQR [11, 12]. Most of them are essentially L2-norm based regularization schemes, which inherently lead to considerable smoothness of the inverse solutions. L1-norm regularization method can overcome this drawback of L2-norm regularization method, which has been applied for epicardial potential reconstruction [13–15]. Although the above-mentioned regularization methods can more or less deal with the geometry and measurement noises for the ECG inverse problems, which depends on the regularization parameters, the robustness of the inverse solution is not always guaranteed. In this paper, without seeking assistance from the regularization techniques, we explore an alternative, more robust approach to solve the inverse ECG problem. The method is called Support Vector Regression (SVR) [16]. To find the solution for the inverse ECG problem, a regression model will be set up with multi-inputs (BSPs) and multioutputs (transmembrane potentials, TMPs). This statistic method based solution will be assessed with the quality of the inversely predicated TMPs from the measured BSPs. Compared with conventional regularization methods (e.g., zero order Tikhonov and LSQR), the SVR method can produce more accurate results in terms of reconstruction of the transmembrane potential distributions on epi- and endocardial surface. In addition, when the PCA and KPCA are adopted to extract useful features from the original inputs for building the SVR model, the SVR method with feature extraction (PCA-SVR and KPCA-SVR) outperforms that without the extract feature extraction (single SVR) in terms of the reconstruction of the TMPs [17].

Compared with using single SVR model, the hybrid models by integrating difference methods show better performance. The self-organizing map (SOM) is an unsupervised and competitive learning algorithm, which can be viewed as clustering techniques [18]. Combining the SOM with SVR or LS-SVM, the proposed hybrid method has the potential to find better inverse solutions than using a single SVR model [19, 20]. Xu et al. [21] proposed the Maximum Margin Clustering (MMC) method, which performs clustering by simultaneously finding the large margin separating hyperplane between clusters. The MMC method has been successfully applied to many clustering problems [22]. However, its efficiency is an issue of concern. Recently, Zhang et al. proposed [23, 24] an efficient approach for solving the MMC via an alternative optimization procedure, which was implemented by using the SVR method with the Laplacian loss in the inner optimization subproblem. The modified MMC algorithm is more accurate, much faster and therefore more practical for solving engineering inverse problems. In this paper, the hybrid model of modified MMC method and SVR is proposed to solve the inverse ECG problem, which is referred to as an MCC-SVR method. The conference version of this submission has appeared in CINC 2011 [25]. This submission has undergone substantial revisions and offers extended experiment results.

The main purpose of this study is to use an MCC-SVR model to investigate the reconstruction capability of TMPs. In this study, based on our previously developed realistic heart-torso model, the equivalent double layer (EDL) source model method was applied to generate the data set for training and testing the SVR model. The proposed algorithm was also compared with a single SVR model for noninvasive ECG imaging.

2. Theory and Methodology

The framework of the proposed MCC-SVR method is shown in Figure 1. The MCC method is used to classify the input data; the SVR is then applied to construct the regression model of each cluster.

Figure 1 

The framework of the proposed MCC-SVR method.

2.1. Maximum Margin Clustering (MMC) Method [23, 24]

The clustering principle is to find a labeling to identify dominant structures in the data and to group similar instances together, so the margin obtained would be maximal over all possible labelings, that is, given a training set {(𝑥_𝑖,𝑦_𝑖)}^𝑛_𝑖=1, where 𝑥_𝑖∈𝜒 is the input and 𝑦_𝑖∈{±1} is the output. The SVM finds a large margin hyperplane to separate patterns of opposite classes by the classify function 𝑓(𝑥) [26]: 𝑓(𝑥)=𝜔^𝑇𝜑(𝑥)+𝑏, (1) where 𝜑(𝑥) denotes the high-dimensional feature space, which is nonlinearly mapped from the input space 𝑥 by the kernel function 𝑘, 𝜔 is the normal vector of the hyperplane, and 𝑏 is the offset of the hyperplane. Computationally, this leads to the following optimization problem [24, 26]: min_{𝜔,𝑏,𝜉}‖𝜔‖²+2𝐶𝜉^𝑇𝑒𝑦subjectto^𝑖𝜔𝜑𝑥_𝑖+𝑏≥1−𝜉_𝑖𝜉_𝑖≥0,𝑖=1,…,𝑛, (2) where 𝜉=[𝜉₁,…,𝜉_𝑛]^𝑇 is the vector of a slack variable for the errors, and 𝐶>0 is the trade-off parameter between the smoothness ‖𝜔‖² and the fitness (𝜉^𝑇𝑒) of the decision function 𝑓(𝑥).

MMC attempts to extend large margin methods to allocate the input data points to different classes, leading to large separation between the different classes. Here, the case with two clusters is considered in this work. Since one could simply assign all the data points to the same class and obtain an unbounded margin, a proper constraint on the class balance needs to be imposed. Xu et al. [21] introduced a class constraint that requires 𝑦 to satisfy−ℓ≤𝑒^𝑇𝑦≤ℓ, (3) where ℓ≥0 is a user-defined constant controlling the class imbalance. Then the margin is maximized with respect to both unknown 𝑦 and unknown SVM parameter (𝜔, 𝑏) as follows: min_𝑦min_{𝜔,𝑏,𝜉}‖𝜔‖²+2𝐶𝜉^𝑇𝑒𝑦subjectto^𝑖𝜔𝜑𝑥_𝑖+𝑏≥1−𝜉_𝑖𝜉_𝑖≥0,𝑦_𝑖∈{±1},𝑖=1,…,𝑛−ℓ≤𝑒^𝑇𝑦≤ℓ. (4) The origin nonconvex MMC problem in (4) can be formulated as a sequence of QPs which can be solved by some efficient QP solvers. However, it suffers from a premature convergence and easily gets stuck in poor local optima. Zhang et al. [23, 24] proposed to replace the SVM by SVR with Laplacian loss, which can lead to a significant improvement in the clustering performance compared to that of iterative SVM procedure. The primal problem of SVR with Laplacian loss can be formulated as min_{𝜔,𝑏,𝜉𝑖,𝜉∗𝑖}‖𝜔‖²+2𝐶^𝑛_𝑖=1𝜉_𝑖+𝜉^∗_𝑖𝑦subjectto^𝑖−𝜔^𝑇𝜑𝑥_𝑖+𝑏≤𝜉_𝑖𝜔^𝑇𝜑𝑥_𝑖+𝑏−𝑦_𝑖≤𝜉^∗_𝑖𝜉for𝑖=1,…,𝑛,_𝑖≥0,𝜉^∗_𝑖≥0, (5) where 𝜉_𝑖 and 𝜉^∗_𝑖 are slack variables. With the obtained labels, the MMC problem based on the iterative SVR with the Laplacian loss becomes min_{𝜔,𝑏,𝜉𝑖,𝜉∗𝑖}‖𝜔‖²+2𝐶^𝑛_𝑖=1𝜉_𝑖+𝜉^∗_𝑖𝑦subjectto^𝑖−𝜔^𝑇𝜑𝑥_𝑖+𝑏≤𝜉_𝑖𝜔^𝑇𝜑𝑥_𝑖+𝑏−𝑦_𝑖≤𝜉^∗_𝑖𝜉_𝑖≥0,𝜉^∗_𝑖𝑦≥0for𝑖=1,…,𝑛_𝑖∈{±1}−ℓ≤𝑒^𝑇𝑦≤ℓ. (6) After 𝜔 is obtained from the optimization of SVR, the problem in (6) is reduced to the form min^𝑛_{𝑦,𝑏𝑖=1}𝜔^𝑇𝜑𝑥_𝑖+𝑏−𝑦_𝑖𝑦subjectto^𝑖∈{±1},𝑖=1,…,𝑛−ℓ≤𝑒^𝑇𝑦≤ℓ. (7) According to Zhang’s proposition [24], for a fixed 𝑏, the optimal strategy to determine the 𝑦_𝑖’s in (7) is to assign all 𝑦_𝑖’s as −1 for those with 𝜔^𝑇𝜑(𝑥_𝑖)+𝑏<0 and assign 𝑦_𝑖’s as 1 for those with 𝜔^𝑇𝜑(𝑥_𝑖)+𝑏>0. The bias 𝑏 can be determined as follows. (i) we sort the 𝜔^𝑇𝜑(𝑥_𝑖)’s and use the set of midpoints between any two consecutive sorted values as the candidates of 𝑏; (ii) from these sorted 𝑏’s, the first and the last (𝑛−ℓ)/2 of them can be dropped, and the middle ℓ can be remained; (iii) for each remaining candidate, we determine the 𝑦_𝑖’s according to the above proposition and compute the corresponding objective value in (7); (iv) finally, we choose the 𝑏 that has the smallest objective. The complete iterative SVR procedure for MCC method is shown in Algorithm 1.

(1) Initialize the labels 𝑦 by using the 𝑘-means; (2) Fix 𝑦, and perform SVR with Laplacian loss; (3) Compute 𝜔 from the Karush-Kuhn-Tucker (KKT) conditions; (4) Compute the bias 𝑏 as described above; (5) Assign the labels as 𝑦𝑖=sign(𝜔𝑇𝜑(𝑥𝑖)+𝑏); (6) Repeat steps 2–5 until convergence.

Algorithm 1 

Iterative SVR procedure for MCC method.

2.2. Support Vector Regression (SVR) Model

The SVR algorithm [26] is only briefly described here; for details, see [16, 26]. As a linear regression model, the SVR algorithm relies on an estimation of a linear regression function: 𝑓(𝑥)=⟨𝜔,𝑥⟩+𝑏,(𝜔,𝑥∈ℜ), (8) where 𝜔 and 𝑏 are the slope and offset of the regression linear, and ⟨⋅,⋅⟩ denotes the dot product in ℜ. The above regression problem can be written as a convex optimization problem: 1min2‖𝜔‖²𝑦subjectto^𝑖−𝜔,𝑥_𝑖−𝑏≤𝜀𝜔,𝑥_𝑖+𝑏−𝑦_𝑖≤𝜀. (9) In (9), an implicit assumption is that a function 𝑓 essentially approximates all pairs (𝑥_𝑖, 𝑦_𝑖) with 𝜀 precision, but sometimes this may not be the case. Therefore, one can introduce two additional positive slack variables 𝜉_𝑖, 𝜉^∗_𝑖 to refine the estimation of variables 𝜔 and 𝑏. Now (9) can be reformulated [16] as 1min2‖𝜔‖²+𝐶^𝑛_𝑖=1𝜉_𝑖+𝜉^∗_𝑖𝑦subjectto^𝑖−𝜔,𝑥_𝑖−𝑏≤𝜀+𝜉_𝑖𝜔,𝑥_𝑖+𝑏−𝑦_𝑖≤𝜀+𝜉^∗_𝑖𝜉_𝑖,𝜉^∗_𝑖≥0, (10) where the constant 𝐶 is a trade-off parameter and 𝑛 denotes the number of samples; 𝜉_𝑖 represents the upper training error, and 𝜉^∗_𝑖 is the lower training error subject to 𝜀 intensive tube. According to the strategy outlined by Vapnik [26], using Lagrange multipliers, the constrained optimization problem shown in (3) can be further restated as the following equation: 𝑓𝑥,𝛼_𝑖,𝛼^∗_𝑖=^𝑛_𝑖=1𝛼_𝑖−𝛼^∗_𝑖𝐾𝑥_𝑖,𝑥+𝑏subjectto^𝑛_𝑖=1𝛼_𝑖−𝛼^∗_𝑖=0,0≤𝛼_𝑖,𝛼^∗_𝑖≤𝐶, (11) where 𝛼_𝑖 and 𝛼^∗_𝑖 are the Lagrange multipliers. The term 𝐾(𝑥_𝑖,𝑥_𝑗) in (11) is defined as the kernel function, whose values are the inner product of two vectors 𝑥_𝑖 and 𝑥_𝑗 in the feature space 𝜑(𝑥_𝑖) and 𝜑(𝑥_𝑗). And bias 𝑏 can be computed as follows: 𝑦𝑏=^𝑖−∑^𝑛_𝑗=1𝛼_𝑖−𝛼^∗_𝑖𝐾𝑥_𝑗,𝑥_𝑖−𝜀for𝛼_𝑖∈()𝑦0,𝐶_𝑖−∑^𝑛_𝑗=1𝛼_𝑖−𝛼^∗_𝑖𝐾𝑥_𝑗,𝑥_𝑖+𝜀for𝛼^∗_𝑖∈(0,𝐶). (12) The kernel function handles any dimension feature space with no explicit calculation of 𝜑(𝑥). In this study, the Gaussian kernel function is chosen as the SVR’s application mapping in this study: 𝐾𝑥_𝑖,𝑥_𝑗−𝑥=exp^𝑖−𝑥_𝑗²2𝜎², (13) where 𝑥_𝑖 and 𝑥_𝑗 are input vector spaces; 𝜎² is the bandwidth of the kernel function.

In this study, an accurate and fast approach based on the GA and the simplex search techniques is presented to determine the optimal hyperparameters of the SVR model [17], as shown in Figure 2. The GA algorithm used here is based on a GA toolbox developed by Chipperfield et al. [27], and the simplex optimization method is implemented using the MATLAB optimization toolbox. The developed SVR model was trained and validated with the software LIBSVM [28].

Figure 2 

GA-Simplex optimization procedure for the parameter selection in the MCC-SVR model.

2.3. Simulation Protocol and Data Set

The SVR model is tested with our previously developed realistic heart-torso model [6, 17]. In this study, an equivalent double layer (EDL) source model is adopted to simulate the cardiac equivalent source, which represents the cardiac electrical activity by means of double layer source on the closed surface (including the endo- and epicardial surface of ventricle). For the ECG inverse problem studies, the ventricular surface TMPs and body surface potentials (BSPs) are evaluated based on the EDL source model. The transfer matrix 𝐴 between TMPs and BSPs is evaluated by the boundary element method (BEM), and it has the dimension of 412×478 and its condition number (the ratio of largest and smallest singular values) is 5.6 × 10¹². As shown in Figure 3, the EDL source method is used to obtain the BSPs 𝜑_𝐵 and the TMPs 𝜑𝑚. For the construction of the training and testing data set Different Action Potentials (APs) for various myocardial cells and the normal Ventricular Excitation Sequence (VES) are used to calculate the TMPs (𝜑𝑚) at different times; from the calculated TMPs, the corresponding BSPs are deduced with the transfer matrix 𝐴.

Figure 3 

The block diagram of the simulation protocol for the construction of the data sets.

In this study, a normal ventricular excitation data set is prepared for the setup of 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 𝜑_𝐵 and TMPs 𝜑𝑚 temporal data sets are numerically recorded; in addition, the 30 dB simulated Gaussian white noise is added into the BSPs 𝜑_𝐵 representing the measurement noises. 60 datasets at times of 3 ms, 9 ms, 15 ms, …, and 357 ms after the first ventricular breakthrough are used as testing samples to evaluate the generalization capacity of the proposed SVR model. The rest 298 in 358 data sets are employed as the training samples for building the SVR model. With the consideration of a wide numerical range of the 𝜑_𝐵 values, for each time, the 𝜑_𝐵 values can be scaled to the range (0, 1): 𝜑_𝐵𝑡𝑁=𝜑_𝐵𝑡−𝜑_𝐵𝑡min𝜑_𝐵𝑡max−𝜑_𝐵𝑡min, (14) where 𝜑_𝐵𝑡 are the body surface potentials at time instant 𝑡, 𝜑_𝐵𝑡max is the maximum value of BSPs at the time 𝑡, and 𝜑_𝐵𝑡min is the minimum value of BSPs at the time 𝑡.

As the TMPs are known in advance in the simulation study, the accuracy of reconstructed TMPs at the testing time 𝑡 can be evaluated by either relative errors (REs): 𝜑RE=^𝑐𝑡−𝜑^𝑒_𝑡𝜑^𝑒_𝑡, (15) or the correlation coefficient (CC), given by ∑CC=^𝑛𝑖=1𝜑^𝑐_𝑡𝑖−𝜑^𝑐_𝑡𝜑^𝑒_𝑡𝑖−𝜑^𝑒_𝑡𝜑^𝑐_𝑡−𝜑^𝑐_𝑡𝜑^𝑒_𝑡−𝜑^𝑒_𝑡, (16) where 𝑛 is the number of nodes on the ventricular surface. 𝜑^𝑒_𝑡 denotes the simulated TMPs distribution at time 𝑡, and 𝜑^𝑐_𝑡 are inversely computed. The quantities 𝜑^𝑐_𝑡 and 𝜑^𝑒_𝑡 are the mean value of 𝜑^𝑐_𝑡 and 𝜑^𝑒_𝑡 over the whole ventricular surface nodes at time 𝑡.

3. Results

According to the MCC method, the above 298 training samples are classified four clusters as shown in Figure 4(a), and the numbers of the four clusters is 80, 74, 70, and 74, respectively. Then the individual SVR model is trained for each cluster, and the hyperparameters are determined using the GA-Simplex method. For 60 testing samples, the MCC method is used to find their corresponding clusters, as shown in Figure 4(b).

(a)

(b)

Figure 4 

On one representative epicardial point, (a) the 298 training samples are classified into four clusters by using MCC method; (b) for the 60 testing samples, the MCC method is used to find their corresponding clusters.

To illustrate the performances of the reconstructed TMPs, four sequential testing time points (3, 15, 27, and 39 ms after ventricle excitation) are presented. The inverse ECG solutions are shown in Figure 5; in contrast to the conventional regularization methods, such as zero order Tikhonov regularization method and LSQR regularization method, the single SVR method can yield rather better results with lower RE and higher CC. Moreover, it can be seen that the MCC-SVR method offers superior performances than the single SVR method, as its solution is more close to the simulated TMPs distributions. The time courses of the simulated TMPs and reconstructions for one representative source point on the heart surface are depicted in Figure 6. It can be found that, in reconstructing the TMPs for one representative source point over all the testing times, the MCC-SVR method offers better solution compared with single SVR method.

(a)

(b)

(c)

(d)

(e)

Figure 5 

The TMPs distribution on the ventricular surface at four sequential testing time points (3, 15, 27, and 39 ms, respectively, after the first ventricular breakthrough). In each subfigure, the upper row shows the TMPs distribution from an anterior view and the lower from a posterior view. (a) The simulated TMPs by using the equivalent double layer (EDL) source model; (b) the reconstructed TMPs by using zero order Tikhonov method; (c) the reconstructed TMPs by using the LSQR method; (d) the reconstructed TMPS by using the single SVR method; (e) the reconstructed TMPs by using the MCC-SVR method.

Figure 6 

The time courses of the TMPs for one representative source point on the heart surface. The reconstruction TMPs over the 60 testing times with SVR method and the MCC-SVR method are all compared with those simulated TMPs.

The RE and CC of the reconstructed TMPs with different methods can be found in Figure 7. In contrast to the single SVR method, the MCC-SVR method can yield improved results with a lower RE and a higher CC over the 60 testing samples.

(a)

(b)

Figure 7 

The performances of the reconstructed TMPs over 60 sampling times by using SVR and MCC-SVR, respectively. (a) The REs of the reconstructed TMPs; (b) the CCs of the reconstructed TMPs.

By dividing the training samples into smaller clusters, the training time for each cluster can be reduced. The training times for each cluster by using the MCC-SVR method are 6715.9 seconds, 6821.1 seconds, 4550.6 seconds, and 5162.4 seconds, respectively, and the total time of the four clusters is 23250 seconds. When using the single SVR model to train the model of the all training samples, it takes 35233.4 seconds.

4. Discussion and Conclusion

In this study, MCC-SVR method is proposed to solve the noninvasive ECG imaging problem. Here, the MMC approach is adopted to cluster the training samples firstly, and then SVR method is applied to construct the model for each cluster. After building different cluster models, for the testing sample, we can find its matched cluster and then use the corresponding SVR model to reconstruct the TMPs. From the reconstructed TMPs as shown in Figures 5 and 6, it can be seen that the MCC-SVR methods offer better solution compared with single SVR method. According to the evaluation indices RE and CC, the performances of the reconstructed TMPs by using the MCC-SVR can constantly converge to a smaller RE and a higher CC on the testing samples than those of the single SVR method, as shown in Figure 7. In terms of the computation efficiency of the training SVR model, for the given training samples, the MCC-SVR method can save about 34% time than the single SVR method. With the increasing of the training samplings, the MCC-SVR method lessens more training time than the single SVR method. Moreover, the training process for each cluster can be implemented simultaneously using parallel computing, therefore further enhance the training efficiency.

In summary, this paper proposed the MCC-SVR method for the inverse solutions of the ECG problem. The new algorithm was tested and compared with single SVR schemes using a realistic heart-torso model. The experimental results show that the MCC-SVR can improve the generalization performance of the single SVR in reconstructing the TMPs, leading to a more accurate reconstruction of the TMPs. In our future work, we plan to improve the MCC-SVR method for solving various nonlinear regression problems in noninvasive ECG imaging.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (30900322, 61070063, and 61101046). This work is also funded by the Zhejiang Provincial Natural Science Foundation of China (Y2090398).