Spatially weighted pooling invariant rank-(l,l,1,1) block term decomposition algorithm for multi-subject fMRI

Abstract

The application discloses a spatial weighted pooling shift invariant rank-(L, L, 1, 1) block term decomposition algorithm suitable for multi-subject fMRI data, and belongs to the field of medical signal processing. On the basis of an alternating least squares (ALS) method, a spatial weighted pooling processing method is proposed to preprocess multi-subject fMRI data, down-sample and smooth fMRI data images, and significantly reduce the fMRI data volume and remove most of the noise. In addition, considering the high space-time difference between subjects, the method combines spatial orthogonalization constraints and time shift invariance, relaxes and compresses the rank-(L, L, 1, 1) BTD model of the fMRI data, and improves the separation performance of the algorithm.

Description

Technical Field

[0001] This invention relates to the field of medical signal processing, and specifically to a method for separating high-dimensional multi-subject functional magnetic resonance imaging (fMRI) data. Background Technology

[0002] fMRI, with its advantages of high spatial resolution and non-invasiveness, has become a crucial tool in brain science research. Multi-subject fMRI data possesses high-dimensionality, including three-dimensional brain space, one-dimensional time, and one-dimensional subject. However, most fMRI data separation methods, such as independence component analysis (ICA) and canonical polyadic decomposition (CPD), reduce the three-dimensional brain space to one dimension, fundamentally destroying the three-dimensional brain spatial structure of the fMRI data. Rank-(L,L,1,1) block term decomposition (BTD) unfolds the three-dimensional brain spatial image into a matrix, yielding four-dimensional multi-subject fMRI data (two-dimensional brain space × time × subject dimension). Separating this data yields two factor matrices containing shared brain spatial information, a shared temporal component, and subject intensity. Therefore, rank-(L,L,1,1) BTD effectively preserves the higher-dimensional brain spatial structure characteristics and has been proven to be more robust than CPD.

[0003] However, during fMRI data acquisition, equipment gradient noise and physiological noise caused by subjects' breathing and heartbeat are unavoidable, resulting in high noise levels in fMRI data. Furthermore, multi-subject fMRI data exhibits spatiotemporal differences among subjects, specifically spatial variations in the location, size, and shape of brain activation, as well as temporal differences due to differences in reaction time. Compared to CPD, the brain space voxels in the four-dimensional multi-subject fMRI data of the rank-(L,L,1,1)BTD model include both intrabrain and extrabrain voxels, and their number is typically much greater than the product of the time dimension and the number of subject dimensions. Additionally, the rank-(L,L,1,1)BTD model includes not only the component number N but also the rank L of the brain space matrix, making it more computationally complex than CPD. Moreover, the algorithm complexity increases with the increase of N and L.

[0004] In summary, the application of rank-(L,L,1,1)BTD to multi-subject fMRI data needs to address the issues of high noise and spatial temporal variability in fMRI data, large computational load, and slow operation. Summary of the Invention

[0005] This invention proposes a spatially weighted pooling method to preprocess multi-subject fMRI data based on alternating least squares (ALS). This method downsamples and smooths the fMRI images, significantly reducing the amount of fMRI data and removing most of the noise. In addition, considering the high spatiotemporal variability among subjects, this method improves the separation performance of the algorithm by combining spatial orthogonalization constraints and time shift invariance to relax the rank-(L,L,1,1)BTD model of the compressed fMRI data.

[0006] This invention first assigns higher weights to intrabrain voxels and lower weights to extrabrain voxels, then downsamples the fMRI data, effectively reducing the amount of original four-dimensional multi-subject fMRI data without affecting the effective estimation of the factor matrix. Furthermore, spatial orthogonality constraints and time shift invariance are added to improve the overall separation performance of the algorithm. The steps are as follows:

[0007] Step 1: Input five-dimensional multi-subject fMRI data Where I1, I2, and I3 represent the number of voxels in the x, y, and z directions of the three-dimensional brain space, respectively; J represents the number of scans; K represents the number of subjects; the number of components is N (N is a positive integer and 0 < N ≤ J); and the rank L of the two-dimensional brain space matrix is ​​L (L is a positive integer and 0 < L ≤ I1).

[0008] Step 2: Spatial weighted pooling processing. Let the spatial index be... satisfy:

[0009]

[0010] α represents the voxel weights in the brain, α ≥ 0.5. Then, the three-dimensional fMRI data for each subject k (k = 1, ..., K) at each time point j (j = 1, ..., J) are... Become

[0011]

[0012] Here, "*" represents dot product. Then, the spatially weighted multi-subject fMRI data is four-dimensional to obtain the fMRI data tensor. Where I 23 =I2×I3. Finally, for Spatial brain images of each subject k (k = 1, ..., K) at each time point j (j = 1, ..., J) Perform pooling with a step size of 3 and a window size of 3×3:

[0013]

[0014] in `average_pooling(·)` performs average pooling. This yields spatially weighted pooled four-dimensional fMRI data.

[0015] Step 3: Initialize the factor matrix. Randomly initialize the factor matrix, which contains shared brain spatial information. and Shared time components and subject intensity Set the iteration number iter = 0, and calculate the iteration error ε of the shift-invariant rank-(L,L,1,1)BTD model. iter :

[0016]

[0017] in The outer product is given by "||·||", which represents the 2-norm of the tensor, i.e., the square root of the sum of the squares of all elements of the tensor. The superscript "T" indicates the matrix transpose. For the time delay of N components of K subjects, Let c represent the time component of the k-th subject, where c n (j-τ k,n ) indicates that for the k-th (k=1,…,K) subject, the n-th (n=1,…,N) time component has a time delay τ k,n .

[0018] Step 4: Update the factor matrix containing shared brain spatial information and Define the joint mixture matrix The element is z j+k(J-1),n =d k,n c n (j-τ k,n If the time-invariant rank-(L,L,1,1)BTD model is then expressed as:

[0019]

[0020]

[0021] Where “⊙” represents the matrix-level Khatri-Rao product, satisfying For Kronecker product, and They are and The row vector and They are and The row vector and They are The one-dimensional and two-dimensional expansion forms. According to equations (5) and (6), ALS is used to minimize the mean square error:

[0022]

[0023]

[0024] Therefore, according to equations (9) and (10), the factor matrix is ​​respectively... and Update them separately:

[0025]

[0026]

[0027] Among them, superscript It is a false rebellion.

[0028] Step 5: Update the pooled shared brain space activation components And standard orthogonalization constraints were added. Shared brain spatial activation components were calculated.

[0029]

[0030] Where vec(·) is the vectorization operation, and then through economical singular value decomposition, the vectorization is further processed. Standard orthogonalization process:

[0031]

[0032] in and These are the left singular matrix and the right singular matrix, respectively, and "esvd(·)" represents the reduced singular value decomposition.

[0033] Step 6: Update the shared time component C. Three-dimensional unfolding form Perform a Fourier transform to convert to the frequency domain. Then the frequency domain form of the shift-invariant rank-(L,L,1,1)BTD model can be expressed as:

[0034]

[0035] in, and yes and The f-th row vector, where f = 1, ..., F, F = J. According to equation (13), ALS is used to minimize the mean square error:

[0036]

[0037] Therefore, the update formula for the shared time component C can be obtained as follows:

[0038]

[0039] Among them, "⊙" c " is a vector-level Khatri-Rao product, satisfying exp{·} is an exponential function.

[0040] Step 7: Update the subject's latency According to the literature " al. Shift-invariant multilinear decomposition of neuroimaging data. NeuroImage 42(4):1439-1450, 2008”, based on the following least squares error for the time delay τ of each subject k and each component n (k=1,…,K,n=1,…,N). k,n Update:

[0041]

[0042] Step 8: Update the subject strength D. By performing a four-dimensional unfolding, we obtain Then the time-invariant rank-(L,L,1,1)BTD model in the frequency domain can be expressed as:

[0043]

[0044] Where, d k , and e (4)k It is D. and The k-th row, where k = 1, ..., K. According to equation (17), ALS is used to minimize the mean square error:

[0045]

[0046] The formula for updating the subject intensity D is as follows:

[0047]

[0048] Step 9: Calculate the iteration error. Let iter = iter + 1, and calculate the iteration error ε according to equation (4). iter .

[0049] Step 10: If ε iter Less than the preset error threshold ε iter_min If the above steps are not executed, proceed to step 12; otherwise, proceed to step 11.

[0050] Step 11: If iter is greater than the preset maximum number of iterations iter max If the above steps are not executed, proceed to step 12; otherwise, proceed to step 4.

[0051] Step 12: Calculate the shared spatial activation component S in the whole brain. Using the previously estimated C, And D, update the joint mixture matrix Let the rank of the two-dimensional whole-brain space be L1, and initialize it randomly. Based on (9) and (10), the following ALS update formulas for factor matrices A and B are obtained:

[0052]

[0053]

[0054] in and They are and The i1th line, and They are B and The i-th 23 Okay, X (1) and X (2) They are The one-dimensional and two-dimensional unfolding forms. After estimating A and B, the brain spatial activation components shared by the entire brain can be obtained.

[0055]

[0056] Step 13: Output factor matrices A and B containing shared brain spatial information, the whole-brain shared brain spatial activation component S, the shared temporal component C, and the subject's latency. With the intensity D of the test subject.

[0057] This invention enables rapid, stable, and efficient estimation of task-related components in multi-subject fMRI data. In the analysis of fMRI data from a finger-tapping task involving 16 participants, compared to the rank-(L,L,1,1)BTD algorithm with orthogonalization constraints, the estimated separation results are significantly improved. Furthermore, the algorithm of this invention requires approximately 8.68 times faster execution per iteration and has a complexity reduction of approximately 21 times. Compared to the ALS algorithm with orthogonal constraints, the average correlation coefficient between task-related shared brain space components and reference brain space components is increased by approximately 36.67%, and the average correlation coefficient between task-related shared temporal components and reference temporal components is increased by approximately 202.81%. This invention can rapidly and effectively extract shared brain functional information from multiple participants, showing promising prospects for future research on large-scale high-dimensional fMRI data and smart healthcare. Attached Figure Description

[0058] Appendix Figure 1 This is a flowchart of the workflow for analyzing multi-subject fMRI data in this invention. Detailed Implementation

[0059] A specific embodiment of the present invention will now be described in detail with reference to the technical solution and accompanying drawings.

[0060] We have fMRI data collected from 16 participants performing a finger-tapping task, i.e., K=16. Each participant underwent J=165 scans, each scan acquiring 53×63×46 whole-brain data. After removing some extrabrain data voxels, we retained voxels of 48×60×45, i.e., I1=48, I2=60, I3=45, to obtain five-dimensional fMRI data.

[0061] Step 1: Input five-dimensional multi-subject fMRI data The number of components is N=35, the rank of the two-dimensional brain space matrix is ​​L=15, and the weight is α=0.6.

[0062] Step 2: Spatial weighted pooling processing. According to equations (1-2), spatially weighted multi-subject fMRI data can be obtained. This spatially weighted multi-subject fMRI data is then four-dimensionalized to obtain... Then, according to equation (3), for Spatial brain images of each subject k (k = 1, ..., 165) at each time point j (j = 1, ..., 16). Perform pooling with a step size of 3 and a window size of 3×3 to obtain... Thus, we obtain the four-dimensional fMRI data after spatial weighted pooling.

[0063] Step 3: Randomly initialize four factor matrices, which contain shared brain space components. Shared time components and subject intensity Set the iteration count iter = 0 and the error threshold ε. min =10 -6 and maximum number of iterations iter max =200, calculate the iteration error ε according to equation (4). iter .

[0064] Step 4: Update the factor matrix containing shared brain spatial information. and According to equations (9) and (10), respectively, the factor matrix is... and Update.

[0065] Step 5: Update the pooled shared brain space activation components According to equation (11), update And add according to equation (12) Standard orthogonalization process.

[0066] Step 6: Update the shared time component C. Perform 3D unfolding Perform a Fourier transform on it to convert it to the frequency domain. According to equation (15), the shared time component C is updated.

[0067] Step 7: Update the subject's latency Based on equation (16), the time delay τ for each subject k and each component n is... k,n (k = 1, ..., 16, n = 1, ..., 35) are updated.

[0068] Step 8: Update the subject strength D. By performing a four-dimensional unfolding, we obtain The factor matrix D is updated according to equation (19).

[0069] Step 9: Calculate the iteration error. Let iter = iter + 1, and calculate the iteration error ε according to equation (4). iter .

[0070] Step 10: If ε iter Less than the preset error threshold ε iter_min =10 -6 If the above steps are not executed, proceed to step 12; otherwise, proceed to step 11.

[0071] Step 11: If iter is greater than the preset maximum number of iterations iter max =200, jump to step 12, otherwise execute step 4.

[0072] Step 12: Calculate the shared spatial activation component S in the whole brain. Using the previously estimated C, And D, update the joint mixture matrix Let the rank of the two-dimensional whole-brain space be L1 = 35, and initialize it randomly. According to equations (21) and (22), update the factor matrix of the whole-brain shared brain space respectively. And B, calculate the shared spatial activation components of the whole brain according to equation (23).

[0073] Step 13: Output factor matrices A and B containing shared brain spatial information, the whole-brain shared brain spatial activation component S, and the time delay. Shared time component C and subject intensity D.

Claims

1. A spatially weighted pool shift-invariant fast rank-(L,L,1,1) block term decomposition algorithm suitable for multi-subject fMRI data. Its features include the following steps: Step 1: Input five-dimensional multi-subject fMRI data ,in , , Let N represent the number of voxels in the x, y, and z directions of the three-dimensional brain space, J represent the number of scans, and K represent the number of subjects; N is the number of components, and L is the rank of the two-dimensional brain space matrix; N is a positive integer and... L is a positive integer and ; Step 2: Spatial weighted pooling processing, assuming the spatial index is... ( , , ),satisfy: (1) ； Indicates the voxel weights in the brain. Then for each subject k at each time point j ( 3D fMRI data Become (2) ； in, ;" "" represents dot product. Then, the spatially weighted multi-subject fMRI data is four-dimensional to obtain the fMRI data tensor. ,in Finally, regarding Spatial brain images of each subject k at each time point j Perform pooling with a step size of 3 and a window size of 3×3. : (3) ； in, , , , For average pooling operation, spatially weighted pooling is used to obtain four-dimensional fMRI data. , , ; Step 3: Initialize the factor matrix. Randomly initialize the factor matrix, which contains shared brain spatial information. , and , Shared time components and subject intensity Set the iteration number iter=0 and calculate the iteration error of the shift-invariant rank (L,L,1,1) BTD model. : (4) ； in," "For the outer product," "T" represents the 2-norm of the tensor, which is the square root of the sum of the squares of all elements of the tensor. The superscript "T" indicates the matrix transpose. For the time delay of N components of K subjects, Let k represent the time component of the k-th subject, where This indicates that for the nth time component of the kth subject, there is a time delay. ; ; ; Step 4: Update the factor matrix containing shared brain spatial information and Define the joint mixture matrix , where the element is Then the shift-invariant rank-(L,L,1,1) BTD model is expressed as: (5) ； (6) ； in," " represents the matrix-level Khatri-Rao product, satisfying , " "For Kronecker product, , and They are , and The row vectors ( ), , and They are , and The row vectors ( ), and They are The one-dimensional and two-dimensional expansion forms, according to equations (5) and (6), are used with ALS to minimize the mean square error: (7) ； (8) ； According to equations (9) and (10), the factor matrix is ​​respectively... and Update them separately: (9) ； (10) ； Among them, the superscript " "This is a false rebellion;" Step 5: Update the pooled shared brain space activation components Furthermore, standard orthogonalization constraints were incorporated to calculate shared brain spatial activation components. : (11) ； in, For vectorization operations, further reduction to singular value decomposition is performed. Standard orthogonalization process: (12) ； in, and These are the left singular matrix and the right singular matrix, respectively. "This is for reducing singular value decomposition; Step 6: Update the shared time component ,Will Three-dimensional unfolding form Perform a Fourier transform to convert to the frequency domain. Then, the frequency domain form of the time-invariant rank-(L,L,1,1) BTD model is expressed as: ， (13) ； in, , and yes , and The row vectors, where F=J, according to equation (13), ALS is used to minimize the mean square error: (14) ； Therefore, we can obtain the shared time component. The updated formula is: (15) ； in," " is a vector-level Khatri-Rao product, satisfying , , It is an exponential function. ; Step 7: Update subject latency Based on the following least square error for the time delay of each subject k and each component n Perform the update for k=1,…,K, n=1,…,N: (16) ； Step 8: Update the subject intensity ,Will By performing a four-dimensional unfolding, we obtain Then, the frequency domain form of the time-invariant rank-(L,L,1,1) BTD model is expressed as: ， (17) ； in, , and yes , and The Okay, among them According to equation (17), ALS is used to minimize the mean square error: (18) ； The intensity of the test subject can be obtained The updated formula is as follows: (19) ； Step 9: Calculate the iteration error. Let iter = iter + 1, and calculate the iteration error according to equation (4). ; Step 10: If Less than the preset error threshold If the above steps are incorrect, proceed to step 12; otherwise, proceed to step 11. Step 11: If Greater than the preset maximum number of iterations If the above steps are incorrect, proceed to step 12; otherwise, proceed to step 4. Step 12: Calculate the shared spatial activation components of the whole brain Using the estimated , ,and Update the joint mixture matrix Let the rank of the two-dimensional whole-brain space be L1, and initialize it randomly. Based on (9) and (10), the following ALS pair factor matrices are obtained. and Update the formula: (21) ； (22) ； in and They are and The OK, and They are and The OK, and They are The one-dimensional and two-dimensional expansion forms are estimated. and Subsequently, the brain space activation components shared by the entire brain were obtained. : (23) ； Step 13: Output a factor matrix containing shared brain spatial information , Whole brain sharing brain space activation components Shared time components Subject delay , with the intensity of the test subject .

Images