Fast coupling tensor decomposition algorithm suitable for conjoint analysis of brain imaging data

Abstract

The invention discloses a fast coupling tensor decomposition algorithm suitable for conjoint analysis of brain imaging data. According to the method, a multi-dimensional tensor joint decomposition model is constructed by adopting a CP model. By introducing low-rank approximation and replacing large-size tensor data with a small-size factor matrix, the calculation complexity of the algorithm is remarkably reduced. In the decomposition process, the coupling constraint is utilized to enable different tested brain imaging data to realize consistency and coordination on a shared factor matrix, so that effective conjoint analysis is carried out. And gradually updating the factor matrix through an APG iterative optimization algorithm, and minimizing a reconstruction error, thereby efficiently separating and extracting common and personalized characteristics in the imaging data of the multiple tested brains. The fast coupling tensor decomposition algorithm provides an efficient and reliable method for conjoint analysis of multi-tested brain imaging data. Through the algorithm, meaningful neural activity modes can be rapidly extracted from complex brain imaging data, and powerful support is provided for neuroscience research and clinical application.

Description

Technical Field

[0001] The present invention belongs to the technical field of biomedical signal processing, and in particular relates to a fast coupled tensor decomposition algorithm suitable for joint analysis of brain imaging data. Background Art

[0002] In scientific or clinical research, it is often necessary to jointly analyze brain imaging data of multiple subjects to extract universal brain activity patterns, which is of great significance to improve our understanding of brain function and the pathology of many brain diseases. Brain imaging data has the characteristics of high dimension, multiple subjects, multi-modality and multiple coupling. How to use the strong complementary information between data to make up for the defects of single data analysis, maintain the intrinsic high-dimensional multilinear structural information of the data, and accurately and quickly extract and represent the components of interest in the data, and then reveal the information flow of brain activity under specific cognitive functions, is the key to the development of brain science research technology today.

[0003] The tensor structure has strong computational characteristics and expressive power. It can well maintain the intrinsic structural information of high-dimensional spatial data and fully explore the interaction between the dimensions within the data. The corresponding tensor decomposition algorithm has been widely used in the field of brain imaging data processing. However, most tensor decomposition algorithms are currently only applicable to single-modal data analysis. For multi-subject data, it cannot effectively utilize the correlation information and complementary information between data. Although joint blind source separation methods such as group independent component analysis can complete the joint analysis of multiple brain imaging data, such two-dimensional methods usually simply splice or stack the original data into a matrix for processing. Such operations will inevitably lose the potential interaction between the dimensions of high-dimensional data and easily lead to overfitting.

[0004] Coupled tensor decomposition can be seen as an extension of tensor decomposition in processing multi-set data. This method uses the coupling information and difference information between data to provide a more accurate and stable solution for the joint analysis of high-dimensional data. In particular, in the joint analysis of multi-subject brain imaging data, coupled tensor decomposition technology has begun to attract people's attention due to the coupling relationship between data. However, research in this area is still in its infancy, with limited research results. In addition, there are many problems that need to be solved in the practical application of mathematical models and algorithms, such as: (1) The existing coupled tensor decomposition model is relatively simple and does not match the actual data well; (2) The existing coupled tensor decomposition algorithm fails to fully exploit the prior constraint information of the data for effective use; (3) The existing algorithm is only applicable to small and medium-sized tensor data. Summary of the invention

[0005] The purpose of the present invention is to provide a fast coupled tensor decomposition algorithm suitable for joint analysis of multi-subject brain imaging data. By introducing fast calculation and coupling mechanisms, it can not only effectively separate and extract common and individual features in multi-subject tensor data, but also reduce the consumption of computing resources and processing time while maintaining data integrity, thereby better supporting brain function research and neuroscience applications.

[0006] The technical solution of the present invention is to use the CANDECOMP / PARAFAC (CP) model to construct a multidimensional tensor joint decomposition model; by introducing low-rank approximation, a small-size factor matrix is ​​used to replace large-size tensor data to reduce the algorithm calculation complexity; in the decomposition process, coupling constraints are introduced to make different tensor data consistent and coordinated on the shared factor matrix, thereby realizing joint analysis; non-negativity constraints are imposed to further extract features with physiological significance; through the APG iterative optimization algorithm, the factor matrix is ​​gradually updated to minimize the reconstruction error and extract the components of interest. The specific implementation includes the following steps:

[0007] (1) Input non-negative brain imaging tensor data s=1,2,…,S, S represents the number of subjects; N represents the number of factor matrices; I n represents the number of components of the corresponding factor matrix;

[0008] (2) Each tensor satisfy in represents the sth kernel tensor, with non-zero elements located on its superdiagonal; the factor matrix A (n,s) It consists of two parts: and 0≤L n ≤R (s) ; Use the CP model to build a multi-dimensional tensor joint decomposition model, that is, assume that all tensors share corresponds to the coupling information between tensors, and Corresponding to the private information in each tensor;

[0009] (3) Based on the above model, the cost function is constructed to obtain the tensor data 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:

[0010]

[0011] (4) The alternating proximal gradient strategy is used to solve the cost function, that is, to alternately update a factor matrix and a core tensor while keeping the others unchanged; for the common components between tensors, coupling constraints are imposed as needed for joint update, and for the private components of the tensors, they are updated separately. (n,s) and the core tensor The update is explained: For the core tensor Use the following cost function to solve

[0012]

[0013] Core Tensor The update can be solved by the following closed-form solution:

[0014]

[0015] in, is the extrapolation point, is the Lipschitz constant, The cost function formula (1) is about The derivative of :

[0016]

[0017] Where "vec" represents the vectorization operation of the tensor, and "ddiag" represents the vector formed by the superdiagonal elements of the tensor;

[0018] For the factor matrix A (n,s) , assuming that it is and It consists of two parts and is solved using the following cost function:

[0019]

[0020] After calculation, and The closed-form solutions are:

[0021]

[0022] in, is the extrapolation point, is the Lipschitz constant, Formula (6) is about The derivative of , can be calculated as:

[0023]

[0024] (5) A low-rank approximation strategy is used for dimensionality reduction, and an unconstrained tensor decomposition algorithm is used for optimization. Available Then use the estimated factor matrix To replace the tensor in the cost function formula (1)

[0025] (6) Therefore, the cost function in step (3) can be converted to:

[0026]

[0027] (7) In formula (4) and formula (7), and The computational complexity is very large. Through low-rank approximation, the two formulas can bypass direct operations on tensors and indirectly use small factor matrices for calculations:

[0028]

[0029] Through the above replacement operation, The computational complexity of Ο(SRΠ n I n ) down to The computational complexity of n I n ) down to

[0030] (8) Parameter setting: Assuming it is the kth iteration, the Lipschitz constant for the k-1th iteration is defined as The weight of the extrapolated point is defined as

[0031]

[0032] Among them, δ ω =0.9999, t0=1, and The extrapolation points are defined as:

[0033]

[0034] (9) The derivation ends here. Next, define the iteration termination condition. Let the number of iterations be k, and calculate the change value l of the relative error between two adjacent iterations. 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 k is reached, the algorithm stops iterating, otherwise it returns to step (4); in this algorithm, 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;

[0035] (10) Output all factor matrices A (n,s) and the core tensor And output the corresponding brain activation map, power spectrum, time domain waveform and other information.

[0036] 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) and the APG coupled tensor decomposition algorithm that does not use low-rank approximation, in the computer simulation problem of coupled tensor decomposition of multiple tensor data, the present invention can maintain the accuracy of tensor and factor matrix decomposition while increasing the operation time by more than 5 times. Among them, for the joint analysis of event-related potential (ERP) data collected from 44 subjects, the tensor fit, operation time, and ERP multi-domain feature correlation coefficient presented after 100 runs of each algorithm are used as evaluation indicators. Compared with the above-mentioned coupled tensor decomposition algorithm that does not adopt an acceleration strategy, the present invention, under the premise that the extracted task-related ERP multi-domain features (including latency, amplitude, spatial distribution and spectrum, etc.) conform to the prior information, has a correlation coefficient that is only about 0.6% lower than the highest performance APG algorithm; at the same time, the present invention can maintain a high precision of 0.8490 in tensor fitting, which is basically consistent with the highest precision of 0.8491 and the APG algorithm that does not adopt acceleration, but in terms of running time cost, the present invention can save 350% to 450% of running time, greatly improving the running efficiency. Therefore, the present invention can effectively perform joint analysis on brain imaging data such as ERP data, and can identify and extract multi-domain features related to the task. The extraction of this information can provide a more effective means for brain science research and neuroscience applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] Figure 1 It is a workflow diagram of the present invention for jointly analyzing ERP data of multiple subjects. DETAILED DESCRIPTION

[0038] A specific example of the present invention is described in detail below in conjunction with the technical scheme and the accompanying drawings.

[0039] There are ERP data collected from 44 child subjects (21 dyslexic children and 19 attention-deficient children) from two groups under the oddball paradigm, that is, S = 44. After a series of preprocessing operations, the data of each subject is organized into a tensor form, and the data dimension is 9 (number of channels) * 71 (frequency points) * 60 (time points), that is, N = 3. According to the prior information, it is assumed that the number of components of the tensor data and the number of coupling components between the data are 36 respectively. The steps of the joint analysis of the non-negative ERP tensor data of 44 subjects using the present invention are shown in Figure 1.

[0040] Step 1: Input 44 non-negative ERP tensor data s=1,2,…,44.

[0041] Step 2: Low-rank approximation reduces data dimension: Specifically, an unconstrained tensor decomposition algorithm is used for optimization. For each tensor data, we get So we can use a small factor matrix to replace large size tensors.

[0042] Step 3: Initialization: Randomly initialize the factor matrix and core tensor, and each element obeys a uniform distribution. The factor matrix corresponding to the channel, frequency, and time is initialized to A (1,s) , A (2,s) , A (3,s) , the core tensor is initialized as

[0043] Step 4: Calculate the estimated tensor using back projection At the same time, the RelErr formula is used to calculate the initial tensor relative error RelErr0.

[0044] Step 5: Coupling decomposition based on non-negative constraints: Set and calculate the Lipschitz constant, extrapolation points and their weights, and use formulas (6) and (3) to update each factor matrix and core tensor alternately. It should be noted that for the common components between tensors, coupling constraints are imposed as needed for joint update, and for the private components of the tensors, they are updated separately.

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

[0046] Step 7: Determine whether the iteration stop condition is met: set the threshold to ε = 10 -8 , the maximum number of iterations is k max =1000. That is, when l k Satisfy k |<ε or the number of iterations reaches k max , the optimization algorithm iteration stops, otherwise it returns to step 5.

[0047] Step 8: Determine the ERP multi-domain features based on prior information, calculate the multi-domain feature correlation coefficient based on 100 runs, and output the corresponding brain activation map, power spectrum, time domain waveform and other information.

Claims

1. A fast coupled tensor decomposition algorithm for joint analysis of brain imaging data, It is characterized in that The following steps are involved: (1) Input non-negative brain imaging tensor data S represents the number of subjects; N represents the number of factor matrices; I n represents the number of components of the corresponding factor matrix; (2) Each tensor satisfy in represents the sth kernel tensor, with non-zero elements located on its superdiagonal; the factor matrix A (n,s) It consists of two parts: and The CP model is used to construct a multi-dimensional tensor joint decomposition model, assuming that all tensors share corresponds to the coupling information between tensors, and Corresponding to the private information in each tensor; (3) Based on the model described in step (2), the cost function formula (1) 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) The alternating proximal gradient strategy is used to solve the cost function formula (1), that is, the factor matrix A is updated alternately. (n,s) and the core tensor While keeping everything else unchanged; (5) A low-rank approximation strategy is used for dimensionality reduction, and an unconstrained tensor decomposition algorithm is used for optimization. have to Then use the estimated factor matrix To replace the tensor in the cost function (6) The cost function in step (3) is converted into formula (8); (7) In formula (4) and formula (7), and The calculation is complex through low-rank approximation. The two formulas bypass the direct operation of tensors and indirectly use the small factor matrix formula (9) for calculation; (8) Parameter setting: Assuming it is the kth iteration, the Lipschitz constant for the k-1th iteration is defined as The weight of the extrapolated point is defined as formula (10); (9) 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 k is reached, the algorithm stops iterating, otherwise it returns to step (4); the maximum number of iterations in the algorithm 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; (10) Output all factor matrices A (n,s) and the core tensor And output the corresponding brain activation map, power spectrum and time domain waveform information.

2. The fast coupled tensor decomposition algorithm for joint analysis of brain imaging data according to claim 1, It is characterized in that The cost function is expressed as follows:

3. The fast coupled tensor decomposition algorithm for joint analysis of brain imaging data according to claim 1, It is characterized in that Step (4) for the core tensor The following cost function is used for solving: Core Tensor The update has the following closed-form solution: in, is the extrapolation point, is the Lipschitz constant, The cost function formula (1) is about The derivative of : Where "vec" represents the vectorization operation of the tensor, and "ddiag" represents the superdiagonal elements of the tensor forming a vector.

4. The fast coupled tensor decomposition algorithm for joint analysis of brain imaging data according to claim 1, It is characterized in that Step (4) For the factor matrix A (n,s) , assuming that it is and It consists of two parts and is solved using the following cost function: After calculation, and The closed-form solutions are: in, is the extrapolation point, is the Lipschitz constant, Formula (6) is about The derivative of , can be calculated as:

5. The fast coupled tensor decomposition algorithm for joint analysis of brain imaging data according to claim 1, It is characterized in that The formula (8) is expressed as:

6. The fast coupled tensor decomposition algorithm for joint analysis of brain imaging data according to claim 1, It is characterized in that The formula (9) is expressed as:

7. The fast coupled tensor decomposition algorithm for joint analysis of brain imaging data according to claim 1, It is characterized in that The formula (10) is expressed as: Among them, δ ω =0.9999, t 0 =1, and The extrapolation points are defined as:

Images