Coupling tensor decomposition algorithm suitable for multi-group high-dimension EEG data analysis

Abstract

The invention discloses a coupling tensor decomposition algorithm suitable for multi-group high-dimension EEG (electroencephalogram) data analysis. Constructing a multi-dimensional tensor joint decomposition model by adopting a CP model; applying linear coupling constraint, and effectively extracting common components and specific components of the data; performing component extraction by adopting an alternating direction multiplier method optimization algorithm; non-negative constraint is applied in the algorithm updating process, and components and components with more representative significance are further extracted. According to the method, while the original high-dimensional structure of the data is kept, the implicit features in different groups of EEG data and the interaction relationship between the data are fully mined, so that the accuracy and effectiveness of EEG data analysis are improved; and a technical basis can be provided for joint processing and analysis of multi-set brain imaging data in the fields of brain science research, brain disease diagnosis, brain-computer interfaces and the like.

Description

Technical Field

[0001] The present invention belongs to the technical field of biomedical signal processing, and specifically relates to a coupled tensor decomposition algorithm suitable for multi-group high-dimensional EEG data analysis. Background Art

[0002] As a relatively objective and effective electrophysiological indicator for observing brain activity, EEG has important practical value in the current diagnosis of neurological diseases and brain science research. Through the analysis of EEG data, we need to explore the high-dimensional characteristics of the data under different subjects, different groups, different time and space scales. At the same time, there is common information and difference information between different subjects and different groups. However, in most EEG studies, people do not consider the intrinsic connection between data and the inherent characteristics within the data at the same time. More and more studies have shown that this information is helpful to reveal the potential internal mechanism between brain function and various cognitive tasks or brain diseases. Therefore, the joint analysis of multi-group EEG data will contribute to the further development of brain science research.

[0003] Coupled tensor decomposition is a method for the joint decomposition of multiple tensor data, which can solve the problem of multiple couplings between components in multi-group EEG data. Through coupled tensor decomposition, under the premise of ensuring the high-dimensional multilinear structural characteristics of multi-group EEG data in time, space and frequency, the common information and difference information of each group in dimensions such as time, space and frequency can be extracted. At present, coupled tensor decomposition technology has become a new research hotspot in the field of signal processing, and has shown greater advantages than tensor decomposition in practical applications such as joint blind source separation, array signal processing and hyperspectral image fusion. However, for the current EEG data with high-dimensional characteristics such as multiple subjects, time, space and groups, the existing coupled tensor decomposition models and methods have problems such as slow convergence speed, low optimization accuracy and inflexible coupling / constraint conditions. At the same time, most algorithms are only applicable to situations with high signal-to-noise ratios, and there are still certain limitations in performance when processing real EEG data with high noise and large model errors. Summary of the invention

[0004] The purpose of the present invention is to provide a coupling tensor decomposition algorithm suitable for the analysis of multi-group high-dimensional EEG data. By introducing non-negativity constraints, flexible linear coupling constraints and the ADMM algorithm based on CP decomposition, many problems existing in the above-mentioned existing coupling tensor decomposition process are solved, and results that are significantly better than existing algorithms are achieved in the joint analysis of multi-group high-dimensional EEG data.

[0005] The technical solution of the present invention is to use the CANDECOMP / PARAFAC (CP) model to construct a multidimensional tensor joint decomposition model; impose linear coupling constraints to extract common components and specific components of the data; use an optimization algorithm based on the alternating direction multiplier method (ADMM) to extract components; impose non-negativity constraints during the algorithm update process to further extract more representative components and components. The specific implementation includes the following steps:

[0006] (1) Input multiple groups of high-dimensional non-negative EEG tensor data s=1,2,…,S; S represents the number of groups, N represents the number of factor matrices or the order of tensors, I n represents the number of components of the corresponding factor matrix;

[0007] (2) Based on the data features, a multidimensional tensor joint decomposition model based on CP decomposition is constructed, which can be expressed as:

[0008]

[0009] in Represents the sth tensor data modulo-n factor matrix The rth column of Represents the sth core tensor, and the non-zero elements r=1,2,…,R (s) Located on its superdiagonal; each factor matrix A (n,s) It consists of two parts: and 0≤L n ≤R (s) All tensor data is shared represents the coupling information of the data, and Corresponding to the private information of each tensor data;

[0010] (3) Based on the above model, the cost function is constructed to obtain the tensor data x (s) Extract the factor matrix A from (n,s) and the core tensor The error between the original tensor and the estimated tensor is minimized using an optimization strategy based on minimization of the Euclidean distance, so the cost function can be expressed as follows:

[0011]

[0012] The alternating optimization strategy is used to solve the cost function, that is, to update a certain factor matrix and core tensor alternately while keeping other factor matrices and core tensors unchanged. For each convex optimization subproblem, the alternating direction method of multipliers (ADMM) algorithm is used to solve it. According to the imposed coupling constraints, the non-coupled factor matrices are updated independently, and the coupled factor matrices are jointly updated. (n,s) and the core tensor The update is explained:

[0013] (4) Calculate the factor matrix A (n,s) , by introducing auxiliary variables The minimization problem of the cost function (1) is transformed into the optimization problem of the augmented Lagrangian function, which can be expressed as follows:

[0014]

[0015] Among them, Λ (n,s) represents the Lagrangian factor (dual variable), the penalty factor ρ (n,s ) is defined as Therefore, the original variable A can be (n,s) , auxiliary variables and the dual variable Λ (n,s) To update, the update method is as follows:

[0016] (5) For the original variable A without coupling components (n,s) , by mathematically calculating and derivatizing equation (2), we can get A (n ,s) The update formula is as follows:

[0017]

[0018] For the original variable A with coupled components (n,s) ,need and The two parts are updated separately, and through mathematical calculation, let The update formula is as follows:

[0019]

[0020] (6) For auxiliary variables without coupling components By taking the derivative of formula (2), we can get the update method as follows: “+” is a truncation operation, and elements less than 0 are set to 0 to impose a non-negative constraint. Similarly, for the original variable auxiliary variable with coupling components The coupled variables and private variables need to be updated separately:

[0021]

[0022] (7) Update the dual variable Λ (n,s) , which is updated as

[0023] (8) Calculate the core tensor Similarly, by introducing the auxiliary variable and the dual variable Γ (s) , the minimization problem of the cost function can be transformed into the optimization problem of the augmented Lagrangian function, which is expressed as follows:

[0024]

[0025] For formula (8), and Γ (s) Derivation, through mathematical calculations, the corresponding solution can be obtained;

[0026] (9) Update the original variables Assume d (s) for The vector of superdiagonal elements of , then d (s) The update method is as follows:

[0027]

[0028] Among them, "vec" represents the vectorization operation of the tensor, and "diag" represents taking the superdiagonal elements of the tensor to form a vector;

[0029] (10) Update auxiliary variables and the dual variable Γ (s) , define "ddiag" as the inverse operation of "diag", and its update method is as follows:

[0030]

[0031] (11) So far, all variable updates have been given. Next, define the iteration termination condition; let the number of iterations be k, and calculate the change in the relative error between two adjacent iterations. when The absolute value of Or k reaches the maximum number of iterations k max When k is reached, the algorithm stops iterating, otherwise it returns to step 6. In this algorithm, the maximum number of iterations is defined as k max =1000, ε is set to 10 -8 ; RelErr is defined as where x (s) and are the original value and estimated value of the tensor data respectively;

[0032] (12) Output all factor matrices A (n,s) and the core tensor

[0033] Beneficial effects of the present invention: Compared with the existing coupled tensor decomposition algorithms based on alternating least squares (ALS), hierarchical alternating least squares (fHALS) and multiplicative updating (MU), the method of the present invention shows obvious advantages in decomposition accuracy on computer simulation data with a signal-to-noise ratio of -10dB-10dB, wherein the factor matrix decomposition accuracy can be improved by 10% to 30%, and the operation time is only half of that based on the MU algorithm. Among them, in the analysis of continuous EEG data collected from 39 subjects including 20 patients with severe depression and 19 matched healthy controls, compared with the traditional independent component analysis (ICA) method and tensor decomposition (TD) method, the present invention can directly extract the common characteristic patterns of brain activity of the two groups of subjects, as well as the independent characteristic patterns corresponding to the healthy control group and the patient group. The results showed that the brain activity states triggered by the patient group and the healthy control group when listening to music were partially similar, but at the same time, the patient group also brought about some unique changes in the characteristics of the brain oscillation network when listening to music. Therefore, the present invention can effectively extract hidden information from multiple groups of EEG data and identify specific changes in brain activity in specific groups. These changes may provide a basis for the clinical diagnosis and treatment of patients with depression. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 This is a workflow diagram for analyzing multiple groups of high-dimensional EEG data in the present invention. DETAILED DESCRIPTION

[0035] In order to facilitate the understanding of the present invention, the present invention will be described more fully below. However, the present invention can be implemented in many different forms and is not limited to the embodiments described herein. On the contrary, the purpose of providing these embodiments is to make the understanding of the disclosure of the present invention more thorough and comprehensive.

[0036] There are continuous EEG data from 39 subjects in total (20 patients with major depression and 19 healthy control subjects) from two groups under natural music stimulation, that is, S=2. The EEG acquisition device has 64 electrodes, and the acquisition data length is 512 seconds. Through preprocessing operations such as short-time Fourier transform, the dimension of each subject's data is 64 (number of channels)*130 (frequency points)*500 (time points). Therefore, for the healthy control group and the depression group, the EEG data is reconstructed into two fourth-order tensors with dimensions of 64 (number of channels)*130 (frequency points)*500 (time points)*19 / 20 (number of subjects). The steps of using the present invention to jointly analyze high-dimensional non-negative EEG data of two groups are as follows: Figure 1 shown.

[0037] Step 1: Input two groups of high-dimensional non-negative EEG tensor data (healthy control group) and (Depression patient group).

[0038] Step 2: Calculate the number of components. The principal component analysis (PCA) method based on interpretable variance is used to calculate the number of components of each tensor. That is, the number of principal components corresponding to the cumulative interpretable variance accounting for 99% is defined as the number of components of the tensor. The calculation yields R (1) =26, R (2) =36. At the same time, correlation analysis is used to determine the number of coupled components between tensor data, and it can be calculated that L1=L2=L3=17, L4=0.

[0039] Step 3: Randomly initialize the tensor factor matrix (for simplicity, the core tensor is assumed to be a unit tensor), and the elements of the factor matrix obey uniform distribution. (1) , the factor matrix is ​​initialized to V (1) , V (2) , V (3) and V (4) ; for x (2) , the factor matrix is ​​initialized to U (1) , U (2) , U (3) and U (4) .

[0040] Step 4: Calculate the estimated tensor using back projection and The RelErr formula is used to calculate the initial tensor relative error RelErr0.

[0041] Step 5: Update the factor matrix. Use formulas (3), (4), (5), (6) and (7) to update the factor matrix V in sequence. (n) and U(n) , by updating one of them and fixing the others, and applying non-negativity constraints. The formula (3) is as follows:

[0042]

[0043] The formula (4) is as follows:

[0044]

[0045] The formula (5) is as follows:

[0046]

[0047] The formula (6) is as follows:

[0048]

[0049] The formula (7) is as follows:

[0050]

[0051] Step 6: Use the RelErr formula to calculate the relative error RelErr of the kth iteration k .

[0052] Step 7: Iteration stop condition determination. Set the minimum threshold to ε = 10 -8 , the maximum number of iterations is k max = 1000. Calculate the change in relative error between two adjacent iterations when The absolute value of Or k reaches the maximum number of iterations k max , the algorithm meets the iteration stop condition, otherwise it returns to step 5.

[0053] Step 8: Extract music features and perform correlation analysis with the final time domain waveform to determine the brain activity stimulated by music. At the same time, extract shared and private features from the final factor matrix, and use the EEGLAB toolbox to outline the shared and private brain topography, power spectrum, time domain waveform and other information of the two groups of subjects.

[0054] The above-mentioned embodiments only express several implementation methods of the present invention, and the descriptions thereof are relatively specific and detailed, but they cannot be understood as limiting the scope of the invention patent. It should be pointed out that, for ordinary technicians in this field, several variations and improvements can be made without departing from the concept of the present invention, and these all belong to the protection scope of the present invention. Therefore, the protection scope of the patent of the present invention shall be subject to the attached claims.

Claims

1. A coupled tensor decomposition algorithm suitable for multi-group high-dimensional EEG data analysis, characterized in that: Use the following steps: (1) Input multiple groups of high-dimensional non-negative EEG tensor data: s=1,2,…,S; where S represents the number of groups, N represents the number of factor matrices or the order of tensors, and I n represents the number of components of the corresponding factor matrix; (2) Construct a multidimensional tensor joint decomposition model based on CP decomposition based on data features; (3) Based on the decomposition model described in step (2), a cost function is constructed to obtain the tensor data Extract the factor matrix A from (n,s) and the core tensor Use an optimization strategy based on Euclidean distance minimization to minimize the error between the original tensor and the estimated tensor; (4) Calculate the factor matrix A (n,s) , by introducing auxiliary variables The minimization problem of the cost function is transformed into the optimization problem of the augmented Lagrangian function; (5) For the original variable A without coupling components (n,s) , by mathematically calculating the derivative of the cost function, we get A (n,s) Update formula (4); for the original variable A with coupling components (n,s) ,need and The two parts are updated separately, and through mathematical calculation, let Formulas (5) and (6) can be updated; (6) For auxiliary variables without coupling components By taking the derivative of the cost function, the update method is: "+" is a truncation operation, and elements less than 0 are set to 0 to impose a non-negative constraint; for the original variable with coupling components, the auxiliary variable Update the coupled variables and private variables separately, and get the update formula (7); (7) Update the dual variable Λ (n,s) , which is updated as (8) Calculate the core tensor By introducing auxiliary variables and the dual variable Γ (s) , the minimization problem of the cost function is transformed into the optimization problem of the augmented Lagrangian function; (9) Update the original variables Assume d (s) for The vector formed by the superdiagonal elements of , we get d (s) Update formula (9); (10) Update auxiliary variables and the dual variable Γ (s) , define "ddiag" as the inverse operation of "diag", and update formula (10); (11) Let the number of iterations be k, and calculate the change in the relative error between two adjacent iterations l k =RelErr k -RelErr k-1 , when l k The absolute value of |l k |<ε or k reaches the maximum number of iterations k max When , the algorithm stops iterating, otherwise it returns to step (5); the maximum number of iterations is defined as k max =1000, ε is set to 10 -8 ; RelErr is defined as in and are the original value and estimated value of the tensor data respectively; (12) Output all factor matrices A (n,s) and the core tensor 2. The coupled tensor decomposition algorithm for multi-group high-dimensional EEG data analysis according to claim 1, characterized in that: The decomposition model in step (2) is expressed as: in Represents the sth tensor data modulo-n factor matrix The rth column of Represents the sth core tensor, and the non-zero elements r=1,2,…,R (s) Located on its superdiagonal; each factor matrix A (n ,s) It consists of two parts: and 0≤L n ≤R (s) ; All tensor data are shared represents the coupling information of the data, and Corresponding to the private information of each tensor data.

3. The coupled tensor decomposition algorithm for multi-group high-dimensional EEG data analysis according to claim 1, characterized in that: The cost function in step (3) is expressed as follows:

4. The coupled tensor decomposition algorithm for multi-group high-dimensional EEG data analysis according to claim 1, characterized in that: The Lagrangian function in step (4) is expressed as follows: Among them, Λ (ns, represents the Lagrangian factor, the penalty factor ρ (n,s) Defined as 5. The coupled tensor decomposition algorithm for multi-group high-dimensional EEG data analysis according to claim 1, characterized in that: The update formula (4) in step (4) is as follows: The update formula (5) is as follows: The update formula (6) is as follows:

6. The coupled tensor decomposition algorithm for multi-group high-dimensional EEG data analysis according to claim 1, characterized in that: The update formula (7) in step (6) is as follows:

7. The coupled tensor decomposition algorithm for multi-group high-dimensional EEG data analysis according to claim 1, characterized in that: The Lagrangian function in step (8) is expressed as follows: For formula (8), and Γ (s) Derivative, through mathematical calculations can get the corresponding solution.

8. The coupled tensor decomposition algorithm for multi-group high-dimensional EEG data analysis according to claim 1, characterized in that: The update formula (9) in step (9) is as follows: Among them, "vec" represents the vectorization operation of the tensor, and "diag" represents taking the superdiagonal elements of the tensor to form a vector.

9. The coupled tensor decomposition algorithm for multi-group high-dimensional EEG data analysis according to claim 1, characterized in that: The update formula (10) in step (10) is as follows:

Images