A L-based 1 -L 2 Sparse decomposition method for single-subject complex fMRI data based on norm

Abstract

A sparse decomposition method for single-subject complex fMRI data based on the L1-L2 norm, belonging to the field of biomedical signal processing. The present invention adopts an L1-L2 norm that is closer to the L0 norm in terms of sparsity, namely the difference between the L1 norm and the L2 norm, to constrain the spatial activation components, realizing high-performance sparse decomposition of single-subject complex fMRI data, so as to further denoise the brain spatial components and extract useful activation maps. Specifically, a rank-1 decomposition model of complex fMRI data and methods such as DCA, ADMM, and augmented Lagrangian multiplier method are used to obtain the update formula and implementation steps of the present invention. Taking the decomposition of resting-state complex fMRI data of a healthy subject as an example, compared with the cDL method, the correlation coefficients of the DMN component and the AUD component of the present invention are increased by 8.8% and 20% respectively, and the effective voxel ratios are increased by 120% and 96.7% respectively, providing a new method for brain cognition and brain disease research.

Description

Technical Field

[0001] The present invention belongs to the field of biomedical signal processing and relates to a method based on L 1 -L 2 A sparse decomposition method for single-subject complex functional magnetic resonance imaging (fMRI) data based on the norm. Background Art

[0002] fMRI functional magnetic resonance imaging is widely used in the study of brain activation characteristics due to its non-invasive and high-precision advantages. It plays an important role in revealing the internal activities of the brain and exploring brain functions in clinical research. Complete fMRI data is essentially complex, including amplitude and phase. In recent years, research on complex fMRI data has also been fruitful at home and abroad. Although the introduction of phase information will bring high noise problems, effective phase denoising methods have been proposed, including the Chinese invention patent "Lin Qiuhua, Yu Mouchuan, Gong Xiaofeng, Cong Fengyu. A post-processing denoising method for ICA analysis of complex fMRI data, CN201410191416" and the paper "MCYu, QHLin, LDKuang, XFGong, F.Cong, VD Calhoun. ICA of full complex-valued fMRI data using phase information of spatial maps. Journal of Neuroscience Methods, vol. 249, pp. 75-91, 2015". Phase denoising method distinguishes signal and noise by phase change. Useful brain activation voxels (i.e. signal voxels) have small spatial source phase changes, distributed within the range of [-π / 4,π / 4]; noise voxels correspond to other source phase changes outside the range of [-π / 4,π / 4]. Its excellent effect is that it can eliminate high noise from high-noise complex fMRI data and extract weak activation; the integrity of brain activation is significantly better than the method that only uses amplitude fMRI data (which only accounts for half of the complex fMRI data).

[0003] In recent years, spatial sparsity has been shown to be consistent with the activation characteristics of the brain. Moreover, compared with the independent component analysis (ICA) method, sparse decomposition does not require compression of fMRI data and will not cause information loss. For the "time × space" structure of single-subject fMRI data, spatial activation maps (SMs) and their corresponding time courses (TCs) components can also be extracted from them, that is, while maintaining the integrity of fMRI data, spatial-temporal brain function information that conforms to the spatial sparsity characteristics of the brain can be obtained. Therefore, the sparse decomposition method is suitable for the analysis of fMRI data and has become one of the important methods for brain function research and brain disease diagnosis.

[0004] In the existing research on sparse decomposition of fMRI data, the commonly used sparse constraints are L 0 Norm, L 1 Norm, L p Norm. L 0 The sparse norm has the best performance, but it is difficult to solve and can only be infinitely approximated. 1 , L p The norm is formally the sum of L 0 Simulation and approximation of norm. Among them, L 1 The L norm has been used for sparse decomposition of complex fMRI, and several algorithms have been developed, such as the complex-valued dictionary learning (cDL) algorithm, see the literature "A. Iqbal, MN Meziane, AK Seghouane and KA Meraim, "Adaptive complex-valued dictionary learning: Application to fMRI data analysis," Signal Processing, vol. 166, pp. 1-14, 2020." However, compared with L 0 Compared with the norm, L 1 , L p The sparsity of the norm does not achieve the desired effect.

[0005] To solve the above problems, the present invention adopts a method that is closer to L in terms of sparsity. 0 Norm of L 1 -L 2 Norm, i.e. L 1 Norm and L 2 The difference of L norm is used to perform sparse decomposition of single-subject complex fMRI data to further denoise the brain spatial components and extract useful activation maps. 1-L 2 The advantages of the norm have been verified in the field of real image denoising, but there are no reports on its application in the complex domain, and it cannot be used for the analysis of complex fMRI data with spatial sparsity. Summary of the invention

[0006] The present invention provides a method based on L 1 -L 2 The spatial activation SM component and the time course TC component are extracted from the single-subject complex fMRI data of the "time × space voxel" two-dimensional structure. 1 Compared with the norm results, the performance of SM component is significantly improved.

[0007] The technical solution adopted by the present invention is as follows:

[0008] make represents the single-subject complex fMRI data, V is the number of voxels in the brain, and T is the number of time points. The rank-1 decomposition model of complex fMRI data is:

[0009]

[0010] Where K is the number of components; is the SM matrix, and the K rows of S contain K spatial activation component vectors is the TC matrix, and the K columns of D contain K time course component vectors

[0011] First, for the kth spatial activation component s k Introducing L 1 -L 2 The sparse constraint of the norm, k = 1, ..., K, constructs the sparse decomposition model of single-subject complex fMRI data as follows:

[0012]

[0013] In the formula, “‖·‖ F ", "‖·‖ 1 ", "‖·‖ 2 ” are L F Norm, L 1 Norm and L 2 norm; λ is the sparsity constraint parameter; E k Defined as:

[0014]

[0015] Next, the Difference of Convex Algorithm (DCA) is used to solve the L of equation (2).1 -L 2 The norm is as follows:

[0016]

[0017] Among them, "<>" represents the inner product, n represents the DCA algorithm loop iteration pointer,

[0018]

[0019] yes‖ s k ‖ 2 exist Subgradient l represents the solution of the Alternating Direction Method of Multipliers (ADMM) algorithm s k d k The loop iteration pointer.

[0020] Application initialization guarantee The ADMM algorithm is used to solve equation (4) as follows:

[0021]

[0022] In the formula, yes The dual variable of .

[0023] Furthermore, the augmented Lagrangian multiplier method is used to obtain the Lagrangian function from equation (6):

[0024]

[0025] In the formula, is the Lagrange multiplier and δ is the penalty parameter.

[0026] Combining equations (3) and (7), we can use algebraic methods to update s k , Use the soft threshold method to update the dual multiplier z k ,use s k and z k Difference update step μ k , update d using algebraic methods k ; The following update formula is obtained:

[0027]

[0028]

[0029] Among them, sign() represents the sign function, ⊙ represents the vector dot product, represents a V-dimensional all-one vector.

[0030]

[0031]

[0032] According to formula (8), the update error of ADMM algorithm is calculated:

[0033]

[0034] make Calculate the DCA algorithm update error:

[0035]

[0036] Combined with formula (1), we can get the updated spatial activation component matrix s = { s 1 , s 2 ,…, s K} H and the time process component matrix D = {d 1 ,d 2 ,…,d K}.

[0037] There are two stopping conditions for iterative updates: one is the calculation error ε iter as follows:

[0038]

[0039] ε iter The iteration can be stopped when the value is less than the preset threshold; the second is based on the maximum number of iterations N total Make a judgment, the number of iterations exceeds N total Then stop the iteration.

[0040] In addition, for E in formula (3) k Perform singular value decomposition (SVD) to obtain the maximum singular value σ and its corresponding left singular vector u and right singular vector v, and initialize s k With d k as follows:

[0041] s k =σv (15)

[0042] d k =u (16)

[0043] In summary, the specific implementation steps of the present invention are as follows (see attached Figure 1 ):

[0044] Step 1: Input data. Input single-subject complex fMRI data X in the form of "time × space voxel";

[0045] Step 2: Parameter setting. Set the number of components K, the value of the sparse constraint parameter λ in formula (2), the penalty parameter δ in formula (7), and the maximum number of iterations N of the overall algorithm. total 、ADMM maximum number of iterations N ADMM 、DCA maximum number of iterations N DCA , minimum iteration error ε min ;

[0046] Step 3: Initialize S, D and parameters. Initialize the SM matrix S and TC matrix D by generating random sparse matrices, and set the iteration parameter iter = 1 and the iteration error ε 0 =1;

[0047] Step 4: Initialize component pointers: set k=1; that is, select k=1, ..., Kth vectors to be updated in sequence.

[0048] Step 5: DCA initialization. Let n = 1. Apply equation (3) to calculate E k , use equation (15) and equation (16) to initialize the SM component TC ingredients make

[0049] Step 6: ADMM initialization. Let the iteration parameter l = 1 and the Lagrange multiplier Dual Multiplier Algorithm iteration error

[0050] Step 7: ADMM update. Apply equations (8)-(11) to update the SM component. Dual Multiplier Lagrange multipliers TC ingredients

[0051] Step 8: ADMM error calculation. Calculate the ADMM algorithm error according to formula (12):

[0052] Step 9: ADMM stop condition judgment. If l>N ADMM or Jump to step 10, otherwise execute l=l+1 and jump to step 7;

[0053] Step 10: Subgradient update.

[0054] Step 11: DCA error calculation. Let Calculate the DCA algorithm error according to Equation (13)

[0055] Step 12: DCA stop condition judgment. If n > N DCA or Jump to Step 13, otherwise execute n = n + 1 and jump to Step 6;

[0056] Step 13: Component pointer judgment. If k < K, execute k = k + 1 and jump to Step 5, otherwise jump to Step 14;

[0057] Step 14: Total error calculation. Calculate the iterative error ε according to Equation (14) iter ;

[0058] Step 15: Stop condition judgment. If the iterative error ε iter is less than the preset error threshold ε min , or iter is greater than the preset maximum number of iterations N total , then jump to Step 16, otherwise execute iter = iter + 1 and jump to Step 4.

[0059] Step 16: Combine all components to form the SM matrix S = { s 1 , s 2 , …, s K} H , the TC matrix D = {d 1 , d 2 , …, d K}.

[0060] Step 17: Output the SM matrix and the TC matrix D.

[0061] The present invention focuses on the sparse decomposition problem of single-subject multiple fMRI data and provides a method based on L 1 -L 2The sparse decomposition method of single-subject complex fMRI data based on the norm can perform high-performance decomposition of single-subject complex fMRI data in the form of "time × space voxel" to obtain rich and low-noise spatial activation components, providing a basis for brain cognition and brain disease research. Taking the resting-state complex fMRI data of a healthy subject as an example, two spatial activation components of interest are selected: the default mode network (DMN) and the auditory network (AUD), and the results obtained by the method of the present invention are compared with those obtained by the cDL method. The performance indicators compared include the Smith2009 reference template The absolute value of the Pearson correlation coefficient |ρ c |, and the effective voxel ratio V ratio , defined as follows:

[0062] |ρ c |=|corr( s ref , s k )| (17)

[0063]

[0064] Where V total express s k The total prime number of V inside express s k Falling s ref The number of voxels within V ratio The larger it is, the more useful voxels there are and the less noise there is. c The larger the |, the more similar it is to the reference template and the higher the component correctness.

[0065] After the spatial activation components are three-dimensionalized, Figure 2-3 The results of the comparison are shown. As can be seen from the figure, compared with cDL, the correlation coefficients of the DMN component and the AUD component of the present invention increased by 8.8% and 20% respectively (DMN: 0.62 vs. 0.57; AUD: 0.54 vs. 0.45), and the effective voxel proportions increased by 120% and 96.7% respectively (DMN: 0.55 vs. 0.25; AUD: 0.59 vs. 0.30). BRIEF DESCRIPTION OF THE DRAWINGS

[0066] Figure 1 The present invention is implemented as a flow chart.

[0067] Figure 2 This is a comparison diagram of the DMN spatial activation components extracted by the present invention and the cDL method.

[0068] Figure 2 A is the present invention, c |=0.62,V ratio =0.55;

[0069] Figure 2 B is the cDL method, |ρ c |=0.57,V ratio =0.25;

[0070] Figure 2 C is the DMN component reference.

[0071] Figure 3 This is a comparison diagram of the AUD spatial activation components extracted by the present invention and the cDL method.

[0072] Figure 3 A is the present invention, c |=0.54,V ratio =0.59;

[0073] Figure 3 B is the cDL method, |ρ c |=0.45,V ratio =0.30;

[0074] Figure 3 C is the AUD component reference. DETAILED DESCRIPTION

[0075] A specific embodiment of the present invention is described in detail below in conjunction with the technical solution.

[0076] There is a resting-state complex fMRI data of 1 subject, containing T = 146 whole-brain scans and V = 62336 voxels in the brain. The following are the specific steps for extracting DMN components and AUD components:

[0077] Step 1: Input Data

[0078] Step 2: Parameter setting. Set K = 80, λ = 120, δ = 1, N total =15, N ADMM =20, N DCA =25,ε min =10 -6 ;

[0079] Step 3: Initialize S, D and parameters. Generate a random sparse matrix. Initialize, set the number of iterations iter = 1, the iteration error ε 0 =1;

[0080] Step 4: Initialization of component pointer: Let k = 1.

[0081] Step 5: DCA initialization. Let n = 1. Calculate E using Equation (3) k , and initialize the SM component using Equations (15) and (16) TC component Let

[0082] Step 6: ADMM initialization. Let the iteration parameter l = 1, the Lagrange multiplier dual multiplier algorithm iteration error

[0083] Step 7: ADMM update. Update the SM component using Equations (8)-(11) dual multiplier Lagrange multiplier TC component

[0084] Step 8: ADMM error calculation. Calculate the ADMM algorithm error according to Equation (12)

[0085] Step 9: ADMM stop condition judgment. If l > N ADMM or go to Step 10, otherwise execute l = l + 1 and go to Step 7;

[0086] Step 10: Subgradient update.

[0087] Step 11: DCA error calculation. Let Calculate the DCA algorithm error according to Equation (13)

[0088] Step 12: DCA stop condition judgment. If n > N DCA or go to Step 13, otherwise execute n = n + 1 and go to Step 6;

[0089] Step 13: Component pointer judgment. If k < K, execute k = k + 1 and go to Step 5, otherwise go to Step 14;

[0090] Step 14: Total error calculation. Calculate the iteration error ε according to Equation (14) iter ;

[0091] Step 15: Stop condition judgment. If the iteration error ε iter is less than the preset error threshold ε min , or iter is greater than the preset maximum number of iterations Ntotal , then jump to step 16, otherwise execute iter=iter+1 and jump to step 4;

[0092] Step 16: Combine all components to form the SM matrix S = { s 1 , s 2 ,…, s K} H , TC matrix D = {d 1 ,d 2 ,…,d K};

[0093] Step 17: Perform phase correction and phase denoising on the SM component, apply the Smith2009 reference template of the DMN component and the AUD component, select the DMN component and the AUD component decomposed by the present invention (see Appendix Figure 2-3 ), and use equations (17) and (18) to calculate the performance index.

Claims

1. A method based on L 1 -L 2 A sparse decomposition method for single-subject complex fMRI data based on the norm. It is characterized in that make represents the complex fMRI data of a single subject, V is the number of voxels in the brain, and T is the number of time points; the rank-1 decomposition model of complex fMRI data is: Where K is the number of components; is the SM matrix, and the K rows of S contain K spatial activation component vectors is the TC matrix, and the K columns of D contain K time course component vectors First, for the kth spatial activation component s k Introducing L 1 -L 2 The sparse constraint of the norm, k = 1, ..., K, constructs the sparse decomposition model of single-subject complex fMRI data as follows: In the formula, "‖·‖ F ”、"‖·‖ 1 ”、"‖·‖ 2 ” are L F Norm, L 1 Norm and L 2 norm; λ is the sparsity constraint parameter; E k Defined as: Next, the Difference of Convex Algorithm (DCA) is used to solve the L of equation (2). 1 -L 2 The norm is as follows: wherein, "<>" represents the inner product, and n represents the DCA algorithm loop iteration pointer yes‖ s k ‖ 2 exist Subgradient l represents the solution of the Alternating Direction Method of Multipliers (ADMM) algorithm s k d k Loop iteration pointer; application initialization guarantees The ADMM algorithm is used to solve equation (4) as follows: In the formula, yes The dual variable of Furthermore, the augmented Lagrangian multiplier method is adopted, and the Lagrangian function is obtained from Equation (6) as follows where is the Lagrange multiplier and δ is the penalty parameter; Combining equations (3) and (7), we can use algebraic methods to update s k , Use the soft threshold method to update the dual multiplier z k ,use s k and s k Difference update step μ k , update d using algebraic methods k ; The following update formula is obtained: Among them, sign() represents the sign function, ⊙ represents the vector dot product, represents a V-dimensional all-1 vector; According to Equation (8), calculate the ADMM algorithm update error make Calculate the DCA algorithm update error: Combined with formula (1), we can get the updated spatial activation component matrix S = { s 1 , s 2 ,…, s K } H and the time process component matrix D = {d 1 ,d 2 ,…,d K }; There are two stopping conditions for iterative updates: one is the calculation error ε iter as follows: ε iter The iteration can be stopped when the value is less than the preset threshold; the second is based on the maximum number of iterations N total Make a judgment, the number of iterations exceeds N total Then stop the iteration; In addition, for E in formula (3) k Perform singular value decomposition (SVD) to obtain the maximum singular value σ and its corresponding left singular vector u and right singular vector v, and initialize s k With d k as follows: s k =σv (15) d k =u (16)。 2. A method based on L according to claim 1. 1 -L 2 A sparse decomposition method for single-subject complex fMRI data based on the norm. It is characterized in that The specific steps are as follows Step 1: Input data; Input the single-subject complex fMRI data X in the form of "time × spatial voxel" Step 2: Parameter setting; set the number of components K, the value of the sparse constraint parameter λ in formula (2), the penalty parameter δ in formula (7), and the maximum number of iterations N of the overall algorithm. total 、ADMM maximum number of iterations N ADMM 、DCA maximum number of iterations N DCA , minimum iteration error ε min ; Step 3: Initialize S, D and parameters; Initialize the SM matrix S and TC matrix D by generating random sparse matrices, set the iteration parameter iter = 1, and the iteration error ε 0 =1; Step 4: Component pointer initialization: Let k = 1; that is, sequentially select the k = 1,..., K vectors to be updated Step 5: DCA initialization; set n = 1; apply equation (3) to calculate E k , use equation (15) and equation (16) to initialize the SM component TC ingredients make Step 6: ADMM initialization Let the iteration parameter l = 1, the Lagrange multiplier Dual Multiplier Algorithm iteration error Step 7: ADMM update: Apply equations (8)-(11) to update the SM component Dual Multiplier Lagrange multipliers TC ingredients Step 8: ADMM error calculation Calculate the ADMM algorithm error according to formula (12): Step 9: ADMM stop condition judgment; if l>N ADMM or Jump to step 10, otherwise execute l=l+1 and jump to step 7; Step 10: Subgradient update; Step 11: DCA error calculation; let According to formula (13), the DCA algorithm error is calculated Step 12: DCA stop condition judgment; if n>N DCA or Jump to step 13, otherwise execute n=n+1 and jump to step 6; Step 13: Component pointer judgment; If k < K, execute k = k + 1 and jump to Step 5, otherwise jump to Step 14 Step 14: Calculate the total error; Calculate the iterative error ε according to formula (14) iter ; Step 15: Stop condition judgment; if the iteration error ε iter Less than the preset error threshold ε min , or iter is greater than the preset maximum number of iterations N total , then jump to step 16, otherwise execute iter=iter+1 and jump to step 4; Step 16: Combine all components to form the SM matrix S = { s 1 , s 2 ,…, s K } H , TC matrix D = {d 1 ,d 2 ,…,d K }; Step 17: Output the SM matrix and the TC matrix D

Images