Efficient Tracking of Sparse Signals via an Earth Movers Distance Dynamics Regularizer

1 Ef ficient T racking of Sparse Signals via an Earth Mo ver’ s Distance Dynamics Re gularizer Nicholas P . Bertrand ∗ , Adam S. Charles ∗ , Member , IEEE, John Lee ∗ , Pa vel B. Dunn, Christopher J. Rozell Senior Member , IEEE Abstract —T racking algorithms such as the Kalman filter aim to impro ve infer ence performance by le veraging the temporal dynamics in streaming observations. Howe ver , the tracking reg- ularizers are often based on the ` p -norm which cannot account for important geometrical relationships between neighboring signal elements. W e propose a practical approach to using the earth mover’ s distance (EMD) via the earth mover’ s distance dynamic filtering (EMD-DF) algorithm f or causally tracking time-varying sparse signals when there is a natural geometry to the coefficient space that should be respected (e.g ., meaningful ordering). Specifically , this paper presents a new Beckmann formulation that dramatically reduces computational complexity , as well as an evaluation of the perf ormance and complexity of the proposed approach in imaging and frequency tracking applications with real and simulated neurophysiology data. I . I N T R O D U C T IO N T racking algorithms aim to improve the performance of statistical inference procedures for time series by incorporating information from a dynamics model that describes how the signal ev olves. W e consider the linear observ ation model y n = A n x n + σ  n , (1) where for each time step n , x n is the underlying signal, A n is a linear observ ation operator , σ  n is Gaussian measurement noise with variance σ 2 , and y n is the resulting measurement vector . The signal e volv es according to a dynamics function g : x n +1 = g n ( x n ) + η n , (2) where η n is a noise vector called the innovations that accounts for inaccurate modeling of the dynamics. When g is linear and the signal, observ ation noise and innov ations are Gaussian, the classical Kalman filter provides an efficient way to compute the optimal (i.e., minimum expected ` 2 error) estimate taking into account all measurements up to the current time [ 1 ]. The estimate produced by Kalman filtering may be expressed as b x n = argmin x k y n − A n x k 2 B + k x − G n − 1 b x n − 1 k 2 C , (3) ∗ Equal contributions. N. P . Bertrand, J. Lee, P . B. Dunn, and C. J. Rozell are with the School of Electrical and Computer Engineering, Georgia Institute of T echnology , Atlanta, GA 30332-0250 USA (email: nbertrand@gatech.edu; john.lee@gatech.edu; pav eldunn@gatech.edu; crozell@gatech.edu). A. S. Charles is with the Princeton Neuroscience Institute, Princeton University , Princeton, NJ 08540 USA (email: adamsc@princeton.edu). This work was partially supported by NSF grant CCF-1409422, James S. McDonnell Foundation grant number 220020399, NIH NRSA Training Grant in Quantitativ e Neuroscience number T32MH065214 and the DSO National Laboratories of Singapore. where G n is the matrix representing the dynamics operator g n , b x n − 1 is the previous signal estimate, and k·k 2 B and k·k 2 C denote Mahalanobis norms weighted appropriately by the cov ariance matrices of the noise, innov ations, and previous signal estimate. The Kalman filter and its extensions [ 2 ] hav e been used e xhausti vely in a multitude of scientific and engineering applications. In addition to these classic models, non-Gaussian sparsity models hav e become increasingly popular due to their state-of- the-art performance in a v ariety of problems (e.g., in image processing [ 3 ] and compressi ve sensing [ 4 ]). Sparse inference problems with static data vectors are well-studied, resulting in many algorithmic advances and performance guarantees [ 4 ], [ 5 ]. In the spirit of the Kalman filter , sparse tracking algorithms hav e sho wn utility for dynamic filtering with time-varying sparse signals [ 6 ]–[ 13 ]. Howe ver , in many applications with discretized domains, commonly used pointwise dynamics regularizers (e.g., the ` p -norm) disproportionately penalize predictions with slight mismatch in the signal support because they do not incorporate knowledge of meaningful geometry (when it exists) into the penalty . Consider, for example, an imaging scenario where we wish to track a single-pixel target moving through a scene. An ` p -norm based regularizer assigns equal penalties to any prediction in which the target is not precisely in the correct support location regardless of the distance between the erroneous pixel and the true position. Similarly , when tracking time varying frequencies, the ordering of the frequencies in the discrete Fourier transform (DFT) matrix results in a geometric relationship among the DFT coefficients which is not ef fectiv ely utilized with ` p -norm regularizers. Although one may consider tracking in coordinate space instead of pixel space, this approach scales poorly in the number of targets. In this work, we propose a practical approach to using the earth mover’ s distance (EMD) via the earth mover’ s distance dynamic filtering (EMD-DF) algorithm for causally tracking time-v arying sparse signals when there is a natural geometry to the coef ficient space that should be respected (e.g., meaningful ordering). Specifically , this paper presents a new Beckmann formulation that dramatically reduces computational comple xity , as well as an e v aluation of the performance and complexity of the proposed approach in imaging and frequency tracking applications with real and simulated neurophysiology data. The nov el algorithmic formulation, performance characterization, and ev aluation on real data introduced in this work represent crucial dev elopments in the practicality of our approach beyond the preliminary explorations presented in [ 14 ], [ 15 ]. 2 I I . B AC K G RO U N D A. Sparse Dynamic F iltering A vector x ∈ C N is said to be sparse if only a few of its elements are non-zero. Suppose y contains noisy observ ations of x through a linear measurement operator A ∈ C M × N . For example, results in the compressed sensing literature show that under certain conditions on A , x may be recovered from y e ven when M  N . Of the many sparse in verse algorithms that exist (e.g., [ 5 ], [ 16 ]), one popular optimization-based approach is Basis-Pursuit Denoising (BPDN): b x = argmin x 1 2 k y − Ax k 2 2 + λ k x k 1 . (4) Recent work has also extended these ideas for static sparse recovery to tracking algorithms for sparse time-varying signals. Early w ork in this area included batch (i.e., non- causal) approaches [ 17 ], [ 18 ] and modifications to the causal Kalman filter [ 19 ]. More recent causal approaches include Basis Pursuit Denoising Dynamic Filtering (BPDN-DF) which provides theoretical con ver gence guarantees, and Reweighted- ` 1 Dynamic Filtering (R WL1-DF) which was found to be more robust to model mismatch [ 12 ]. BPDN-DF modifies standard BPDN with the addition of a tracking regularizer: b x n = argmin x 1 2 k y n − Ax k 2 2 + λ k x k 1 + γ k x − e x n k 2 2 , (5) where e x n = g ( b x n − 1 ) is the prediction produced using the dynamics function g . This additional term encourages solutions which adhere to the dynamics model. Similarly , R WL1-DF modifies R WL1 by injecting dynamics into the recov ery process via second order statistics. While both BPDN-DF and R WL1- DF can improve performance in the recov ery of time varying signals, each of these algorithms injects dynamics information in a point-wise f ashion and thus fail to capture the geometric relationship between neighboring signal elements. B. Earth Mo ver’ s Distance The earth mo ver’ s distance (EMD) is a metric which gre w out of the optimal transport (O T) literature initiated by Monge [ 20 ]. The EMD has recently been increasingly used in a variety of applications such as image and histogram comparison [ 21 ], [ 22 ], as well as for sparse inv erse problems [ 23 ]–[ 25 ]. Intuitively , if we visualize the first signal as being composed of piles of dirt and the second as holes, the EMD computes the minimum amount of work needed to fill the holes with dirt 1 . A key property to note is that the EMD is inherently aw are of the geometric relationship between signal elements via a user- defined distance matrix ( R ij ) which describes the cost to transport mass along the signal support. This is in stark contrast to ` p metrics, and is the primary moti vation for its use as a tracking regularizer . The traditional EMD formulation in volves solving for O  N 2  flow variables, which has the potential to be com- putationally prohibiti ve for large problems. For applications 1 W e refer the reader to [ 15 ] for mathematical details of the traditional EMD formulation. where the cost matrix R represents Euclidean distances (e.g., video with uniform ly gridded pixels), geometric structure can be exploited to also reduce the optimization variable complexity in exact EMD solutions from O  N 2  to O ( N ) via the Beckmann pr oblem [ 26 ], [ 27 ]. This formulation poses the EMD as a minimum flux problem of a fluid flowing between a source and a sink (i.e., the input arguments of the EMD problem): d emd ( x , y ) = min M k M k 2 , 1 subject to div( M ) + y − x = 0 , (6) where the diver gence operator is defined as div( M )[ i, j ] = ( M x [ i, j ] − M x [ i − 1 , j ]) + ( M y [ i, j ] − M y [ i, j − 1]) , (7) the ro ws of M contain points in a D -dimensional vector field, k M k 2 , 1 := P N i =1 k m i k 2 denotes the sum of their Euclidean norms and zero-flux boundary conditions are enforced (i.e., M [ i, j ] = 0 whenev er i or j falls outside the support). I I I . E A RT H M OV E R ’ S D I S TA N C E D Y N A M I C F I LT E R I N G In earth mover’ s distance dynamic filtering (EMD-DF), the causal estimate of the signal at time n is gi ven by: b x n = argmin x 1 2 k y n − Ax k 2 2 + λ k x k 1 + γ d emd ( x , e x n ) . (8) EMD-DF has a similar structural form as BPDN-DF at first glance [ 12 ], [ 14 ], [ 15 ], b ut the use of an EMD penalty instead of an ` 2 dynamics regularizer is non-trivial because the ev aluation of the EMD itself requires the solution of an optimization program. Using the traditional EMD formulation for general cost distances, the EMD-DF optimization program ( 8 ) in volves solving N signal v ariables and an additional N 2 flow variables. Thus, EMD regularization incurs a potentially prohibitiv e increase in computational complexity compared to algorithms such as BPDN or R WL1. Howe ver , in the common case when the transport costs are Euclidean (i.e., R ij = k z i − z j k 2 2 where z i and z j are the support locations of signal elements x i and x j ), we can exploit Beckmann’ s formulation of the EMD ( 6 ) to reduce the number of EMD v ariables from O  N 2  to O ( N ) . This critical reduction in computational complexity enables pre- viously intractable applications. Preliminary explorations [ 14 ], [ 15 ] have pre viously addressed traditional EMD limitations to account for signals that are signed or complex valued. The Beckmann EMD requires that the signals have unit mass (i.e., k x k 1 = k y k 1 = 1 ), meaning that we cannot simply apply the pre-existing method. In the following, we outline how a reformulation of the Beckmann problem for unequal total masses [ 28 ] may be incorporated into the EMD-DF program. T o allow input arguments with unequal total mass, we introduce slack variables w , v to bound the flux from the original source x and sink y . The modified EMD program is then: d emd ( x , y ) = min M , w , v k M k 2 , 1 subject to div( M ) + v − w = 0 , 0 ≤ w ≤ x , 0 ≤ v ≤ y , k w k 1 = k v k 1 = min( k x k 1 , k y k 1 ) , (9) 3 where w , v are nonnegati ve v ectors with the same dimensions as x , y . This optimization searches for the minimal vector field configuration that describes, via the first constraint, its flux to be tra veling between a source w and a sink v . The second constraint describes the source and sink as nonne gati ve slack variables that are bounded above by their proxies x and y respectiv ely; this constraint is analogous to the mass preservation constraints in the traditional EMD formulation. The last constraint states that the induced flux must be bounded by the total mass of the smaller operand signal. This formulation has N ( D + 2) variables, where D is the dimensions of the vector field (e.g., D = 2 for images). Applying this EMD formulation and employing the same strategy as [ 15 ] to replace the min term with a slack variable u , ( 8 ) becomes b x n = argmin x , M ,u, v , e v 1 2 k y n − Ax k 2 2 + λ k x k 1 + γ k M k 2 , 1 − µu subject to div( M ) + e v − v = 0 , 0 ≤ v ≤ x , 0 ≤ e v ≤ e x , k v k 1 = k e v k 1 = u, u ≤ k x k 1 , u ≤ k e x k 1 . (10) The complex variant of EMD-DF in [ 15 ] can also be trivially con verted to adopt this formulation, though it is not shown here for the sake of brevity . As we mentioned in Section II-B , ( 10 ) enjoys a reduction in v ariable complexity from O  N 2  to O ( N ) while preserving the benefits of partial O T (in contrast to the traditional balanced OT Beckmann formulation) and av oiding the approximation error associated with methods such as Sinkhorn iterations. Finally , we note that other recent works [ 29 ]–[ 31 ] incorporate optimal transport regularizers in in verse problems using the Sinkhorn algorithm [ 32 ], [ 33 ]. Howe ver , the work presented here is distinct in two subtle but important ways. First, our proposed partial Beckmann formulation provides an alternativ e numerical approach that offers attractiv e linear v ariable com- plexity (similar to Sinkhorn methods) but without sacrificing accuracy . Sinkhorn approaches use entropic regularization to trade off accuracy vs. speed, limiting their utility for finding sparse solutions ov er time due to mass diffusion across neighboring support. In contrast, the Beckmann formulation reflects the true optimal transport distance (subject only to negligible discretization errors). Second, our partial transport formulation results in a linearly constrained quadratic program that is easily implementable with of f-the-shelf solvers (e.g., CVX, Gurobi, Mosek). In contrast, an unbalanced Sinkhorn approach would necessitate a custom solver (e.g., alternating optimization [ 34 ], [ 35 ]) that would be non-tri vial to implement. I V . R E S U L T S W e demonstrate the utility and performance of EMD-DF through a series of simulations on synthetic and real data. First, we study the problem of tracking time v arying frequencies in a 1-D time series of real neurophysiology data. Next, we consider the problem of tracking a wav efront in synthetically generated data motiv ated by the phenomenon of trav eling wa ves which appear in electrophysiology data. Finally , we demonstrate the significant numerical speed up of EMD-DF due to Beckmann’ s formulation. Throughout these simulations, we use the CVX software package [ 36 ] which employs interior point methods to carry out the EMD-DF optimization. Hyperparameters are tuned manually for real data, while direct search [ 37 ] is used for synthetic data where ground truth is known. A. T racking Neural Oscillations In this section, we apply EMD-DF to the problem of spectrum estimation in neurophysiology recordings. Oscillatory behavior is prominent in a variety of neural recording settings and there is great interest in the neuroscience community to understand the functional role of these oscillations [ 38 ], [ 39 ]. In many studies, the tools used for spectral analysis of neural recordings are based on the classical short-time Fourier transform (STFT). The time and frequency resolution of such techniques is thus limited by the uncertainty principle which prev ents simultaneously achieving high frequency and time resolution. Here, we study how higher time-frequency (TF) resolution may be obtained by imposing a sparsity model on the data and using EMD-DF for recovery in an overcomplete DFT dictionary . EMD-DF may be used as a causal alternativ e to well- established spectral sharpening methods [ 40 ]–[ 43 ] which also aim to improve the resolution of time-frequency representations, but are restricted to operate on batch data [ 15 ]. Causal algorithms are crucial in online applications such as closed-loop control. Here we employ EMD-DF to estimate the spectrum in a segment of real electrophysiology data recorded from a tetrode in rat hippocampus [ 44 ]. 2 W e take the measurement matrix A to be a 5 times o vercomplete DFT matrix and track the top two frequencies by setting the remaining frequencies in the prediction to zero via the dynamics function g . Figure 1 shows TF plots produced by the spectrogram and EMD-DF . Because EMD-DF utilizes the overcomplete DFT matrix for recov ery , it produces a TF plot with vastly improved frequency resolution. Additionally , the spectrogram suffers from se vere leakage in the 5–10Hz frequenc y band, an artifact which is not present in the sparse TF representation. Finally , the improv ed resolution of the sparse TF plot rev eals more subtle oscillatory dynamics that cannot be observed in the spectrogram. B. T racking T raveling W aves T ra veling wav es are another form of neural oscillation pattern of interest in the neuroscience community . For example, wa ve propagation has been sho wn to correlate to events and performance in tasks in volving neurosurgical patients [ 45 ]. W e generate synthetic traveling wav e data using the phase-coupled Kuramoto oscillator model which has been used to study the properties of traveling wav es in neuronal activity [ 46 ], [ 47 ]. The Kuramoto model describes the instantaneous phase of each node in an array of linked oscillators via a system of differential equations. Motiv ated by the work in [ 46 ], we simulate a 40 × 40 oscillator array where the instantaneous phase of oscillator 2 The authors would like to thank C. Kemere for the tetrode recording data. 4 Fig. 1. T ime-frequency plots for a single channel of tetrode data recorded from the rat hippocampus. Data is sampled at 250 Hz and an analysis window length of 72 samples is used for both plots. The spectrogram (left), which is produced using the traditional STFT with a hamming window , yields lower frequency resolution and severe leakage in the lower frequencies. The TF plot on the right is produced by EMD-DF with an 5 times overcomplete DFT matrix, resulting in high enough frequency resolution to smoothly track subtle changes in frequency . ( i, j ) with intrinsic frequency ω ij depends on its four nearest neighbors via the equation θ 0 ij = ω ij + 300 P 4 k =1 sin( θ k ) . W e choose ω ij = 2 + 0 . 103 i + 0 . 359 j to produce a linear frequency gradient which has been observed to yield traveling wa ve solutions. After solving for the θ ij , we threshold the oscillator voltage sin( θ ij ( t )) to extract the wa vefront x ij ( t ) . The goal of these simulations is to reco ver the wav efront from noisy linear measurements y ( t ) = Φ x ( t ) +  . W e measure upper bound tracking performance by providing the ground truth pre vious frame as the prediction. Figure 2 shows how EMD-DF enables nearly perfect recovery compared to BPDN and BPDN-DF which are unable to ef fecti vely use the prediction. C. Computational Scalability Finally , we ev aluate the runtime of Beckmann EMD-DF compared to previously dev eloped versions which use the traditional EMD formulation with general ground costs R ij . W e conduct a stylized tar get tracking simulation in which a sparse collection of targets mov e to adjacent support locations with equal probability . W e take random Gaussian measure- ments and scale the problem between state sizes of 12 × 12 ( N = 144) and 48 × 48 ( N = 2304) . For each state size, the sparsity level is fixed at 5% and 10 trials are run on a personal computer with a 3.5 Ghz Intel Core i7 processor . The Beckmann formulation yields solutions significantly faster than the general EMD formulation, especially for large problem sizes (Figure 3 ). W e also note that the discrepancy between solutions obtained using each method (measured in root mean squared error: k x − y k 2 2 / √ N ) is ne gligible — on the order of 10 − 3 throughout the entire range of problem sizes. V . S U M M A RY EMD-DF can be a more effecti ve dynamics regularizer for sparse dynamic filtering compared to traditional point- wise methods (e.g., the ` p -norm), but potentially suffers from computational complexity that limits problem sizes. The results presented here demonstrate that EMD-DF is a practical sparse tracking method when used with a Beckmann formulation that scales linearly in the number of optimization variables Fig. 2. Single step recov ery of wavefronts from a 40 × 40 Kuramoto oscillator array via linear Gaussian measurements. (T op) A linear gradient in the oscillators’ intrinsic frequencies results in traveling waves across the array . (Middle) Examples of wav efronts recovered using several methods. BPDN produces poor recovery performance due to undersampled ( M / N = 0 . 15 ) and noisy ( σ = 0 . 08 ) measurements. Dynamics regularization in BPDN-DF is of little help due to its inability to effectiv ely utilize the prediction which has little support overlap with the ground truth signal. EMD regularization is robust to the mismatch in pixel location between the ground truth and prediction and thus enables successful recovery of the wavefront. (Bottom row) Recovery performance averaged over 10 trials. Error bars indicate α = 0 . 01 confidence intervals. (Bottom-left) EMD-DF produces superior performance for various values of the noise standard deviation σ . (Bottom-right) The ` 2 - norm dynamics regularization in BPDN-DF actually degrades performance compared to BPDN for moving targets. In contrast, EMD-DF is substantially more robust to support location mismatch caused by target movement. Fig. 3. Evaluation of computational speed up. W e compare the runtime and the difference in solutions for two formulations of EMD-DF: EMD-DF (General) which adopts generic distance costs, and EMD-DF (Beckmann) which assumes Euclidean distance costs. The left plot demonstrates that EMD- DF (Beckmann) significantly outperforms EMD-DF (General) in runtime, and in the right plot, the differences in solutions are shown to be negligible. instead of quadratically , enabling substantial performance gains for large-scale problems. Finally , we demonstrate the utility of EMD-DF by tracking wav efronts in a Kuramoto oscillator network and by tracking sparse frequencies in real electrophysiology data. 5 R E F E R E N C E S [1] R. E. Kalman, “ A ne w approach to linear filtering and prediction problems, ” Journal of Basic Engineering , v ol. 82, no. 1, pp. 35–45, Mar . 1960. [2] S. S. Haykin, Kalman Filtering and Neural Networks . New Y ork: W iley, 2001, oCLC: 52366672. [3] M. Elad, M. A. T . Figueiredo, and Y . Ma, “On the role of sparse and redundant representations in image processing, ” Pr oceedings of the IEEE , vol. 98, no. 6, pp. 972–982, Jun. 2010. [4] R. G. Baraniuk, “Compressi ve sensing [lecture notes], ” IEEE Signal Pr ocessing Magazine , vol. 24, no. 4, pp. 118–121, Jul. 2007. [5] J. A. T ropp and A. C. Gilbert, “Signal recovery from random mea- surements via orthogonal matching pursuit, ” IEEE T ransactions on Information Theory , vol. 53, no. 12, pp. 4655–4666, Dec. 2007. [6] M. S. Asif and J. Romberg, “Dynamic updating for ` 1 minimization, ” IEEE Journal of Selected T opics in Signal Pr ocessing , vol. 4, no. 2, pp. 421–434, Apr . 2010. [7] M. S. Asif, A. Charles, J. Romberg, and C. Rozell, “Estimation and dynamic updating of time-varying signals with sparse variations, ” in 2011 IEEE International Conference on Acoustics, Speech and Signal Pr ocessing (ICASSP) , May 2011, pp. 3908–3911. [8] J. Ziniel, L. C. Potter , and P . Schniter , “Tracking and smoothing of time-varying sparse signals via approximate belief propagation, ” in 2010 Confer ence Record of the F orty F ourth Asilomar Conference on Signals, Systems and Computers , Nov . 2010, pp. 808–812. [9] N. V aswani and W . Lu, “Modified-CS: Modifying compressive sensing for problems with partially known support, ” IEEE T ransactions on Signal Pr ocessing , vol. 58, no. 9, pp. 4595–4607, Sep. 2010. [10] A. Charles, M. S. Asif, J. Romberg, and C. Rozell, “Sparsity penalties in dynamical system estimation, ” in 2011 45th Annual Conference on Information Sciences and Systems , Mar . 2011, pp. 1–6. [11] E. C. Hall and R. M. W illett, “Dynamical models and tracking regret in online con ve x programming, ” in Proceedings of the 30th International Confer ence on Machine Learning , S. Dasgupta and D. McAllester , Eds. PMLR, Feb. 2013, pp. 579–587. [12] A. S. Charles, A. Bala voine, and C. J. Rozell, “Dynamic filtering of time-varying sparse signals via ` 1 minimization, ” IEEE T ransactions on Signal Pr ocessing , vol. 64, no. 21, pp. 5644–5656, Nov . 2016. [13] A. Balav oine, C. J. Rozell, and J. Romberg, “Discrete and continuous- time soft-thresholding for dynamic signal recovery , ” IEEE Tr ansactions on Signal Processing , vol. 63, no. 12, pp. 3165–3176, 2015. [14] A. S. Charles, N. P . Bertrand, J. Lee, and C. J. Rozell, “Earth-mover’ s distance as a tracking regularizer , ” in 2017 IEEE 7th International W ork- shop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP) , Dec 2017, pp. 1–5. [15] N. P . Bertrand, J. Lee, A. S. Charles, P . Dunn, and C. J. Rozell, “Sparse dynamic filtering via earth mover’ s distance regularization, ” in Pr oceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP), Calgary , Alberta, Canada , Apr 2018. [16] E. Candès and J. Romberg, “ ` 1 -magic: recovery of sparse signals via con ve x programming, ” California Institute of T echnology , T ech. Rep., 2005. [17] A. C. Sankaranarayanan, P . K. Turag a, R. G. Baraniuk, and R. Chellappa, “Compressiv e acquisition of dynamic scenes, ” in Computer V ision – ECCV 2010 , ser . Lecture Notes in Computer Science. Springer , Berlin, Heidelberg, Sep. 2010, pp. 129–142. [18] M. S. Asif, L. Hamilton, M. Brummer, and J. Romberg, “Motion-adaptive spatio-temporal regularization for accelerated dynamic MRI, ” Magnetic Resonance in Medicine , vol. 70, no. 3, pp. 800–812, Sep. 2013. [19] A. Carmi, P . Gurfil, and D. Kanevsky , “Methods for sparse signal recovery using Kalman filtering with embedded pseudo-measurement norms and quasi-norms, ” IEEE T ransactions on Signal Processing , vol. 58, no. 4, pp. 2405–2409, Apr . 2010. [20] G. Monge, Mémoire sur la théorie des déblais et des remblais . Paris: De l’Imprimerie Royale, 1781, oCLC: 51928110. [21] Y . Rubner, C. T omasi, and L. J. Guibas, “The earth mover’ s distance as a metric for image retrieval, ” International Journal of Computer V ision , vol. 40, no. 2, pp. 99–121, Nov . 2000. [22] H. Ling and K. Okada, “ An efficient earth mover’ s distance algorithm for robust histogram comparison, ” IEEE Tr ansactions on P attern Analysis and Machine Intelligence , vol. 29, no. 5, pp. 840–853, May 2007. [23] R. Gupta, P . Indyk, and E. Price, “Sparse recov ery for earth mover distance, ” in 2010 48th Annual Allerton Conference on Communication, Contr ol, and Computing (Allerton) , Sep. 2010, pp. 1742–1744. [24] L. Schmidt, C. Hegde, and P . Indyk, “The constrained earth mov er distance model, with applications to compressive sensing, ” in 10th Inter- national Conference on Sampling Theory and Applications (SAMPT A) , 2013. [25] D. Mo and M. F . Duarte, “Compressive parameter estimation with earth mover’ s distance via K-median clustering, ” in Pr oc. SPIE 8858, W avelets and Sparsity XV , vol. 8858, 2013, pp. 88 581P–8858–10. [26] M. Beckmann, “ A continuous model of transportation, ” Econometrica: Journal of the Econometric Society , pp. 643–660, 1952. [27] W . Li, E. K. Ryu, S. Osher , W . Y in, and W . Gangbo, “ A parallel method for earth mo ver’ s distance, ” Journal of Scientific Computing , vol. 75, no. 1, pp. 182–197, 2018. [28] E. K. Ryu, W . Li, P . Y in, and S. Osher , “Unbalanced and partial L 1 Monge-Kantorovich problem: a scalable parallel first-order method, ” Journal of Scientific Computing , pp. 1–18, Nov . 2017. [29] J. Karlsson and A. Ringh, “Generalized Sinkhorn iterations for regular- izing inv erse problems using optimal mass transport, ” SIAM Journal on Imaging Sciences , vol. 10, no. 4, pp. 1935–1962, Jan. 2017. [30] H. Janati, T . Bazeille, B. Thirion, M. Cuturi, and A. Gramfort, “Group lev el MEG/EEG source imaging via optimal transport: Minimum W asserstein estimates, ” in Information Processing in Medical Imaging , ser . Lecture Notes in Computer Science, A. C. S. Chung, J. C. Gee, P . A. Y ushke vich, and S. Bao, Eds. Cham: Springer International Publishing, 2019, pp. 743–754. [31] H. Janati, M. Cuturi, and A. Gramfort, “W asserstein regularization for sparse multi-task regression, ” arXiv:1805.07833 [cs, stat] , May 2018. [32] R. Sinkhorn, “ A relationship between arbitrary positi ve matrices and doubly stochastic matrices, ” The Annals of Mathematical Statistics , vol. 35, no. 2, pp. 876–879, 1964. [33] M. Cuturi, “Sinkhorn Distances: Lightspeed Computation of Optimal T ransport, ” in Advances in Neural Information Processing Systems 26 , C. J. C. Bur ges, L. Bottou, M. W elling, Z. Ghahramani, and K. Q. W einberger , Eds. Curran Associates, Inc., 2013, pp. 2292–2300. [34] H. Janati, T . Bazeille, B. Thirion, M. Cuturi, and A. Gramfort, “Group lev el meg/ee g source imaging via optimal transport: minimum wasserstein estimates, ” in International Conference on Information Processing in Medical Imaging . Springer , 2019, pp. 743–754. [35] H. Janati, M. Cuturi, and A. Gramfort, “W asserstein regularization for sparse multi-task regression, ” arXiv preprint , 2018. [36] M. Grant and S. Boyd, “CVX: MA TLAB software for disciplined conve x programming, version 2.1, ” Mar . 2014. [37] T . K olda, R. Lewis, and V . T orczon, “Optimization by direct search: New perspectiv es on some classical and modern methods, ” SIAM Review , vol. 45, no. 3, pp. 385–482, Jan. 2003. [38] L. Cornelissen, S.-E. Kim, P . L. Purdon, E. N. Brown, and C. B. Berde, “ Age-dependent electroencephalogram (EEG) patterns during sev oflurane general anesthesia in infants, ” eLife , vol. 4, p. e06513, Jun. 2015. [39] J. Aru, J. Aru, V . Priesemann, M. Wibral, L. Lana, G. Pipa, W . Singer , and R. V icente, “Untangling cross-frequency coupling in neuroscience, ” Curr ent Opinion in Neurobiolo gy , vol. 31, pp. 51–61, Apr. 2015. [40] T . J. Gardner and M. O. Magnasco, “Sparse time-frequency represen- tations, ” Pr oceedings of the National Academy of Sciences , vol. 103, no. 16, pp. 6094–6099, Apr . 2006. [41] D. J. Thomson, “Multitaper analysis of nonstationary and nonlinear time series data, ” in Nonlinear and Nonstationary Signal Pr ocessing , W . J. Fitzgerald, Ed. Cambridge ; New Y ork: Cambridge University Press, 2000, pp. 317–394, oCLC: ocm45339400. [42] F . Auger and P . Flandrin, “Improving the readability of time-frequency and time-scale representations by the reassignment method, ” IEEE T ransactions on Signal Pr ocessing , vol. 43, no. 5, pp. 1068–1089, May 1995. [43] S. A. Fulop and K. Fitz, “ Algorithms for computing the time-corrected instantaneous frequency (reassigned) spectrogram, with applications, ” The Journal of the Acoustical Society of America , vol. 119, no. 1, pp. 360–371, Jan. 2006. [44] C. Kemere, M. F . Carr , M. P . Karlsson, and L. M. Frank, “Rapid and continuous modulation of hippocampal network state during exploration of new places, ” PLOS ONE , vol. 8, no. 9, p. e73114, Sep. 2013. [45] H. Zhang, A. J. W atrous, A. Patel, and J. Jacobs, “Theta and alpha oscillations are trav eling waves in the human neocortex, ” Neur on , vol. 98, no. 6, pp. 1269–1281.e4, Jun. 2018. [46] G. B. Ermentrout and D. Kleinfeld, “Tra veling Electrical W aves in Cortex: Insights from Phase Dynamics and Speculation on a Computational Role, ” Neur on , vol. 29, no. 1, pp. 33–44, Jan. 2001. [47] Y . Kuramoto, “Rhythms and turbulence in populations of chemical oscillators, ” Physica A: Statistical Mechanics and its Applications , vol. 106, no. 1, pp. 128–143, Mar . 1981.