Exact Dynamic Support Tracking with Multiple Measurement Vectors using Compressive MUSIC
aa r X i v : . [ c s . I T ] O c t Exact Dynamic Support Tracking withMultiple Measurement Vectors usingCompressive MUSIC
Jong Min Kim, Ok Kyun Lee and Jong Chul Ye
Abstract
Dynamic tracking of sparse targets has been one of the important topics in array signalprocessing. Recently, compressed sensing (CS) approaches have been extensively investigated asa new tool for this problem using partial support information obtained by exploiting temporalredundancy. However, most of these approaches are formulated under single measurement vectorcompressed sensing (SMV-CS) framework, where the performance guarantees are only in aprobabilistic manner. The main contribution of this paper is to allow deterministic tracking oftime varying supports with multiple measurement vectors (MMV) by exploiting multi-sensordiversity. In particular, we show that a novel compressive MUSIC (CS-MUSIC) algorithm withoptimized partial support selection not only allows removal of inaccurate portion of previoussupport estimation but also enables addition of newly emerged part of unknown support. Numericalresults confirm the theory.
Index Terms
Compressed sensing, joint sparsity, time varying signal, compressive MUSIC, optimized partialsupport selection
Correspondence to:Jong Chul Ye, Ph.D. Associate ProfessorDept. of Bio and Brain Engineering, KAIST373-1 Guseong-dong Yuseong-gu, Daejon 305-701, KoreaEmail: [email protected]: 82-42-350-4320Fax: 82-42-350-4310
I. I
NTRODUCTION
Dynamic target tracking problem that addresses the estimation of time varying support of movingtarget has been one of the important classical topics in array signal processing including radar,communication, and medical imaging applications. For example, in electroencephalography (EEG)or magnetoencephalography (MEG) source localization problems, it has been shown that theposition of the dipole moments during epileptic activities varies according to time and we areinterested in their spatio-temporal dynamics [1]. Dynamic MRI problem that tracks the motion ofhearts also belongs to this class of problem.Recently, there have been renewed interests for this problem with the help of a modern math-ematical tool called compressed sensing [2], [3]. These approaches try to exploit knowledges ofpartial support information obtained at the previous time point. More specifically, consider thefollowing time varying support estimation problem: min x ( t ) k x ( t ) k , subject to b ( t ) = A x ( t ) , t = 0 , , · · · , (1)where b ( t ) ∈ R m , and x ( t ) ∈ R n are noiseless measurement vector, and sparse signal at time t .Assuming that the support is assumed to change slowly, theoretical results [4] have demonstratedthat we can reduce the required sampling in compressed sensing reconstruction if we have partiallyknown support from the prior estimation results. For example, Vaswani and Lu proposed modified-CS algorithm [4] which addresses the exact reconstruction of noiseless case with partially knownsupport: min x ( t ) k ( x ( t )) I ( t − c k , subject to b ( t ) = A x ( t ) , t = 1 , , · · · , (2)where I ( t − is the previously estimated support, and ( x ( t )) I ( t − c denotes a subvector afterremoving the elements that correspond to the index set I ( t − . Suppose, furthermore, k = | supp x ( t ) | , u = | I ( t ) \ I ( t − | , and e = | I ( t − \ I ( t ) | . Then, if the restricted isometryconstant (RIP) for the sensing matrix A satisfies δ k + e + u < , (3)then the solution x ( t ) of Eq. (2) is the unique solution [4]. This is much weaker than ≤ δ k < for the original SMV-CS problem [5], in case of slowly time varying support with u ≪ k and e ≪ k . They further showed an l convex relaxation of Eq. (2) can provide the same l solutionof Eq. (2), if the following RIP condition is satisfied: δ u + δ u + δ k + e − u + δ k + e + 2 δ k + e + u < , (4)which is again relaxed sampling requirement than that of original CS problem δ k < √ − [5]. Therefore, exploiting the temporal redundancy has significant impact for reducing samplingrequirement for dynamic support tracking.Rather than solving the tracking problem Eq. (1) at each time, batch type approaches such asT-SBL (temporal sparse Bayesian learning) [1] collect the multiple snapshot data (for example, { b ( t ) } Nt =1 ) and process them together to estimate the dynamic varying support. Note that if thesupport changes slowly over time, then the resulting collection of problem becomes an multiplemeasurement vector problem. Accordingly, T-SBL converts the resulting MMV problem into ablock-sparse SMV problem, after which each block statistics are modeled using a specific Gaussianform temporal correlation structure. The update rule using the expectation-maximization (EM)method and its accelerated version can be then used to solve the resulting Bayesian problem [1].However, these approaches for dynamic support tracking is with SMV-CS framework and theirperformance guarantees is in a probabilistic sense. In practice, there are many situations wherewe can obtain multiple measurement vector information for time varying objects. For example, insingle-input multiple-output (SIMO) multiple access channel (MAC), multiple antenna can observelinear combination of individual codewords multiplied by the unknown channel gain from theindividual user [6]. In parallel MR cardiac imaging, multiple coils simultaneous obtain k-spacemeasurements of temporally varying hearts with distinct coil sensitivities. In EEG/MEG sourcelocalization problem, the dipole moments can be assumed relatively stationary during a short timewindow from which multiple snapshot of the sensor measurement can be obtained. All theseexamples acquire multiple measurement of the unknown signal vectors that share the same supportwith different weighting through identical sensing matrices.A fundamental question under this setup is what kind of diversity gain we can obtain overSMV-CS support tracking. To our knowledge, we are not aware of any prior investigation in thisregard. One of the main contributions of this paper is to show that a multiple measurement vector(MMV) framework not only extend the SMV counterpart, but also provides a unique advantageof “deterministic” support tracking for slow varying support estimation. Recall that MMV can measure multiple information of a set of vector that share the same sparsity pattern through theidentical sensing matrix. This paper shows that this joint sparsity pays off significantly in dynamicsupport tracking by relaxing probabilistic guarantee to a deterministic guarantee. The feasibility ofthe exact support tracking has significant impacts in practice.The breakthrough is based on our novel compressive multiple signal classification (CS-MUSIC)algorithm in MMV compressed sensing problem [7], in which a part of supports are foundprobabilistically using the conventional CS, after which the remaining supports are determineddeterministically using the generalized MUSIC criterion. In addition, CS-MUSIC allows us to findall k support as long as at least k − r + 1 support out of any k -support estimate are correct [8],where r denote the rank of the measurement matrix. This result provides an important clue fordeterministic and exact dynamic support tracking under MMV setup, in which the probabilisticcompressed sensing support estimation step is replaced by the support estimate from the previoussnapshots, after which the CS-MUSIC algorithm eliminates the incorrect portion of previous timepoint support estimation and then add newly updated support deterministically. This update schemeguarantees the exact support tracking in noiseless case under an appropriate sampling condition.Other contributions of our method include that the support error does not propagate along timedue to the self-correction step. Furthermore, using large system model, we can derive conditionswith which the proposed algorithm correct track the time varying support even in noisy cases.We believe that with these noticeable advantages of our algorithm we may find many importantapplications in radar, communication as well as biomedical application.This paper consist of following. Section II reviews the compressive MUSIC and support correc-tion criterion for MMV setup. In Section III, we derive our main theoretical results on samplingcondition for deterministic support tracking. Numerical results are given in Section IV, which isfollowed by conclusion in Section V. A. Notations and Mathematical Preliminaries
Throughout the paper, x i and x j correspond to the i -th row and the j -th column of matrix X ,respectively. When S is an index set, X S , A S corresponds to a submatrix collecting correspondingrows of X and columns of A , respectively. The following definitions are also used throughout thepaper. Definition 1: [9] The rows (or columns) in R n are in general position if any n collection ofrows (or columns) are linearly independent. Definition 2: [10]
Spark( A ) denotes the smallest number of linearly dependent columns of amatrix A . Definition 3 (Restricted Isometry Property (RIP)):
A sensing matrix A ∈ R m × n is said to havea k -restricted isometry property (RIP) if there exist left and right RIP constants < δ Lk , δ Rk < such that (1 − δ Lk ) k x k ≤ k A x k ≤ (1 + δ Rk ) k x k for all x ∈ R n such that k x k ≤ k . A single RIP constant δ k = max { δ Lk , δ Rk } is often referred toas the RIP constant.II. MMV C OMPRESSIVE S ENSING USING C OMPRESSIVE
MUSIC: A R
EVIEW
Let m , n and r be a positive integers ( m < n ) that represents the number of sensor elements,the ambient space dimension, and the number of snapshots, respectively. Suppose that we aregiven a multiple-measurement vector B ∈ R m × r , X = [ x , · · · , x r ] ∈ R n × r , and a sensing matrix A ∈ R m × n . A canonical form MMV problem [7] is given by the following optimization problem: minimize k X k (5) subject to B = AX, where k X k = | supp X | = k , supp X = { ≤ i ≤ n : x i = 0 } , and the measurement matrix B isfull rank, i.e. rank( B ) = r ≤ k X k .Recall that every MMV problem can be converted to a canonical form MMV using a singularvalue decomposition and dimension reduction as described in [7]. Now, We can easily expect thatthe diversity due to the joint sparsity can improve the recovery performance over SMV compressedsensing. Indeed, Chen and Huo [11], Feng and Bresler [12] and recently Davies and Elder [13]showed that X ∈ R n × r is the unique solution of AX = B if and only if k X k < spark( A ) + rank( B ) − ≤ spark( A ) − . (6)Note that we can expect rank( B ) / gains over SMV thanks to the MMV diversity. Furthermore,Feng and Bresler [12] showed that the noiseless l bound in Eq. (6) is achievable using MUSICalgorithm as long as r = rank( B ) = k . More specifically, suppose that the columns of a sensing matrix A ∈ R m × n are in general position. Then, according to [12], [14], for any j ∈ { , · · · , n } , j ∈ supp X if and only if Q ∗ a j = 0 , (7)where Q ∈ R m × ( m − r ) consists of orthonormal columns such that Q ∗ B = 0 so that R ( Q ) ⊥ = R ( B ) , which is often called “noise subspace”. Using the compressive sensing terminology, Eq. (7)implies that the recoverable sparsity level by MUSIC (with a probability for the noiselessmeasurement case) is given by k X k < m = spark( A ) − , (8)where the last equality comes from the definition of the spark . Therefore, the l bound (6) can beachieved by MUSIC bound in (8) when r = k [12].However, for any r < k , the MUSIC condition (7) does not hold. This is a major drawbackof MUSIC compared to CS algorithms that allow perfect reconstruction with a extremely largeprobability by increasing the sensor elements m . One the other hand, even thought the conventionalCS algorithms for MMV such as simultaneous OMP (S-OMP), p -thresholding [15], [16] have goodrecovery performance when r ≪ k , but they exhibit performance saturation as r increases andnever achieve the l bound with finite snapshot even in noiseless case. Recently, we showed thatthis drawback of the existing approaches can be overcome by the following generalized MUSICcriterion [7]. Theorem 1: [7] Assume that A ∈ R m × n , X ∈ R n × r , and B ∈ R m × r satisfy AX = B .Furthermore, we assume that k X k = k and A satisfies the RIP condition with the left RIP constant < δ L k − r +1 < . If we are given I k − r ⊂ supp X with | I k − r | = k − r and A I k − r ∈ R m × ( k − r ) ,which consists of columns whose indices are in I k − r , then for any j ∈ { , · · · , n } \ I k − r , a ∗ j h P R ( Q ) − P R ( P R ( Q ) A Ik − r ) i a j = 0 (9)if and only if j ∈ supp X .In [7], we demonstrate that the condition < δ L k − r +1 < for generalized MUSIC is equivalent to l bound (6), which implies that a computational expensive combinatorial optimization problem isnow reduced to | I k − r | support estimation from the original | I k | support estimation . Furthermore, When r = k , the condition (9) is the same as the MUSIC criterion (7) and no combinatorial algorithm is necessary. by Theorem 1, we can develop a computationally tractable relaxation algorithm called CompressiveMUSIC (CS-MUSIC) that relaxes the combinatorial optimization step of finding I k − r support usingthe conventional MMV-CS algorithms [7]. The algorithm can be stated as following: • (Step 1: compressed sensing step) Find k − r indices of supp X by any MMV compressivesensing algorithms such as 2-thresholding or SOMP. Let I k − r be set of selected indices and S = I k − r . • (Step 2: generalized MUSIC step) For j ∈ { , · · · , n } \ I k − r , calculate the quantities η ( j ) = a ∗ j [ P R ( Q ) − P R ( P R ( Q ) A Ik − r ) ] a j for all j / ∈ I k − r . Make an ascending ordering of η ( j ) , j / ∈ I k − r and choose indices that correspond to the first r elements and put these indices into S .In compressive MUSIC, we determine k − r indices of supp X with CS-based algorithms suchas 2-thresholding or S-OMP rather than l optimization, where the exact identification of k − r indices is a probabilistic matter. After that process, we recover remaining r indices of supp X witha generalized MUSIC criterion, which is given in Theorem 1, and this reconstruction process isdeterministic. This hybridization makes the compressive MUSIC applicable for all ranges of r ,outperforming all the existing methods. Similar observation have been made independently by Leeand Bresler [17] in their subspace augmented MUSIC (SA-MUSIC) algorithm.To analyze the performance of the compressive MUSIC, we should find the number of mea-surements with which we can identify the support of X . Due to the reduction of uncertainty from | I k | to | I k − r | , we can expect more relaxed sampling condition. In [7], we derived the samplingrequirements when subspace S-OMP or 2-thresholding is used as a compressed sensing step forcompressive MUSIC. The results can be summarized as following. The number of measurementsfor subspace S-OMP for partial support recovery exhibits two distinct characteristics dependingon the number of the measurement vectors. First, if the number of multiple measurement vectors r is sufficiently small, then the number of samples for S-OMP is reciprocally proportional to thenumber of multiple measurement vectors. On the other hand, we have sufficiently large number ofsnapshots such that lim n →∞ (log n ) /r is close to 0, then the number of measurements for S-OMPvaries from k to k according to the ratio of r and k so that the log n is not necessary. In particular,if the number of snapshots approaches the sparsity k , then we can identify the indices of supp X with only k measurements, which is equivalent to the required number of multiple measurementvectors for the success of conventional MUSIC. Furthermore, we demonstrated that the required SNR for the success of support recovery can be reduced and when the asymptotic ratio of thenumber of snapshots and the sparsity level (that is, lim n →∞ r/k ) is nonzero in the large systemlimit, only finite SNR is required, which is significant improvement over SMV-CS.In the original form of CS-MUSIC, the performance is, however, very dependent on the selectionof k − r correct indices of the support of X . In practice, even though the consecutive k − r stepsof S-OMP may not be correct, there are chances that among the estimates of k -sparse solution,part of the supports could be correct. Hence, if we have a mean to identify k − r correct supportin any order out of any k -sparse, then we can expect that the performance of the compressiveMUSIC will be improved. Of course, when (cid:0) kk − r (cid:1) is small, we may apply the exhaustive search,but if both k − r and r are not small, then the exhaustive search is hard to apply so that we have tofind some alternative method to identify the correct indices from the estimate of supp X . Indeed,the following support selection criterion can address the problem [8]. Theorem 2: [8] Assume that we have a canonical MMV model AX = B where A ∈ R m × n , X ∈ R n × r , k X k = k and r < k < m < n . If there is an index set I k ⊂ { , · · · , n } such that | I k | = min { k, spark( A ) − r } and | I k ∩ supp X | ≥ k − r + 1 , then for any j ∈ I k , j ∈ supp X ifand only if P Q k,j a j = , (10)where Q k,j is the orthogonal complement for R ([ B A I k \{ j } ]) , A I k \{ j } consists of columns of A whose index belongs to I k \{ j } and P ⊥ R ([ B A Ik \{ j } ]) is the orthogonal projection on R ([ B A I k \{ j } ]) ⊥ .In particular, if the columns of A are in general position, then we can take index set I k with | I k | = min { k, m − r + 1 } . Also, if A has an RIP condition with < δ k < , then we can take | I k | = k since r ≤ k .Theorem 2 informs us that we only require the success of partial support recover out of k -sparseestimate, rather than k − r consecutive correct CS step [7]. Accordingly, the compressive MUSICwith optimized partial support is then performed by following procedure. • [Step 1: compressed sensing] Estimate k indices of supp X by any MMV compressivesensing algorithm. Let I k be the set of indices which are taken in step 1. • [Step 2: support deletion] For j ∈ I k , calculate the quantities ζ ( j ) = k P Q k,j a j k . Makean ascending ordering of ζ ( j ) , j ∈ I k and choose indices that corresponds the first k − r elements and put these indices into S and remove the remaining ones. • [Step 3: support addition] For j ∈ { , · · · , n } \ S , calculate the quantities η ( j ) = a ∗ j [ P R ( Q ) − P R ( P R ( Q ) A Ik − r ) ] a j . Make an asending ordering of η ( j ) , j / ∈ S and choose indices that correspond to the first r elements and put these indices into S .The step in the above algorithm need not to be greedy so that we can also apply the convexoptimization algorithm such as l , minimization [18] or belief propagation [19].III. D ETERMINISTIC S UPPORT T RACKING USING C OMPRESSIVE
MUSIC
A. Noiseless Cases
In this section, we will show how the compressive MUSIC with optimized partial support canbe used for dynamic support tracking, whose joint support supp X ( t ) changes slowly along timeas illustrated in Fig.1. First, we define a canonical form of dynamic MMV problem. Definition 4:
A canonical form of noiseless dynamic MMV problem is given by set of MMVproblem with time varying k -sparse vectors X ( t ) ∈ R n × r that satisfies Y ( t ) = AX ( t ) as describedin following formulation: min X ( t ) k X ( t ) k , subject to B ( t ) = AX ( t ) , t = 0 , , · · · , (11)where supp X ( t ) = { ≤ i ≤ n : x ( t ) i = 0 } and | supp X ( t ) | = k ( t ) , the measurement matrix B ( t ) is full rank, i.e. rank( B ( t )) ≤ k ( t ) . Here we assume that rank( B ( t )) is constant so that welet r := rank( B ( t )) .Note that the canonical form MMV has the additional constraints that the measurement matrix isfull rank and rank( B ( t )) = r ≤ k ( t ) . This is not problematic since every dynamic MMV problemcan be converted into a canonical form using the following dimension reduction similar to [7] . • Suppose we are given the following linear sensor observations: B ( t ) = AX ( t ) where A ∈ R m × n and X ∈ R n × l satisfies k X ( t ) k = k ( t ) . • Compute the SVD as B ( t ) = U D r V ∗ , where D r is an r × r diagonal matrix, V ∈ C l × r consists of right singular vectors, and r = rank( B ) , respectively. • Reduce the dimension as B SV ( t ) = B ( t ) V and X SV ( t ) = X ( t ) V . • The resulting canonical form MMV becomes B SV ( t ) = AX SV ( t ) . Fig. 1. MMV problem for slowly time varying sparsity pattern.
We can easily show that rank( B SV ) = r ≤ k ( t ) and full rank and the sparsity k ( t ) := k X ( t ) k = k X SV ( t ) k with probability 1. Therefore, without loss of generality, the canonical form of dynamicMMV in Definition 4 is assumed throughout the paper.For such dynamic support tracking, we can apply our CS-MUSIC algorithm. However, if thenumber of snapshots is not sufficient, the amount of support estimation that need to be done byCS step is significantly larger than those recovered by the deterministic generalized MUSIC step.Since CS step allows the support recovery in a probabilistic sense, it is more prone to error; so weare interested in finding a deterministic algorithm that significantly outperform the existing one.The following Theorem 3 shows that if we have a correct estimation for the initial support I (0) of X (0) and the support changes are sufficiently small and the sparsity k ( t ) is fixed for all timepoint, then we can recursively identify the support of time-varying input signals in a deterministic manner. Theorem 3:
Suppose a noiseless canonical form of dynamic MMV problem satisfies | supp X ( t ) \ supp X ( t − | ≤ r − , (12)for all t = 1 , , · · · . Furthermore we assume that r ≤ k ( t ) ≤ k max for a positive integer k max and ≤ δ k max ( A ) < . Then, if we have a correct initial support estimation for X (0) , then we canidentify the correct support for all t > by applying the following procedure recursively: • [Initial support estimation] Let I ( t − be the support estimation of X ( t − ; • [Support deletion] Find an index set I ( t ) a ⊂ I ( t − such that I ( t ) a := { j ∈ I ( t −
1) : a ∗ j P Q ( t ) k,j a j = 0 } , where Q ( t ) k.j is the orthogonal complement for R [ B ( t ) , A I ( t − \{ j } ] ; • [Support addition] Find an index set I ( t ) such that I ( t ) = { j : a ∗ j [ P R ( Q ( t )) − P R ( P R ( Q ( t )) AI ( t ) a ) ] a j = 0 } , where Q ( t ) ∈ R m × ( m − r ) consists oforthonormal columns such that Q ( t ) ∗ B ( t ) = 0 ; • Set ˆ k ( t ) := | I ( t ) | be the sparsity estimate for X ( t ) and I ( t ) be the support estimate for X ( t ) . Proof:
See Appendix A.In Theorem 3, we assume the RIP condition ≤ δ L k max ( A ) < , instead of ≤ δ L k max − r +1 ( A ) < .If we assuming the RIP condition ≤ δ L k max − r +1 ( A ) < , when r > k max − k ( t ) / , we mayhave | I k | < k ( t ) . However, we can modify the support deletion procedure in Theorem 3 as thefollowing, under the condition | supp X ( t ) \ supp X ( t − | ≤ k max / . • [Support deletion] Find an index set I ( t ) a ⊂ I ( t − such that I ( t ) a := { j ∈ I ( t −
1) : a ∗ j P Q ( t ) k,j a j = 0 } , where Q ( t ) k.j is the orthogonal complement for R [ ˜ B ( t ) , A I ( t − \{ j } ] and ˜ B ( t ) consists of k max ] columns of B ( t ) . B. Noisy Cases
In practice, the measurements are noisy, so the theory we derived for noiseless measurementshould be modified. In the noisy case, when the sparsity are known a priori and does not changealong time, we can apply the following procedure. • Let t = 0 and let I (0) be the support estimation of X (0) . • For all t = 1 , , · · · , do – Let I ( t ) = ∅ . – For all j ∈ I ( t − , calculate the quantities ζ ( j ) = k P Q ( t ) j,k a j k . – Make an ascending ordering of ζ ( j ) and choose indices that correspond to the first k − r elements and put these indices into I ( t ) . – For j ∈ { , · · · , n }\ I ( t ) , calculate the quantities η ( j ) = a ∗ j h P R ( Q ( t )) − P R ( P R ( Q ( t )) AI ( t ) ) i a j . – Make an ascending ordering of η ( j ) , j / ∈ I ( t ) and choose indices that correspond to thefirst r indices and add these indices to I ( t ) . – I ( t ) is the estimation of supp X ( t ) and let t = t + 1 .However, if the sparsity changes along time, in the noisy cases, some of the steps in Theorem3 should be modified as follows: • [Support deletion] Set ǫ > and find an index set I ( t ) a such that I ( t ) a = { j ∈ I ( t −
1) : a ∗ j P Q ( t ) k,j a j < ǫ } where Q ( t ) k,j is the orthogonal complement for R [ Y ( t ) A I ( t ) \{ j } ] , where I ( t ) ⊂ I ( t − such that nrank[ Y ( t ) A I ( t ) ] = nrank[ Y ( t ) A I ( t − ] = r + | I ( t ) | , where nrank( A ) denotes the numerical rank of A . • [Support addition] Set ǫ > and find an index set I ( t ) b such that I ( t ) b = { j / ∈ I ( t ) a : a ∗ j P R ([ Y ( t ) A I t ) ]) ⊥ a j < ǫ } , where an index set I ( t ) ⊂ I ( t ) a such that nrank[ Y ( t ) A I ( t ) ] = nrank[ Y ( t ) A I ( t ) a ] = r + | I ( t ) | . In this section, we derive sufficient conditions for the threshold values and signal to noise ratiothat guarantee the correct identification of time varying support. For CS-MUSIC [7], we derived anexpression of SNR and the minimum number of sensor elements. Even though these derivation isbased on a large system model with a Gaussian sensing matrix, it has provided very useful insight.Therefore, we employed a large system model to derive a sufficient condition for the success ofproposed algorithm.
Definition 5:
A large system noisy canonical form of dynamic MMV is defined as an estimationproblem of k ( t ) -sparse vectors X ( t ) ∈ R n × r that shares a common sparsity pattern through multiplenoisy snapshots Y ( t ) = AX ( t ) + N ( t ) using the following formulation: minimize k X ( t ) k (13) subject to Y ( t ) = AX ( t ) + N ( t ) , where A ∈ R m × n is a random matrix with i.i.d. N (0 , /m ) entries, N = [ n , · · · , n r ] ∈ R m × r is an additive noise matrix, m → ∞ , k → ∞ as n → ∞ and rank( AX ( t )) = r ( t ) ≤ k ( t ) = k X ( t ) k . Here, we assume that ρ := lim n →∞ m/n > and γ = lim n →∞ k max /m > , α :=lim n →∞ r/k max ≥ exist and α ≤ − ǫ for some < ǫ < .Under the large system model, we have the following theorem. Theorem 4:
Consider the large system model dynamic MMV in Definition 5. Suppose a mini-mum
SNR satisfies
SNR min ( Y ( t )) := σ min ( B ( t )) k N k > κ ( B ( t )) + 1)1 − γ (1 + α ) , (14) where σ min ( B ( t )) is the minimum singular value for B ( t ) , k N k is the spectral norm of N ∈ R m × r and B ( t ) is the noiseless measurements, and α = lim n →∞ r/k max , γ = lim n →∞ k max /m . Then,for the noisy canonical form dynamic MMV problem for slowly time varying pattern that satisfiesEq. (12), the threshold values for support deletion and addition criterion to the correct partialsupport for X ( t ) are given by ǫ := (1 − γ (1 + α )) / , ǫ := (1 − γ ) / . (15) Proof:
See Appendix B. IV. NUMERICAL RESULTSThe first simulation is to demonstrate the performance of the proposed method to solve thetime varying MMV problem in Eq. (11) for different number of changes in supports at eachtime. We declared the algorithm as a success if the estimated support is the same as the true supp X , and the success rates were averaged for experiments. The simulation parameterswere as follows: m = 40 , n = 100 , r = 9 , and k ∈ { , , · · · , } , respectively. Elements ofsensing matrix A were generated by i.i.d. Gaussian random variable √ m N (0 , , and Gaussiannoise of SNR = 40 dB was added to each measurement vectors. At each time point, X ( t ) supp X ( t ) isgenerated by N (0 , . Fig.2 shows the recovery rates of time varying MMV problem using supporttracking method for t = 1 , , · · · , when the number of changed supports are , , , and ateach time point for Fig.2(a) ∼ (d), respectively. We used CS-MUSIC algorithm with S-OMP andthen applied optimized partial support selection at t = 1 , and time varying supports are estimatedby support tracking method recursively from t = 2 to t = 5 . In Fig.2, we can observe that theperformance gracefully decreases as the number of changes in supports increases. An interestingobservation is that the performance of the proposed method rather improves over time in Fig.2(a)and (b). However, the recovery ratio is getting lower but converges over time when the number ofchanges in supports is close to the upper bound r − for perfect recovery in noiseless case.Next, we applied the proposed algorithm to target tracking problem in D image and comparedit to MUSIC algorithm. The first row of Fig.3 indicates the original targets moving toward thedirection of red arrows over time. Each column (from left to right) indicates the sampled imageat t = 1 , , , and t = 41 , respectively. The simulation setting is the same with the previousone except m = 50 , n = 900 , t = 1 , , · · · , , and each target have a chance to move with Fig. 2. Recovery rates of time varying MMV problem using support tracking method when m = 40 , n = 100 , r = 9 ,SNR = 40 dB, and t = 1 , , · · · , . The number of changes in supports at each time point is (a) , (b) , (c) , and (d) . probability of r − k at each time point. The number of target k is . Here, we considered thenumber of measurement vectors is in the resting state, and used MUSIC algorithm to findsupports at t = 0 . The second and third row of Fig.3 indicate the results of support trackingmethod and MUSIC algorithm, respectively. Note that the proposed method successfully followsthe movement of original targets. V. CONCLUSIONThis paper expanded the sparse recovery with partially known supports in single measurementvector problem to multiple measurement vector problem with joint sparsity and proposed thesupport tracking algorithm to recover the slowly time varying supports. It is based on the recentlydeveloped compressive MUSIC algorithm with optimized partial support selection. The estimatedsupports at previous time can be used in optimized partial support selection to recover partialsupports at current time and it can be used in generalized MUSIC criterion to find remainingsupports. We also provided the maximum allowable number of changes in supports with support Fig. 3. The results of the target tracking problem in D image when m = 50 , n = 900 , k = 24 , and SNR = 40 dB. Weset r = 50 when t = 0 , and r = 9 for t > . The first row indicates the original targets moving toward the direction ofred arrows over time. The second and third row indicate the results of support tracking method and MUSIC algorithm,respectively. Each column (from left to right) indicates the sampled image at t = 1 , , , and t = 41 , respectively. tracking algorithm for exact reconstruction in noiseless case. Numerical results demonstrated thatthe proposed algorithm reliably reconstructs the time varying supports for various level of changesand successfully solves the target tracking problem in D image.A PPENDIX A Proof:
We only need to show that if we have a correct support for X ( t − , then we can alsoobtain a correct support estimation for X ( t ) by the support selection criterion and the generalizedMUSIC criterion. By the assumption, we have m ≥ k max ≥ k ( t −
1) + r so that if we have | supp X ( t ) ∩ I ( t − | ≥ k ( t ) − r + 1 , then by Theorem 2 we have for any j ∈ I ( t − , j ∈ supp X ( t ) if and only if a ∗ j P Q ( t ) k,j a j = 0 where Q ( t ) k,j is the orthogonal complement of R ([ Y ( t ) A I ( t − \{ j } ) . Since we have a noiselessMMV problem with slowly time varying pattern, we have | supp X ( t ) \ supp X ( t − | ≤ r − so that we have | supp X ( t ) ∩ I ( t − | ≥ k ( t ) − r + 1 and we can identify the correct partial supportof X ( t ) which has at least k ( t ) − r + 1 elements. Then, if we let I ( t ) a be the set of indices suchthat I ( t ) a = { j ∈ I ( t −
1) : a ∗ j P Q ( t ) k,j a j = 0 } , we have I ( t ) a ⊂ supp X and R ([ Y ( t ) A I ( t ) a ]) ⊂ R ( A supp X ( t ) ) . On the other hand, if we take aset I ( t, r ) ⊂ I ( t ) a ⊂ supp X such that | I ( t, r ) | = k ( t ) − r , we have R ([ Y ( t ) A I ( t ) a ]) ⊃ R ([ Y ( t ) A I ( t,r ) ]) = R ( A supp X ( t ) ) which implies R ([ Y ( t ) A I ( t ) a ) = R ( Y ( t ) A I ( t,r ) ) . Since ≤ δ k ( t ) − r +1 ( A ) ≤ δ k max − r +1 ( A ) < ,we can apply the generalized MUSIC criterion with I ( t, r ) ⊂ supp X where | I ( t, r ) | = k ( t ) − r .For j ∈ I ( t, r ) , we can easily see that a ∗ j h P R ( Q ( t )) − P R ( P R ( Q ( t )) AI ( t ) a ) i a j = a ∗ j P R ([ Y A I ( t ) a ]) ⊥ a j = a ∗ j P R ([ Y A I ( t,r ) ]) ⊥ a j = 0 . On the other hand, for j / ∈ I ( t, r ) , by the generalized MUSIC criterion, we have j ∈ supp X ( t ) ifand only if a ∗ j h P R ( Q ( t )) − P R ( P R ( Q ( t )) AI ( t ) a ) i a j = a ∗ j h P R ( Q ( t )) − P R ( P R ( Q ( t )) AI ( t,r ) ) i a j = 0 . Since I ( t, r ) ⊂ supp X , we have j ∈ supp X if and only if a ∗ j h P R ( Q ( t )) − P R ( P R ( Q ( t )) AI ( t ) a ) i a j = 0 . Hence, | I ( t ) | = k ( t ) and I ( t ) = supp X ( t ) .A PPENDIX B Proof:
Here, we let B ( t ) = AX ( t ) , σ min ( B ( t )) (or σ min ( B ( t )) ) be the minimum (or themaximum) nonzero singular value of B ( t ) . Then Y ( t ) = B ( t ) + N ( t ) is also of full column rankif k N ( t ) k < σ min ( B ( t )) . By [7], for such an N ( t ) , we have k P R ( Y ( t )) − P R ( B ( t )) k ≤ σ max ( B ( t )) + σ min ( B ( t ))] k N ( t ) k σ min ( B ( t ))( σ min ( B ( t )) − k N ( t ) k ) . (16)By the projection update rule, we have P R ([ B ( t ) A I t ) \{ j } ]) = P R ( A I t ) \{ j } ) + P R ( P ⊥ R ( AI t ) \{ j } ) B ( t )) (17) and P R ([ Y ( t ) A I t ) \{ j } ]) = P R ( A I t ) \{ j } ) + P R ( P ⊥ R ( AI t ) \{ j } ) Y ( t )) . (18)Since [ B ( t ) A I ( t ) \{ j } ] and [ Y ( t ) A I ( t ) \{ j } ] are of full column rank, by applying (17) and (18) asdone in [17], we have k P ⊥ R ([ B ( t ) A I t ) \{ j } ]) − P ⊥ R ([ Y ( t ) A I t ) \{ j } ]) k = k P R ( P ⊥ R ( AI t ) \{ j } ) B ( t )) − P R ( P ⊥ R ( AI t ) \{ j } ) Y ( t )) k≤ k P R ( B ( t )) − P R ( Y ( t )) k . (19)Then for any j ∈ I ( t ) \ supp X , we have a ∗ j P ⊥ R ([ Y ( t ) A I t ) \{ j } ]) a j = a ∗ j P ⊥ R ([ B ( t ) A I t ) \{ j } ]) a j + a ∗ j h P ⊥ R ([ Y ( t ) A I t ) \{ j } ]) − P ⊥ R ([ B ( t ) A I t ) \{ j } ]) i a j (20) ≥ min j / ∈ supp X a ∗ j P ⊥ R ([ B ( t ) A I t ) \{ j } ]) a j − max ≤ j ≤ n k a j k k P R ( Y ( t )) − P R ( B ( t )) k . Here, for each ≤ j ≤ n , m k a j k is a chi-square random variable with degree of freedom m sothat we have by Lemma 3 in [20], lim n →∞ max ≤ j ≤ n k a j k = 1 since lim n (log n ) /m = 0 . Further-more, for any j / ∈ supp X , a j is independent of P ⊥ R ([ Y ( t ) A I ( t − \{ j } ]) , so that m a ∗ j P ⊥ R ([ Y ( t ) A I t ) \{ j } ]) a j is a chi-squared random variable whose degree of freedom is at least m − k ( t ) − r + 1 since P ⊥ R ([ Y ( t ) A I t ) \{ j } ]) is a projection operator onto the orthogonal couplement of R ([ Y ( t ) A I ( t ) \{ j } ]) .Since lim n →∞ (log ( n − k ( t ))) / ( m − k ( t ) − r + 1) = 0 , again by Lemma 3 in [20], we have lim n →∞ min j / ∈ supp X m a ∗ j P ⊥ R ([ Y ( t ) A I t ) \{ j } ]) a j m − k ( t ) − r + 1 ≥ so that lim n →∞ min j / ∈ supp X a ∗ j P ⊥ R ([ Y ( t ) A I t ) \{ j } ]) a j ≥ − γ (1 + α ) (21)since k ( t ) ≤ k max for all t = 0 , , · · · . On the other hand, if we use the definition of SNR min ( Y ( t )) and the definition of the condition number of B ( t ) on (16), i.e. κ ( B ( t )) = ( σ max ( B ( t ))) / ( σ min ( B ( t ))) ,we have k P R ( Y ( t )) − P R ( B ( t )) k ≤ κ ( B ( t )) + 1) SNR min ( Y ( t )) − < − γ (1 + α )2 , (22)by the condition (14). Combining (20), (21) and (22), we have for any j ∈ I ( t − \ supp X , wehave a ∗ j P ⊥ R ([ Y ( t ) A I t ) \{ j } ]) a j > − γ (1 + α )2 . On the other hand, for j ∈ I ( t − ∩ supp X ( t ) , we have a ∗ j P ⊥ R ([ B ( t ) A I t ) \{ j } ]) a j = 0 by the supportselection criterion. Then, by the similar reasoning as above, we have for any j ∈ I ( t ) ∩ supp X ( t ) ,we have a ∗ j P ⊥ R ([ Y ( t ) A I t ) \{ j } ]) a j < − γ (1 + α )2 . This completes the proof for the threshold values for support selection criterion. The proof forthe threshold values for generalized MUSIC are the same except that m a ∗ j P ⊥ R ([ Y ( t ) A I t ) ]) a j is achi-squared random variable whose degree of freedom is m − k ( t ) for j / ∈ supp X ( t ) .A CKNOWLEDGMENT
This work was supported by the Korea Science and Engineering Foundation (KOSEF) grantfunded by the Korea government (MEST) (No.2010-0000855).R
EFERENCES [1] Z. Zhang and B.D. Rao, “Sparse signal recovery with temporally correlated source vectors using joint sparsebayesian learning,”
IEEE J. of Selected Topics in Signal Processing , vol. 5, no. 5, pp. 912–926, 2011.[2] D. L. Donoho, “Compressed sensing,”
IEEE Trans. on Information Theory , vol. 52, no. 4, pp. 1289–1306, April2006.[3] E. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highlyincomplete frequency information,”
IEEE Trans. on Information Theory , vol. 52, no. 2, pp. 489–509, Feb. 2006.[4] N. Vaswani and W. Lu, “Modified-CS: modifying compressive sensing for problems with partially known support,”
IEEE Trans. Signal Process. , vol. 58, no. 9, pp. 4595–4607, 2010.[5] E. Candes and T. Tao, “Decoding by linear programming,”
IEEE Trans. on Information Theory , vol. 51, no. 12,pp. 4203–4215, Dec. 2005.[6] Y. Jin and B.D. Rao, “Support recovery of sparse signals in the presence of multiple measurement vectors,” arXivpreprint , 2011, http://arxiv.org/pdf/1109.1895.[7] J. M. Kim, O. K. Lee, and J. C. Ye, “Compressive MUSIC: a missing link between compressive sensing and arraysignal processing,” to appear in IEEE Trans. Inf. Theory , 2011.[8] J. M. Kim, O. K. Lee, and J. C. Ye, “Compressive MUSIC with optimized partial support for joint sparse recovery,” in Proc. IEEE Int. Symp. Inf. Theory (ISIT) , 2011.[9] D. L. Donoho, “Neighborly polytopes and sparse solution of underdetermined linear equations,”
Tech. report,Department of Statistics, Stanford University , 2005.[10] D. L. Donoho and M. Elad, “Optimally sparse representation in general (non-orthogonal) dictionaries via l minimization,” Proceedings of the National Academy of Sciences of the United States of America , vol. 100, no. 5,pp. 2197–2202, 2003.[11] J. Chen and X. Huo, “Theoretical results on sparse representations of multiple measurement vectors,”
IEEE Trans.on Signal Processing , vol. 54, no. 12, pp. 4634–4643, 2006. [12] P. Feng, Universal minimum-rate sampling and spectrum-blind reconstruction for multiband signals , Ph.D.dissertation, University of Illinois, Urbana-Champaign, 1997.[13] M. E. Davies and Y. C. Eldar, “Rank awareness for joint sparse recovery,” preprint , 2010,http://arxiv.org/PS cache/arxiv/pdf/1004/1004.4529v1.pdf.[14] R. Schmidt, “Multiple emitter location and signal parameter estimation,”
IEEE Trans. on Antennas and Propagation ,vol. 34, no. 3, pp. 276–280, 1986.[15] R. Gribonval, H. Rauhut, K. Schnass, and P. Vandergheynst, “Atoms of all channels, unite! Average case analysisof multi-channel sparse recovery using greedy algorithms,”
Journal of Fourier Analysis and Applications , vol. 14,no. 5, pp. 655–687, 2008.[16] Y.C. Eldar and H. Rauhut, “Average case anlysis of multichannel sparse recovery using convex relaxation,”
IEEETrans. on Information Theory , vol. 56, pp. 505–519, 2010.[17] K. Lee and Y. Bresler, “Subspace-augmented music for joint sparse recovery,” preprint , 2010,http://arxiv.org/PS cache/arxiv/pdf/1004/1004.3071.pdf.[18] D. Malioutov, M. Cetin, and AS Willsky, “A sparse signal reconstruction perspective for source localization withsensor arrays,”
IEEE Trans. on Signal Processing , vol. 53, no. 8, pp. 3010–3022, 2005.[19] J.M. Kim, W.H. Chang, B.C. Jung, D. Baron, and J.C. Ye, “Belief propagation for joint sparse recovery,” arXivpreprint , 2011, http://arxiv.org/pdf/1102.3289.[20] S. Rangan A.K. Fletcher and V.K. Goyal, “Necessary and sufficient conditions for sparsity pattern recovery,”