CN103870710B - Tensor grouping method for multi-subject fMRI data analysis - Google Patents

Abstract

The invention discloses a tensor grouping method for multi-subject fMRI data analysis and belongs to the field of fMRI data analysis. The tensor grouping method is characterized in that according to the difference between minimized subjects, and all subject fMRI data which form a large tensor originally are divided into multiple subgroups with the subjects as units so that maximum cross-correlation can be formed between time course components and between brain space activation region components of all the subjects in each subgroup. Compared with the original large tensor, subgroup tensors constructed by the subgroup fMRI data are better matched with a CP model of tensor decomposition; as a result, the time course components and the brain space activation region components shared by the multiple subjects can be obtained through decomposition, wherein the performance of the time course components and the brain space activation region components is improved, the coefficient of correlation between the task-related time course components and the priori task stimulation process can be increased by about 0.1, the prime number of noise bodies of the task-related brain space activation region components can be reduced by about 23%, and the prime number of expected activation bodies is basically unchanged. The tensor grouping method serves as a reference to solving the problem of mismatching between other types of high dimensional data and CP models.

Description

A kind of packet Tensor Method for how tested fmri data analysis

Technical field

The present invention relates to a kind of analysis method of how tested fmri data, more particularly to one kind how tested fmri data is divided The tensor resolution method of analysis.

Background technology

Functional mri (functional magnetic resonance imaging, fmri) is referred to as observing The valid window of brain because fmri technology can collect tested complete certain particular task (task, such as vision, the sense of hearing, Motion etc.) when brain function data, and there is not damaged and high spatial resolution advantage.By using blind source separating (blind Source separation, bss), the data such as independent component analysis (independent component analysis, ica) Drive analysis method, need not any prior information, just can estimate multiple (generally several under particular task from fmri data Ten) active region (spatial activations) composition and its time course (time courses) composition between Naokong are brain Functional analysis and clinical diagnosis provide full and accurate foundation.

Fmri data is high dimensional data.Wherein, single tested fmri data is 4 dimensions, such as 165 × 53 × 63 × 46, including 3-dimensional Full brain data (53 × 63 × 46) and 1 dimension scanning times (165)；How tested fmri data is 5 dimensions, such as 165 × 53 × 63 × 46 × 16, in the single tested fmri data basis of 4 dimensions, increased the tested number of 1 dimension (16) again.In general, people can be by 3-dimensional Full brain data is launched into one-dimensional voxel (voxels) data (53 × 63 × 46=153594), how tested fmri data at this moment There is as many as 3-dimensional (16 × 153594 × 165).Therefore, it is suitable to process the tensor resolution method of higher-dimension (>=3-dimensional) data how tested Fmri data analysis aspect has very big potentiality.

At present, the cp(canonical polyadic of tensor resolution) model has been used for fmri data analysis.Cp model Assume that tensor data has parallel factor structure, for the how tested fmri data being launched into 3-dimensional, that is, assume that How tested active region composition and its time course composition between identical Naokong is shared with different intensity, specific as follows:

In formula,Represent the tensor being made up of whole how tested fmri data, a_jAnd s_jRepresent tested shared when Between active region composition between procedure component and Naokong, j is the number of composition, and cj represents strength difference, and ε represents residual error.J time Procedure component a₁,a₂,…,a_jActive region ingredient s and between j Naokong₁,s₂,…,s_jIn, there is a higher task of signal to noise ratio Related (task-related) composition, we are designated as a_*And s_*.Below with a_*And s_*As a example how tested fmri data and cp mould are described The mismatch problems of type.Assume total m tested, cp model hypothesis $a_{*}^{\left(1\right)}=a_{*}^2=...=a_{*}^{\left(m\right)},$ $s_{*}^{\left(1\right)}=s_{*}^{\left(2\right)}=...=s_{*}^{\left(m\right)}.$ However, now there are some researches show, even if under identical task, different tested time courses become Divide and also can assume larger difference between the composition of active region and between Naokong, cause $a_{*}^{\left(1\right)}\&notequal;a_{*}^2\&notequal;...\&notequal;a_{*}^{\left(m\right)},$ $s_{*}^{\left(1\right)}\&notequal;s_{*}^{\left(2\right)}\&notequal;...\&notequal;s_{*}^{\left(m\right)},$ That is, how tested fmri data and cp model have mismatch problems, therefore tensor resolution Method is difficult to play performance advantage in current how tested fmri data analysis.

Content of the invention

It is an object of the invention to, a kind of scheme solving how tested fmri data and cp model mismatch is provided, minimizes The otherness of tested, lifts performance in how tested fmri data analysis for the tensor resolution method.

The technical scheme is that, the big tensor x that m tested fmri data is built, with quilt Try to be divided into k subgroup for unit Each subgroup comprises n_kIndividual tested, 2≤n_k＜ m；Tested cross-correlation of k subgroup Coefficient mean value sumM kind packet scheme is maximum,For each in subgroup k Cross-correlation coefficient mean value between the composition of active region, k=1 ..., k between tested task correlation time procedure component and Naokong； Each subgroup tensor is carried out respectively with cp decomposition:

Obtain the shared time course composition a of tested of each subgroup_j(k)Active region ingredient s and between Naokong_j(k), tested Strength difference c_j(k)And residual error_(k).Fundamental block diagram is as shown in Figure 1.Still with the related time course composition of task and brain spatial activation As a example area's composition, due to tested time course composition each in subgroup k between and Naokong between there is between the composition of active region maximum Correlation, then $a_{*\left(k\right)}^{\left(1\right)}\≈a_{*\left(k\right)}^{\left(2\right)}\≈...\≈a_{*\left(k\right)}^{\left(n_k\right)},$ $s_{*\left(k\right)}^{\left(1\right)}\≈s_{*\left(k\right)}^{\left(2\right)}\≈...\≈s_{*\left(k\right)}^{\left(n_k\right)},$ Set up.It can be seen that, subgroup tensor dataWith the matching degree of cp model than original big tensor dataIt is improved, so, the son that subgroup tensor resolution obtains Group is shared active region composition between time course composition and Naokong and will be got a promotion in performance, is analyzed and processed by further, As average it is possible to obtain that performance is obviously improved, for active region between all tested common time course compositions and Naokong Composition.

Need in packet to use active region composition between the time course composition of each tested fmri data and Naokong, the present invention Separated using ica method and obtain.As described above, ica need not any prior information, just can estimate specific from fmri data Active region composition and its time course composition between the Naokong under business.Specifically, what cross-correlation coefficient calculated is tested to liking m The task correlation time procedure component of fmri dataAnd) active region composition between NaokongAsk forWith)Covariance matrix r_aAnd r_s:

$r_a=\left[\begin{array}{cccc}1& r_{a_{*}^{\left(1\right)}{a}_{*}^{\left(2\right)}}& \·\·\·& r_{a_{*}^{\left(1\right)}{a}_{*}^{\left(m\right)}}\\ r_{a_{*}^{\left(2\right)}{a}_{*}^{\left(1\right)}}& 1& \·\·\·& r_{a_{*}^{\left(2\right)}{a}_{*}^{\left(m\right)}}\\ \·& \·& \·& \·\\ \·& \·& 1\·& \·\\ \·& \·& \·& \·\\ r_{a_{\text{*}}^{\left(m\right)}{a}_{*}^{\left(1\right)}}& r_{a_{*}^{\left(m\right)}{a}_{*}^{\left(2\right)}}& \·\·\·& 1\end{array}\right]$

Obtain different tested p, q, p, q ∈ { 1,2 ... m } is mutual between the composition of active region between time course composition and Naokong Coefficient correlationFor reducing in subgroup active region composition between tested time course composition and Naokong simultaneously Otherness, by covariance matrix r_aAnd r_sBuild synthetical matrix r=0.5 | r_a|+0.5|r_s|, the element in r is equal toBetween andBetween positive correlation coefficient mean value.

The basic operation of packet is to find mostly concerned on active region composition between time course composition and Naokong two Tested, referred to as pairwise correlation tested { p, q } find: started with initially tested p, by pth row element zero setting in r (i.e. find except p it Outer other are tested), in pth column element, the line number corresponding to maximizing, is designated as q, then find mostly concerned with tested p Newly tested q；Repeat { p, q } and find (n_k- 1) secondary, initially tested p=is not grouped tested middle numbering minimum to find order for the first time, Find afterwards and make p=q, then build and complete a subgroup tested containing nk.Such as n_k=4, that finds for the first time is initially tested For p, then pass through 3 times { p, q } and find: { p, q } → { q, u } → { u, v }, can build containing { p, q, u, v } four tested subgroups. If w is not to be grouped tested middle numbering minimum, the then initially tested p=w of next subgroup, repeats subgroup and build operation, you can build New subgroup, such as { w, x, y, z }.Repeat subgroup to build k time, build and complete k subgroup, then a kind of packet scheme constructses complete.

For a kind of k subgroup of packet scheme, still it cannot be guaranteed that the optimality of packet.For this reason, for m tested feelings Condition, the present invention builds m kind packet scheme.Specific practice is to change initially tested point that the first subgroup interior { p, q } for the first time is found Wei 1,2 ..., m, then calculate tested cross-correlation coefficient mean value sum of k subgroup of every kind of packet scheme For tested task correlation time procedure component each in subgroup k and brain spatial activation Cross-correlation coefficient mean value between area's composition, k=1 ..., k, that is, calculate every kind of packet scheme k subgroup covariance square Battle array r₁,r₂,…,r_kMean value sum α of upper triangle element₁,α₂,…,α_m, the maximum packet scheme of value preset is final packet knot Really.

The effect that the present invention is reached and benefit are to be used for building a tensor by all how tested fmri data with existing Method compare, the present invention can improve the matching degree of each subgroup fmri data and cp model, and then is obviously improved tensor resolution The performance of method., than original big taking total (common) the task Related Component extracting how tested fmri data as a example Amount method, the task correlation time procedure component that the present invention obtains improves about 0.1 with the coefficient correlation of priori task stimulating course Left and right, and the noise number of voxels of the task related brain spatial activation area composition obtaining have dropped about 23% it is desirable to activate number of voxels then It is basically unchanged.Additionally, this method to solve other types high dimensional data (as sensor array data, biomedical data etc.) with The mismatch problems of cp tensor model have reference role.

Brief description

Fig. 1 is the fundamental block diagram of the present invention.

Fig. 2 is that the present invention analyzes how tested fmri data workflow figure.

Fig. 3 is the flow chart of task Related Component cross-correlation calculation.

Specific embodiment

With reference to technical scheme and accompanying drawing, describe a specific embodiment of the present invention in detail.

Assume 16 tested fmri data of collection under existing Motor stimulation, each tested scanning times is 165.The present invention Be divided into 4 subgroups, each subgroup comprise 4 tested, and then constitute 4 subgroup tensorsWarp The task Related Component that each subgroup is shared is obtained respectively, including time course composition a after tensor resolution_*(1),a_*(2),a_*(3),a_*(4) Active region ingredient s and between Naokong_*(1),s_*(2),s_*(3),s_*(4).Concrete steps are as shown in Figure 2.

The first step, tested for tested 1- 16 fmri data is carried out pca(principle component respectively Analysis) dimensionality reduction.According to the component number estimated result to this Motor stimulation fmri data set for the existing document, will be each using pca Tested fmri data is down to 20 dimensions from 165 dimensions, that is, assume that fmri data contains 20 independent elements.

Second step, carries out ica to each tested fmri data after pca dimensionality reduction.Performance preferably information maximization can be adopted Algorithm, then selects task Related Component, including time course composition from detached 20 compositions of icaWith Active region composition between Naokong

3rd step, the time course composition to task correlationActive region composition and between Naokong6 points) do not carry out cross-correlation calculation, flow chart is as shown in Figure 3:

(1) ask forCovariance matrix r_aAnd r_s, calculate synthetical matrix r=0.5 | r_a| +0.5|r_s|.

(2) make packet scheme pointer m=1.

(3) make child pointer k=1, initially tested p=m, repeat pairwise correlation tested { p, q } find 3 times, obtain { 1,15 } → { 15,9 } → { 9,7 }, subgroup 1:{ 1,15,9,7 } build and finish.

(4) finish if subgroup is unstructured, i.e. k < 4, make k=k+1, enter next subgroup and build: initially tested p=is not grouped for order The minimum of a value of tested numbering, even p=2, repeats the subgroup structure operation that pairwise correlation tested { p, q } finds 3 times, acquisition subgroup 2: {2,11,10,12}.Build subgroup 3:{ 3,16,6,8 in the same manner } and subgroup 4:{ 4,5,13,14.So far, k=4, is unsatisfactory for k < 4, M kind packet scheme constructses finish, and comprise 4 subgroups: { 1,15,9,7 }, { 2,11,10,12 }, { 3,16,6,8 }, and 4,5,13, 14}.Make m=m+1, enter next and be grouped scheme constructses.

(5) judge whether m < 16 sets up, if so, repeat (3), (4), can build and complete next packet scheme.For example, if Determine second and be grouped the initially tested of scheme first subgroup (k=1) 4 new subgroups: { 2,15,9,7 } can be obtained for p=2, Isosorbide-5-Nitrae, 11,10 }, { 3,16,6,8 }, { 5,12,13,14 }.By that analogy, until obtaining 16 kinds of packet schemes.

(6) calculate under 16 kinds of packet schemes, tested cross-correlation coefficient mean value sum of 4 subgroupsFor tested task correlation time procedure component each in subgroup k and brain spatial activation Cross-correlation coefficient mean value between area's composition, k=1 ..., 4, that is, calculate every kind of 4 subgroup covariance squares of packet scheme Battle array r₁,r₂,…,r_kMean value sum α of upper triangle element₁,α₂,…,α_m, the maximum packet scheme of value preset is final packet knot Really.Here, what value preset was maximum is packet scheme 7:{ 7,10,12,9 }, { 1,15,4,11 }, { 2,3,16,6 }, { 5,8,13,14 }, Final group result as the present invention.

4th step, according to above-mentioned group result, by tested for 16 after pca dimensionality reduction fmri data be divided into four subgroups 7,10, 12,9 }, { 1,15,4,11 }, { 2,3,16,6 }, { 5,8,13,14 }, and then build 4 subgroup tensors

5th step, using tensor resolution algorithm, such as comfac (complex parallel factor analysis) calculates 4 subgroup tensors are decomposed by method, respectively obtain active region composition between 20 shared time course compositions and Naokong, from Middle selection task Related Component, including time course composition a_*(1),a_*(2),a_*(3),a_*(4)Active region ingredient s and between Naokong_*(1), s_*(2),s_*(3),s_*(4).If carrying out respectively to it averagely, 16 tested common task correlation time procedure components and brain being obtained Spatial activation area composition.

Claims (1)

1. a kind of packet Tensor Method for how tested fmri data analysis, is characterized in that, by m tested fmri data The big tensor buildingIt is divided into k subgroup in units of tested Each subgroup comprises n_kIndividual tested, 2 ≤n_k<m；Tested cross-correlation coefficient mean value sum of k subgroupThis mean value sum M kind packet scheme is maximum,For tested task correlation time procedure component each in subgroup k and brain spatial activation Cross-correlation coefficient mean value between area's composition, k=1 ..., k；Each subgroup tensor is carried out respectively with cp decomposition:Obtain the shared time course composition a of tested of each subgroup_j(k)And Naokong Between active region ingredient s_j(k), the strength difference c of tested_j(k)And residual epsilon_(k)；J is the number of composition；

Cross-correlation coefficient calculate to as if m tested fmri data task correlation time procedure componentAnd active region composition between NaokongSeparated by ica method and obtain；Ask forWithCovariance matrix r_aAnd r_s, obtain different tested p, q, p, q ∈ { 1,2 ... m } Cross-correlation coefficient between the composition of active region between time course composition and NaokongWithCalculate synthetical matrix r= 0.5|r_a|+0.5|r_s|, carry out pairwise correlation tested { p, q } and find: started with initially tested p, by pth row element zero setting in r, Line number corresponding to maximizing in pth column element, is designated as q, then find the newly tested q mostly concerned with tested p；Repeat { p, q } finds n_k- 1 time, initially tested p=is not grouped tested middle numbering minimum to find order for the first time, finds afterwards and makes p=q, Then build and complete containing n_kAn individual tested subgroup；Repeat subgroup to build k time, build and complete k subgroup, then a kind of packet side Case builds and completes；

When m is tested, change in the first subgroup that { p, q } for the first time find initially tested be respectively 1,2 ..., m, Build altogether m kind packet scheme, then calculate every kind of packet scheme k subgroup covariance matrix r₁,r₂,…,r_kUpper triangle element Mean value sum α₁,α₂,…,α_m, the maximum packet scheme of value preset is final group result.