On the optimality of extended maximal length linear feedback shift register sequences

Statistics and Probability Letters 83 (2013) 1479–1483

Contents lists available at SciVerse ScienceDirect

Statistics and Probability Letters journal homepage: www.elsevier.com/locate/stapro

On the optimality of extended maximal length linear feedback shift register sequences Ming-Hung Kao ∗ School of Mathematical and Statistical Sciences, Arizona State University, Tempe, AZ 85287, United States

highlights • • • •

We provide sufficient conditions for functional MRI designs to be optimal. We present and prove the optimality of an extended m-sequence. An extended m-sequence is universally optimal for studies with one stimulus type. An extended m-sequence is D-optimal for cases with two or more stimulus types.

article

info

Article history: Received 27 January 2013 Received in revised form 15 February 2013 Accepted 15 February 2013 Available online 1 March 2013

abstract Maximal length linear feedback shift register sequences (or m-sequences) are widely used as experimental designs for functional magnetic resonance imaging (fMRI) experiments. However, the statistical optimality of these designs has not been proven. Here, we present and study an extended version of m-sequences and show that these designs are optimal. © 2013 Elsevier B.V. All rights reserved.

Keywords: de Bruijn sequences Optimal experimental designs Functional MRI m-sequences

1. Introduction Maximal length linear feedback shift register (LFSR) sequences, also known as m-sequences, are a certain type of sequences over a finite field. These sequences are widely used in, e.g., information and communication systems. To achieve efficient statistical inference, Buracas and Boynton (2002) advocated the use of m-sequences as experimental designs for functional magnetic resonance imaging (fMRI), a pioneering, noninvasive technology for studying functions of human brains. In a typical fMRI experiment, mental stimuli such as pictures or sounds of one or more types are sequentially presented to an experimental subject while an MR scanner scans the subject’s brain to collect fMRI data for statistical analysis (Lazar, 2008). As described in Buracas and Boynton (2002), an m-sequence can be used to determine the presentation order and timings of the stimuli to allow experimenters to collect informative data. Following this pioneering work, m-sequences and related designs have gained much popularity in practice. However, to our knowledge, there exist no analytical results that rigorously prove the statistical optimality of these designs through optimal design theory. In this paper, we provide some insightful analytical results on optimal fMRI designs. We derive sufficient conditions for fMRI designs to be optimal and show that an extended version of m-sequences, which will be defined in the next section, satisfies these sufficient conditions. Specifically, extended m-sequences are universally optimal in the sense of Kiefer (1975)

∗

Tel.: +1 7062015071. E-mail address: [email protected].

0167-7152/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.spl.2013.02.012

1480

M.-H. Kao / Statistics and Probability Letters 83 (2013) 1479–1483

for experiments with one stimulus type. This ensures that these designs are optimal under a very general class of optimality criteria, including the A- and D-optimality criteria which are widely considered in fMRI. A-optimal designs minimize the average variance of parameter estimates. D-optimal designs minimize the volume of a simultaneous confidence region of the parameters of interest at any specified confidence level, or equivalently, they maximize the determinant of the information matrix. For studies with Q (> 1) stimulus types, results in, e.g., Kao et al. (2009) indicate that m-sequences can be suboptimal under A-optimality. On the other hand, an empirical study by Maus et al. (2010) suggests that these designs might be D-optimal. Here, we are able to prove that extended m-sequences are D-optimal for experiments with two or more stimulus types. Such analytical results are currently unavailable. In the next section, we briefly introduce an extended version of m-sequences. Background information relevant to our study is also provided. Our main results are presented in Section 3. The paper is concluded in Section 4 with a discussion. 2. Background A sequence {dn } over a finite field GF (p) is called an LFSR sequence if it satisfies a homogeneous linear recurrence relation dn+k =

k 

ai dn+k−i ,

(1)

i=1

where k is a given integer, and a1 , . . . , ak ∈ GF (p) are given coefficients. Such a sequence is periodic and repeats itself after ℓ elements; ℓ is the least period and is the smallest integer making dn+ℓ = dn for every n ≥ 1. It can be seen that ℓ ≤ pk − 1 since there are (pk − 1) non-zero k-tuples. An LFSR sequence having the maximum least period is called an m-sequence, which is typically trimmed to have a length of (pk − 1). One way to construct an m-sequence is by a primitive polynomial of GF (p). For example, f (x) = x2 − 2x − 1 is a primitive polynomial over GF (3) and can be used to define the linear recurrence relation: dn+2 = 2dn+1 + dn . In particular, the highest degree of f (·) and the coefficient of xi respectively determine the values of k and ai in (1). An m-sequence can then be generated with a specified initial k-tuple, (d1 , . . . , dk ). For example, with dn+2 = 2dn+1 + dn and (d1 , d2 ) = (1, 1), we have the m-sequence 11012202. We note that primitive polynomials can be found in, e.g., Tables 3.6 and 3.7 of Golomb and Gong (2005) and refer readers to Lidl and Niederreiter (1994), and MacWilliams and Sloane (1977) for further information about m-sequences. A well-known property of m-sequences is that every non-zero k-tuple appears exactly once as we move a ‘window’ of length k along the sequence, treating d1 , . . . , dk−1 as the next k − 1 elements of dpk −1 . Here, we follow Kiefer and Wynn (1984) to insert an additional 0 to a (k − 1)-tuple of zeros to include a zero k-tuple. In the previous example, such a design can be 110012202 or 110122002. We will call the resulting sequences as ‘‘extended m-sequences’’. They belong to a class of de Bruijn sequences (Golomb and Gong, 2005; Aguirre et al., 2011) and are optimal experimental designs for fMRI as will be shown in Section 3. In an fMRI experiment, a mental stimulus to be presented to an experimental subject can possibly occur every τISI seconds after the outset of the experiment; here, τISI is a pre-specified time. The elements of an m-sequence can be used to determine the types of the stimuli to be presented at these time points. Specifically, the nth element dn corresponds to time (n − 1)τISI . When dn = q is a positive integer, a qth-type stimulus appears briefly at the corresponding time point, q = 1, . . . , Q (= the total number of stimulus types). If dn = 0, there is no stimulus presentation, and the subject may be asked to view a visual fixation from the offset of the previous stimulus to the onset of the next stimulus. While the subject is cognitively engaging with the stimuli, an MR scanner repeatedly scans the subject’s brain to collect fMRI time series for making statistical inference about brain functions. A particular interest lies in the MR signal change following each stimulus. This signal change reflects the underlying brain activity evoked by the stimulus and is typically described by a function of time called hemodynamic response function (HRF). Studying the HRF helps us to understand the effects of the stimuli to the brain; see also Lazar (2008) and Lindquist (2008) for further details. Let {y(t ) | t = t1 , . . . , tN } be the fMRI time series acquired by the MR scanner from a brain voxel (volume unit) at regular time points t1 , . . . , tN . The first valid MR scan is at time t1 = 0 and (tn+1 − tn ) = τTR is a pre-specified constant time, n = 1, . . . , N. A popular model for analyzing the time series is (e.g. Zhang and Yu, 2008; Josephs et al., 1997) y(t ) =

Q   q =1

t

xq (t − τ )hq (τ )dτ + s(t ) + ε(t ),

(2)

0

where the indicator function xq (t ) is 1 if a qth-type stimulus appears at time t and is 0 otherwise; hq (t ) is the HRF evoked by the qth-type stimuli; the function s(t ) allows for a drift/trend of the time series; and ε(t ) is noise. One common objective is to estimate the HRF hq (t ). Hypothesis tests for identifying activated brain voxels may also be formulated based on the HRFs; see, e.g., Huettel (2012), Lazar (2008), and Zhang and Yu (2008). The heights of the HRF hq (t ) that contribute to an fMRI time series occur at regular, discrete time points (see also, Kao et al., 2012). As described in Zhang and Yu (2008), only the first few, say K (< N ), HRF heights are considered in analysis since the HRF typically comes back to baseline in a few seconds after the stimulus onset. Model (2) can thus be represented as y = Xh + S γ + ε.

(3)

M.-H. Kao / Statistics and Probability Letters 83 (2013) 1479–1483

1481

Here, y is an N-by-1 vector with elements yi = y(ti ). X = [X1 , . . . , XK ] is the N-by-QK , zero–one design matrix. Specifically, Xk = [x1,k , . . . , xQ ,k ] is the design matrix for the kth heights of the Q HRFs; k = 1, . . . , K . The QK -by-1 vector h = (hT1 , . . . , hTK )T represents the unknown heights of the discretized HRFs, where hk = (h1,k , . . . , hQ ,k )T and hq,k denotes the kth height of the discretized hq (t ). The vector S γ is the nuisance term with a specified S and an unknown parameter (vector) γ . The ith element of S γ corresponds to s(ti ) in Model (2). The N-by-1 vector ε represents the noise with elements εi = ε(ti ) and cov(ε) = 6. When deriving our results, we consider Model (3) with the following assumptions: (i) τISI = τTR ; (ii) S = jN is the N-by-1 vector of ones; and (iii) 6 = IN is the N-by-N identity matrix. For simplicity, we assume that var(εi ) = 1, but our results also hold when var(εi ) = σ 2 ̸= 1. A pre-scanning period before the first valid MR scan is also assumed; this pre-period is almost always needed in practice to allow the MR machine to reach a steady state. This allows us to consider a circular setting in which the last K − 1 stimuli of an fMRI design are presented to the subject in the pre-period; see Section 4 for a further discussion. Similar assumptions are also considered in the previous studies (e.g., Maus et al., 2011; Liu and Frank, 2004). When some of these assumptions are violated, one may resort to a computational method such as the algorithm proposed by Kao et al. (2009) for obtaining (near-)optimal designs. Under the previously described assumptions, an fMRI design can be represented as a sequence of length N. For a design d = {d1 , . . . , dN }, the nth element of the vector xq,1 described under Model (3) is

((xq,1 ))n =



1, 0,

if dn = q; otherwise;

q = 1, . . . , Q . Under a circular setting, the last K −1 elements dN −K +2 , . . . , dN are presented sequentially in the pre-scanning period. We thus have Xk = U k−1 X1 , where

 U =

0T I N −1

1 0



is a permutation matrix; k = 2, . . . , K − 1. For such a design d, the information matrix for h of Model (3) is Md = X T (IN − N −1 JN )X , where JN is the N-by-N matrix of ones. After some algebra, we have



C0

C1

C1T Md =  .

C0



..

CKT −1

..

. ···

··· .. . .. . C1T

CK − 1

.. .

C1 C0

    − JK ⊗ R ,  

(4)

where C0 = diag(n1 , . . . , nQ ) with nq = #{n | dn = q, n = 1, . . . , N } being the number of replicates of the qth-type stimuli (pq)

(pq)

in the design d; Cr = ((nr ))p,q=1,...,Q with nr = #{n | (dn−r , dn ) = (q, p), n = 1, . . . , N }, where dn−r is set to dN +n−r when (n − r ) ≤ 0, r = 1, . . . , K − 1; ⊗ is the Kronecker product; and R = ((np nq /N ))p,q=1,...,Q . 3. Results For fMRI experiments with one stimulus type (Q = 1), we consider the notion of universally optimal designs of Kiefer (1975). Among a class D of fMRI designs, a design d ∗ will be said to be universally optimal if Φ (Md∗ ) = mind∈D Φ (Md ) for every real-valued function Φ which is convex, permutation invariant (i.e., Φ (P T Md P ) = Φ (Md ) for any permutation matrix P), and nonincreasing (i.e., Φ (Md1 ) ≤ Φ (Md2 ) whenever Md1 − Md2 is nonnegative definite); see also Kiefer and Wynn (1984) and Sinha and Mukerjee (1982). This class of Φ -functions includes a wide range of commonly used optimality criteria. Proposition 1′ of Kiefer (1975) states that a design d ∗ with information matrix Md∗ is universally optimal if (i) Md∗ is a multiple of an identity matrix, and (ii) tr(Md∗ ) ≥ tr(Md ) for any other competing design d. This powerful tool allows us to develop Theorem 1, in which we set the design class D to contain all possible designs and use the term ‘circular designs’ to indicate designs to be used in a circular setting. Theorem 1. For Model (3) with Q = 1, τISI = τTR , S = jN , and uncorrelated noise, a circular design d = {d1 , . . . , dN } is universally optimal if (11)

(a) nr = (n1 )2 /N for all r = 1, . . . , K − 1; (b) n1 = N /2. Proof. With Q = 1, (a) and (b) imply that Md in (4) is a multiple of IK and tr(Md ) = K (n1 − (n1 )2 /N ) is maximized. The theorem thus follows from Proposition 1′ of Kiefer (1975).  For a binary m-sequence of length 2K − 1, each non-zero k-tuple appears 2K −k times and the zero k-tuple occurs 2K −k − 1 times for 1 ≤ k ≤ K ; see Proposition 5.2 of Golomb and Gong (2005). It can also be seen that, for an extended m-sequence, (11) each of the 2k k-tuple occurs 2K −k times, 1 ≤ k ≤ K . Consequently, n1 = 2K −1 = N /2 and nr = 2r −1 × 2K −(r +1) = 2K −2 . Conditions (a) and (b) of Theorem 1 are thus satisfied, and we have the following corollary.

1482

M.-H. Kao / Statistics and Probability Letters 83 (2013) 1479–1483

Corollary 1. Binary extended m-sequences of length 2K are universally optimal under the settings described in Theorem 1. For Q > 1, empirical results suggest that m-sequences and extended m-sequences are not A-optimal. Nevertheless, the previous studies indicate that these designs might be D-optimal although analytical results are unavailable; see also Maus et al. (2010). Here, we provide sufficient conditions for fMRI designs to be D-optimal. We then show that extended m-sequences satisfy these sufficient conditions. Theorem 2. For Model (3) with Q > 1, τISI = τTR , S = jT , and uncorrelated noise, a circular fMRI design d = {d1 , . . . , dN } is D-optimal if (pq)

(a) nr = np nq /N for all p, q = 1, . . . , Q , and r = 1, . . . , K − 1; (b) nq = N /(Q + 1), q = 1, . . . , Q . Proof. From (4) and the Fischer’s inequality, we have det(Md ) ≤ {det(C0 − R )}K . The equality holds when (Cr − R ) = O is a matrix of zeros for r = 1, . . . , K − 1. In addition, by using the fact that C0 − R = X1T (IN − N −1 JN )X1 and the formulas for the determinant of a partitioned matrix (e.g., Theorem 13.3.8 of Harville, 1997), one can show that

 det(C0 − R ) =

Q 

 nq

N−

q =1

Q  q=1

N

 nq

,

which is maximized when nq = N /(Q + 1) for q = 1, . . . , Q . Consequently, for a design d ∗ satisfying (a) and (b), we have det(Md∗ ) = max{[det(C0 − R )]K } ≥ det(Md ) for any other design d. That is, d ∗ is D-optimal.



We note that the conditions (a) and (b) in Theorem 2 do not guarantee a design maximizing tr(Md ) = K [tr(C0 − R )], and Proposition 1′ of Kiefer (1975) does not apply here. For an extended m-sequence of length N = (Q + 1)K , each k-tuple occurs (Q + 1)K −k times for 1 ≤ k ≤ K . It can thus be seen that the information matrix for these designs is (pq) IK ⊗ [(Q + 1)K −1 IQ − (Q + 1)K −2 JQ ], and nq = (Q + 1)K −1 = N /(Q + 1) and nr = (Q + 1)K −2 = np nq /N. We have the following result. Corollary 2. Extended m-sequences with elements 0, . . . , Q are D-optimal under the settings described in Theorem 2. 4. Discussion In this article, we provide theoretical results on the statistical optimality of extended m-sequences that are obtained by inserting a zero to a (K − 1)-tuple of zeros in an m-sequence of length (Q + 1)K − 1. To our knowledge, such results are currently unavailable, and our results help us to provide useful insights into fMRI experimental design problems. When deriving these results, we consider a circular setting where the last few elements of an fMRI design are presented in the prescanning period. Such a pre-period is typically required in real fMRI experiments. Nevertheless, we can always circularly shift the positions of the elements in an extended m-sequence so that the last few elements are all zeros. In this case, no stimuli are to be presented in the pre-period. The resulting designs remain optimal. In addition, under some assumptions, we are able to derive insightful analytical results. When some of these assumptions are violated, one may consider computer-generated designs such as those obtained by the approach of Kao et al. (2009). Nevertheless, developing similar theoretical results by relaxing one or more of these assumptions is still of great interest. Another direction of future research is to derive optimality results for designs with a length shorter than the extended m-sequences considered here. References Aguirre, G.K., Mattar, M.G., Magis-Weinberg, L., 2011. De bruijn cycles for neural decoding. NeuroImage 56, 1293–1300. Buracas, G.T., Boynton, G.M., 2002. Efficient design of event-related fMRI experiments using m-sequences. NeuroImage 16, 801–813. Golomb, S.W., Gong, G., 2005. Signal Design for Good Correlation for Wireless Communication, Cryptography, and Radar. Cambridge University Press, Cambridge, New York. Harville, D.A., 1997. Matrix Algebra from a Statistician’s Perspective. Springer, New York. Huettel, S.A., 2012. Event-related fMRI in cognition. NeuroImage 62, 1152–1156. Josephs, O., Turner, R., Friston, K., 1997. Event-related fMRI. Human Brain Mapping 5, 243–248. Kao, M.H., Mandal, A., Lazar, N., Stufken, J., 2009. Multi-objective optimal experimental designs for event-related fMRI studies. NeuroImage 44, 849–856. Kao, M.H., Mandal, A., Stufken, J., 2012. Constrained multi-objective designs for functional MRI experiments via a modified nondominated sorting genetic algorithm. Journal of the Royal Statistical Society: Series C-Applied Statistics 61, 515–534.

M.-H. Kao / Statistics and Probability Letters 83 (2013) 1479–1483

1483

Kiefer, J., 1975. Construction and optimality of generalized youden designs. In: Srivastava, J.N. (Ed.), A Survey of Statistical Designs and Linear Models. North-Holland, Amsterdam, pp. 333–353. Kiefer, J., Wynn, H.P., 1984. Optimum and minimax exact treatment designs for one-dimensional autoregressive error processes. The Annals of Statistics 12, 431–450. Lazar, N.A., 2008. The Statistical Analysis of Functional MRI Data. In: Statistics for Biology and Health, Springer, New York. Lidl, R., Niederreiter, H., 1994. Introduction to Finite Fields and their Applications, Rev. ed. Cambridge University Press, Cambridge, New York. Lindquist, M.A., 2008. The statistical analysis of fMRI data. Statistical Science 23, 439–464. Liu, T.T., Frank, L.R., 2004. Efficiency, power, and entropy in event-related fMRI with multiple trial types, part I: theory. NeuroImage 21, 387–400. MacWilliams, F.J., Sloane, N.J.A., 1977. The Theory of Error Correcting Codes. North-Holland Pub. Co., Amsterdam, New York. Maus, B., van Breukelen, G.J.P., Goebel, R., Berger, M.P.F., 2010. Robustness of optimal design of fMRI experiments with application of a genetic algorithm. NeuroImage 49, 2433–2443. Maus, B., van Breukelen, G.J.P., Goebel, R., Berger, M.P.F., 2011. Optimal design of multi-subject blocked fMRI experiments. NeuroImage 56, 1338–1352. Sinha, B.K., Mukerjee, R., 1982. A note on the universal optimality criterion for full rank models. Journal of Statistical Planning and Inference 7, 97–100. Zhang, C.M., Yu, T., 2008. Semiparametric detection of significant activation for brain fMRI. Annals of Statistics 36, 1693–1725.