Sparse low-rank coupling tensor decomposition method suitable for multi-frequency dynamic function network analysis

Abstract

The invention discloses a coupling tensor decomposition method based on sparse low-rank constraint, which is used for characteristic decomposition of a multi-frequency dynamic function network connection tensor in resting state function magnetic resonance imaging data. According to the algorithm, on the basis of the traditional coupling canonical factorization (CCPD), an optimization model of sparse and low-rank constraint is constructed, and the sparse and low-rank constraint optimization model is constructed by the algorithm. On the spatial connectivity dimension, redundant function connection is reduced through an L1 sparse penalty term, and the spatial specificity of the key brain network is enhanced; and in time and frequency band dimensions, low-rank regularization constraint is adopted to improve discrimination of cross-subject time sequence characteristics. Generally speaking, the method can effectively extract connectivity characteristics with statistical significance and time states of different frequency bands from dynamic function network connection tensors of multiple frequency bands, thereby effectively identifying functional connection heterogeneity characteristics between schizophrenia patients and healthy control groups.

Description

Technical Field

[0001] The present invention relates to the field of medical signal processing and analysis, and more specifically, to the study of multi-band dynamic functional network connectivity (dFNC) in resting-state functional magnetic resonance imaging (rs-fMRI) data, and in particular to a coupled tensor decomposition method and its application in multi-band dFNC analysis. Background Art

[0002] Resting-state fMRI data can provide information about the brain's spontaneous activity when not engaged in a specific task. Multi-band dFNC analysis can reveal dynamic interaction patterns in brain networks at different frequencies. In recent years, dFNC methods have been widely used to investigate functional connectivity changes in various brain diseases, study the dynamic regulatory mechanisms of brain networks, and understand the developmental evolution of the nervous system. While traditional dFNC analysis primarily focuses on connectivity patterns within a single frequency band, multi-band dFNC methods can capture more comprehensive dynamic information, helping to characterize the regulatory properties of brain functional networks at different time scales.

[0003] Current dFNC analysis methods based on time sliding windows can generate high-dimensional, multi-subject dFNC tensors, for example, a four-dimensional tensor consisting of the number of brain network components × the number of brain network components × the number of subjects × the number of time windows. However, most dFNC methods fail to consider the high-dimensional structure of the dFNC tensor. Tensor decomposition methods, such as canonical polyadic decomposition (CPD) and Tucker decomposition, have been widely used in high-dimensional data analysis and are gradually being introduced into dynamic functional analysis research. For example, the paper "Kuang LD, Li HQ, Zhang J, Gui Y, Zhang J. Dynamic functional network connectivity analysis in schizophrenia based on aspatiotemporal CPD framework. Journal of Neural Engineering. 2024Feb26; 21(1): 016032." introduced a dFNC analysis method based on the spatiotemporal CPD framework, applying sparse low-rank CPD to a three-dimensional dFNC tensor (number of brain components × number of brain components × windows of all subjects), decomposing it into a set of shared dynamic spatial module basis matrices and time-varying weight matrices to capture the dynamic changes of brain functional networks in time and space, revealing that the dynamic functional network connectivity of patients with schizophrenia has significant differences in different states, and compared with the healthy control group, their brain functional network conversion patterns and connection strengths are also different.

[0004] However, existing dFNC analysis methods based on tensor decomposition only consider the temporal and spatial dimensions, without factoring in the effects of frequency. Consequently, they cannot fully reveal the dynamic characteristics of functional connectivity across multiple frequency bands. Coupled tensor decomposition allows different tensors to share one or more factor matrices, while allowing for differences in other dimensions. This approach is suitable for the joint analysis of dFNC tensors across different frequency bands. Therefore, developing multi-band dFNC analysis methods based on coupled tensor decomposition is of great significance. Summary of the Invention

[0005] The present invention proposes a sparse low-rank coupled canonical polyadic decomposition (CCPD) method suitable for multi-band dFNC, the purpose of which is to extract significant spatial connectivity and temporal functional connectivity difference information in different frequency bands, and provide new biomarkers for brain cognition and brain disease research. The technical solution of the present invention is to construct a three-dimensional dFNC tensor with N dimensions of (spatial connection pairs × product of time and frequency × subjects) on the basis of CCPD, introduce L1 regularization sparsity constraints in spatial connections, and impose low-rank constraints on time and frequency dimensions to extract stable and statistically significant dynamic functional network modules. Among them, the L1 regularization sparsity constraint can reduce spatial redundancy and retain important spatial connections, and the low-rank constraint can increase the difference in temporal state changes of time-varying weights at different frequencies. The method of the present invention can obtain dynamic spatial modules and temporal states with significant statistical differences, which is significantly improved compared with the original CCPD decomposition model. The specific steps of the method of the present invention are as follows:

[0006] The first step is to input N three-dimensional multi-subject multi-band dFNC tensors The number of dynamic modules R and the number of groups N. Among them, the number of combinations of brain network components is I, and the product of the time window and the number of frequency bands of the nth group is J n , the number of subjects in group n is K n ; Randomly initialize the shared space connection matrix Group-specific time and frequency weight matrices and the group-specific subject strength matrix The CCPD model of the present invention can be used to (n) Decomposed into a shared spatial connection matrix A, group-specific time and frequency weights B (n) and group-specific subject intensity C (n) The CCPD model with L1 regularization, sparse constraints and low-rank constraints is expressed as:

[0007]

[0008] in, Both λ and γ are constants, representing the penalty coefficients of the sparse constraint and the low-rank constraint, respectively. “||g|| F " and "||g||1" represent the Frobenius norm and L1 norm respectively. Let the number of iterations iters = 0, and according to formula (1), the initial error can be calculated And a preset error threshold ε1 is given;

[0009] Step 2: Update the shared space connection matrix A with L1 norm sparsity constraint. According to the idea of alternating least squares (ALS) algorithm, when matrix B (n) and C (n) When fixed, the shared spatial connectivity matrix A with L1-norm sparsity constraint is updated as follows:

[0010]

[0011] in,"(·) T " represents the transpose of the matrix, "(·) -1 " represents the inverse of the matrix, the operation "*" represents the Hardmard product, and λ is a constant. is the tensor X (n) The i-th one-dimensional slice matrix,

[0012] Step 3: Update the group time and frequency weight matrix B with low-rank approximation constraints (n) The specific steps are as follows:

[0013] (3.1) Initialize the relaxation matrix and the Lagrange multiplier matrix Hyperparameter σ (n) and hyperparameters Iteration error threshold ε min and the maximum number of iterations iter max Initialize Q (n) =0, Y (n) =0,σ (n) , Here, the relaxation matrix Q (n) Used to approximate the matrix B (n) , now put A and C (n) Treat it as a constant and update the time-varying weight matrix B at different frequencies (n) ,So:

[0014]

[0015] Among them, “||g|| * ” denotes the nuclear norm, “e” denotes the Khatri-Rao product;

[0016] (3.2) Since rank(B (n) ) is a non-convex term, this method uses B (n) The nuclear norm is used to approximate B (n) The low rank of ||B (n) || * . Update the relaxation matrix Q (n)and the Lagrange multiplier matrix Y (n) , as shown in formula (4) and formula (5):

[0017]

[0018] Y (n) =Y (n) +σ (n) (Q (n) -B (n) ) (5)

[0019] You can get B (n) The update is as follows:

[0020]

[0021] (3.3) Hyperparameter σ (n) The update is as follows:

[0022]

[0023] (3.4) Iteration number iter = iter + 1, calculate the iteration error ε iter :

[0024]

[0025] (3.5) If ε iter Greater than the preset error threshold ε min And the number of iterations iter is less than the maximum number of iterations iter max , execute (3.2).

[0026] Step 4: Update the intensity matrix C of the specific group of subjects (n) . Factor matrix C (n) The conditional least squares update for is as follows:

[0027]

[0028] The matrix

[0029]

[0030] is a tensor X (n) The k-th one-dimensional slice matrix;

[0031] Step 5: Calculate the error ε of the current iteration iters , if the iteration termination condition is not met, that is, the iteration error ε iters If the error is greater than the preset error threshold ε1 and the number of iterations iters is less than the preset maximum number of iterations, jump to the second step. If the iteration termination condition is met, jump to the sixth step.

[0032] Step 6: Output the shared spatial connection matrix A and the specific group time and frequency weight matrix B (n) and the group-specific subject intensity matrix C (n) .

[0033] The main effect and advantage achieved by the present invention lies in its ability to effectively remove redundant information between dynamic spatial modules and improve the statistical significance of temporal states. When analyzing resting-state fMRI data from 74 healthy subjects and 71 patients with schizophrenia in the publicly available COBRE dataset, the sparse low-rank constrained CCPD algorithm proposed in this invention demonstrated superior decomposition performance, while maintaining identical initial conditions and hyperparameter settings. First, among the five functional network modules decomposed by this method, three to four modules showed significant statistical differences between healthy subjects and patients with schizophrenia. In contrast, the CCPD decomposition method without the sparsity constraint only detected significant differences between the two groups in one or two spatial modules, demonstrating that the present method can more effectively extract spatial information that distinguishes group characteristics and improve the ability to discriminate between different subjects in the spatial dimension. Second, the CCPD algorithm of the present invention detected significant statistical differences between temporal states in all three frequency bands. In frequency band 1, which contains the majority of the energy of the temporal process, SLRCCPD produced more significantly different state pairs than the traditional CCPD decomposition method without the low-rank constraint. This result shows that low-rank constraints can effectively improve the classification ability of temporal states and make the differences between different temporal states more prominent. In summary, the sparse low-rank constrained CCPD algorithm proposed in this invention can extract key information that helps to identify the characteristics of different groups in spatial, temporal and frequency dimensions. Compared with the traditional CCPD decomposition method, it can more accurately capture the differences in dynamic functional network connectivity between schizophrenia patients and healthy subjects. The introduction of this method provides a new tool and idea for the subsequent multi-band dFNC analysis of large-scale fMRI data, which has important research and application value. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 This is a workflow diagram of the present invention suitable for multi-frequency dynamic functional network connection analysis. DETAILED DESCRIPTION

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

[0036] This study collected resting-state functional magnetic resonance imaging (fMRI) data from 74 healthy controls and 71 patients with schizophrenia, totaling 145 samples. Dynamic functional network connectivity analysis was performed on all subjects using a sliding time window technique and filter bank connection method, constructing three-dimensional multi-frequency dFNC tensors for healthy and schizophrenia patients. The dimensional parameters of these two tensors were C×S1×K1 and C×S2×K2, respectively (where C=820 represents the connectivity between shared brain components, S1=S2=348 represents the product of the time window and the number of frequency bands for healthy and schizophrenia patients, respectively, K1=74 represents the number of healthy subjects, and K2=71 represents the number of schizophrenia patients). The resulting tensor dimensions were 820×348×74 and 820×348×71, respectively. As shown in the accompanying figure, the novel tensor decomposition algorithm designed in this study primarily includes the following processing steps:

[0037] The first step is to input the three-dimensional dFNC tensor and The number of dynamic modules R = 5 and the number of combinations N = 2; the shared space connection matrix is randomly initialized Group-specific time and frequency weight matrices and the group-specific subject strength matrix Set the number of iterations iters = 0 and the error threshold ε1 = 10 -8 , calculate the initial iteration error ε0 according to formula (1);

[0038] In the second step, according to formula (2), the update formula of the shared space connection matrix A with L1 norm sparsity constraint is obtained. By defining the matrix Simplify the calculation and finally normalize each factor of A to obtain the shared space connection matrix

[0039] The third step is to update the group time and frequency weight matrix B with low-rank approximation constraints. (n) The specific steps are as follows:

[0040] (3.1) Initialize the relaxation matrix and the Lagrange multiplier matrix Hyperparameter σ (n) =10 -6 , and hyperparameters Iteration error threshold ε min =10 -8 , iter = 1 and the maximum number of iterations iter max =1000;

[0041] (3.2) According to formula (4) and formula (5), the relaxation matrix Q is updated respectively (n) and the Lagrange multiplier matrix Y(n) , according to formula (6), update the group time and frequency weight matrix B under low rank approximation (n) ;

[0042] (3.3) According to formula (7), update the hyperparameter σ (n) ;

[0043] (3.4) Let the number of iterations iter = iter + 1, and calculate the iteration error ε according to formula (8) iter ;

[0044] (3.5) If ε iter Greater than the preset error threshold ε min And the number of iterations iter is less than the maximum number of iterations iter max =1000, execute (3.2).

[0045] The fourth step is to update the strength matrix C of the specific group of subjects according to formula (9): (n) ;

[0046] Step 5: Calculate the error ε of the current iteration iters , if the iteration condition is not met, that is, the iteration error ε iters If the error is greater than the preset error threshold ε1 and the number of iterations iters is less than the preset maximum number of iterations 200, jump to the second step. If the iteration conditions are met, jump to the sixth step.

[0047] Step 6: Output the shared space connection matrix Group-specific time and frequency weight matrices and Group-specific subject intensity matrix and

Claims

1. A sparse low-rank constrained coupled tensor decomposition method suitable for multi-frequency dynamic functional network connectivity analysis, characterized by: In the spatial connection dimension, L1 norm constraint is used to achieve sparse representation of functional connection topology, reducing the redundancy of cross-brain functional coupling in resting-state magnetic resonance imaging data; nuclear norm constraint is imposed on time and frequency dimensions to enhance the discriminability of temporal features. The following steps are involved: The first step is to input N three-dimensional multi-subject multi-band dFNC tensors The number of dynamic modules R and the number of groups N, where the number of combinations of brain network components is I, and the product of the time window and the number of frequency bands F of the nth group is J n , the number of subjects in group n is K n , randomly initialize the shared space connection matrix Group-specific time and frequency weight matrices and the group-specific subject strength matrix Calculate initial error The second step is to update the shared space connection matrix A with the L1 norm sparsity constraint added; The third step is to update the group time and frequency weight matrix B with low-rank approximation constraints. (n) ; Step 4: Update the intensity matrix C of the specific group of subjects (n) ; Step 5: Calculate the error ε of the current iteration iters If the iteration condition is not met, jump to the second step. If the iteration condition is met, jump to the sixth step. Step 6: Output the shared space connection matrix A and the specific group time and frequency weight matrix B (n) and the group-specific subject intensity matrix C (n) .

2. The method of claim 1, wherein the shared spatial connection matrix A with L1-norm sparsity constraint is updated in the sparse low-rank constrained coupled tensor decomposition method suitable for multi-frequency dynamic functional network connectivity analysis, wherein: According to the idea of alternating least squares (ALS) algorithm, when the matrix B (n) and C (n) When fixed, the shared spatial connectivity matrix A with L1-norm sparsity constraint is updated as follows: in,"(·) T " represents the transpose of the matrix, "(·) -1 " represents the inverse of the matrix, the operation "*" represents the Hardmard product, and λ is a constant. is the tensor X (n) The i-th one-dimensional slice matrix, 3. Updating the group time and frequency weight matrix B with low-rank approximation constraints in the sparse low-rank constrained coupled tensor decomposition method suitable for multi-frequency dynamic functional network connectivity analysis according to claim 1 (n) , which is characterized by The specific steps include: The first step is to initialize the relaxation matrix and the Lagrange multiplier matrix Hyperparameter σ (n) and hyperparameters Iteration error threshold ε min and the maximum number of iterations iter max ; Initialize Q (n) =0, Y (n) =0,σ (n) , Here, the relaxation matrix Q (n) Used to approximate the matrix B (n) , now put A and C (n) Treat it as a constant and update the time-varying weight matrix B at different frequencies (n) ,So: Among them, "||g|| * " represents the nuclear norm, "e" represents the Khatri-Rao product, It's X (n) The modulo 2 expansion matrix of ; In the second step, due to rank (B (n) ) is a non-convex term, this method uses B (n) The nuclear norm is used to approximate B (n) The low rank of ||B (n) || * , update the relaxation matrix Q (n) and the Lagrange multiplier matrix Y (n) , as shown in formula (3) and formula (4): AND (n) =And (n) +σ (n) (Q (n) -B (n) ) (4) You can get B (n) The update is as follows: The third step is to use the hyperparameter σ (n) The update is as follows: Step 4: Iteration number iter = iter + 1, calculate the iteration error ε iter : Step 5: If ε iter Greater than the preset error threshold ε min And the number of iterations iter is less than the maximum number of iterations iter max , proceed to step 2.

4. Updating the strength matrix C of a specific group of subjects in the sparse low-rank constrained coupled tensor decomposition method suitable for multi-frequency dynamic functional network connectivity analysis according to claim 1 (n) , characterized in that, Factor matrix C (n) The least squares update formula for is as follows: The matrix is a tensor X (n) The k-th one-dimensional slice matrix of .

Images