The recovery of linear combination of damped exponentials from few of their measurements is a classical problem in signal processing, starting with the seminal work of Prony 
. Early work, which focused on uniform sampling, relied on the existence of a filter that annihilates the uniform measurements. Recently, several researchers have extended this idea to the non-uniform sampling setting, which is a more efficient measurement strategy. These methods pose the estimation as the completion of a structured matrix (e.g Hankel or Toeplitz), whose entries are the signal samples, from few of its measurements. Specifically, the annihilation filters are the null space vectors of the Toeplitz matrix, or equivalently the matrix is low-rank.
The above exponential estimation problem is of high significance in several medical imaging applications, including MR spectroscopic imaging and MR parameter mapping. These methods perform a pixel-by-pixel estimation of the exponential parameters (e.g. frequencies, decay rate, amplitudes) from an image time series. The parameters provide valuable clues about the abnormal metabolic activity and tissue micro-structural changes, which are early markers for neurological disorders. The parameters often vary smoothly in space since they depend on the underlying tissue microstructure. Unfortunately, the acquisition of the image time series at high spatial resolution results in prohibitively long scan time. A common approach to accelerate the acquisition is to acquire under sampled data, followed by reconstruction using low rank and sparsity priors.
We introduce a novel formulation that directly exploits the exponential structure of the time series and the spatial smoothness of the exponential parameters. Specifically, we couple the time series estimation problems at all the pixels into a single annihilation relation, involving a spatially smooth annihilation filter. This approach exploits the spatial smoothness naturally; we do not require additional spatial regularization priors. The annihilation relations translate to a low-rank penalty on a Toeplitz matrix, whose entries are the space samples of the time series.
A challenge with the above formulation is the large size of the Toeplitz matrix. The direct implementation of the structured low-rank problem will be prohibitive for high resolution applications. We generalize our generic iteratively reweighted annihilating filter (GIRAF) framework, originally introduced for piecewise smooth images 
, to accelerate the computations. In particular, we use an iterative reweighted formulation, where we exploit the Toeplitz structure of the matrix to avoid its direct computation and storage. The GIRAF scheme approximates the linear convolutions resulting from the Toeplitz multiplications with circular convolutions, which allowed their evaluation using Fast Fourier Transform (FFT). These approximations were enabled by the fast decay of the image Fourier coefficients towards the edges of the observed region. Unfortunately, this approach breaks down in our setting since we rely on annihilation relations inspace, where the signal has considerable intensity at the first image. Hence, we modify the GIRAF steps using a hybrid circular-linear convolution strategy. These generalizations provide a fast algorithm that can be readily applied to large-scale exponential estimation problems. Our validations using MRI data clearly demonstrate the benefits offered by the proposed algorithm, which can simultaneously exploit the exponential structure and the spatial smoothness of the parameters.
This work has similarities to recent works on structured low-rank priors for the estimation of piecewise smooth signals and exponential signals. For example,  exploited the smoothness of the signal using a low-rank prior on a wavelet transform weighted Hankel matrix. This work does not exploit the exponential structure of the signal. In addition, the recovery of each slice is decoupled into a separate problem to keep the computational complexity reasonable . This decoupled approach is less constrained and cannot handle efficient 2-D sampling patterns; Cartesian patterns with fully sampled lines were considered in . In , the linear predictability of the time series at each pixel is exploited individually. Since this decoupled strategy is less constrained, the authors rely on additional low-rank and wavelet sparsity regularization penalties; the optimization of several regularization parameters is also challenging. While the piecewise smoothness and sparsity has also been exploited by several researchers using the structured low-rank setting , they do not exploit the exponential signal structure.
2 Problem formulation
We focus on the recovery of a 3-D volume , consisting of a linear combination of damped exponentials at every pixel, from noisy and undersampled measurements denoted by . We will now introduce the annihilation relations, and the resulting structured low-rank priors.
2.1 Annihilation of spatially smooth exponentials
We model the temporal signal at the spatial location as :
Note that the temporal signal is a linear combination of exponentials, are the amplitudes and are the exponential parameters that are dependent on the physiological process. For example, in mapping applications, the exponential parameters , where is the time between two image frames and is the relaxation parameter of the tissue component (e.g. gray matter, CSF, white matter) at the voxel indexed by .
where (2) represents a 1-D convolution and is a valid index set. The coefficients of at are specified by . Computing the 2-D Fourier transform (along the spatial dimension, denoted by ) of (2), we obtain
Here, are the spatial (2-D) Fourier coefficients of the measurements and denotes 3D convolution. Similarly, denotes the spatial (2-D) Fourier coefficients of . We assume that the parameters describing the physiological process are bandlimited functions of the spatial variable ; this implies that is a 3D FIR filter. In particular, the extent of the filter in the spatial frequency () dimension controls the spatial smoothness of the filter , while its extent in the temporal dimension () depends on the number of exponentials in the model (1). We express (3) in a matrix form as
where is a linear operator that maps a 3D volume to a lifted matrix . Here the lifted matrix has a multi-fold Toeplitz structure, such that corresponds to the 3-D convolution between and the FIR filter . The number of columns of the matrix is equal to the assumed support of the filter. Likewise, if the measurements are restricted to a cube shaped volume in , the rows correspond to the valid convolutions.
In practice, the support of the filter is unknown. Let us assume that (3) is satisfied by a minimal filter of support . In this case, we observe that
where and is any FIR filter. Note that the support of , denoted by is larger than ; i.e, . Hence, if we over-estimate the support of the filter , there will be many linearly independent annihilating filters in the nullspace of , or equivalently is low rank.
To avoid oversmoothing of the parameter maps, we choose the spatial support of the 3D filter to be large. This implies that the lifted matrix is rectangular (more columns than rows) in the parameter mapping setting.
2.2 Recovery using structured low-rank matrix priors
In this paper, we focus on the recovery of weighted MR images from their undersampled multichannel encoded measurements, denoted by
where is a linear operator representing Fourier under sampling and multiplication of coil sensitivities with . We pose the recovery of the signal (1) as the structured low rank matrix recovery problem:
where is a regularization parameter and is the Schatten norm, defined as ;
are the singular values of. Here, denotes the structured multifold Toeplitz matrix, whose entries are the samples of ; corresponds to the 2-D Fourier coefficients of .
3 Optimization Algorithm
We rely on the iterative re-weighted least squares (IRLS) based algorithm to solve (7). The basic idea is to use the identity , where . We use an alternating minimization algorithm, which alternates between the following two sub-problems, to solve (7):
where is added to stabilize the inverse. We will now focus on an efficient implementation of the subproblems.
3.1 Least squares-Update
Let the rows of be denoted by . Substituting for in (9), we obtain
Note that the term is the linear convolution between the 3-D sequences and . In the GIRAF algorithm , we relied on the approximation of the linear convolution by circular convolutions to accelerate its computations using FFT. We exploited the fast decay of the Fourier coefficients in GIRAF, which made the approximations valid. Here we rely on the annihilation relations in the domain, where the signal does not decay rapidly in the temporal dimension. Hence we use a hybrid approach that consists of performing linear convolutions along time and circular convolutions along the spatial dimension.
Let and consist of and frames respectively and let denote . Let each frame of the filter be of dimension . Now can be simplified as,
In the above equation, represents a Toeplitz matrix formed from ; the product corresponds to the 2-D convolution between the frame of the filter and frame of . We safely approximate each of the 2-D convolutions by circular convolutions as in GIRAF and compute them efficiently using Fast Fourier transforms.
We consider the Gram matrix in (8) that has the following structure:
with column and row partitions and is a matrix block of dimension . We obtain the general expression for the matrix block as.
We observe that the Toeplitz matrix can be expressed as , where is a matrix that restricts the convolution onto a valid index set represented by , is a circulant matrix formed from and
is a matrix representing zero-padding outside the filter support; here are the spatial dimensions of . Now can be simplified as
where the entries of are generated from the array given by , conj denotes the conjugate operation and denotes point-wise multiplication.
We evaluate the weight matrix as
where is the eigen decomposition of . Hence, one choice of the matrix square root is .
A fully sampled axial 2-D dataset was acquired on a Siemens 3T Trio scanner using 12 coils and a turbo spin echo sequence. The following scan parameters were used in the acquisition: Matrix size:128128, FOV: 2222 , TR = 2500 ms and slice thickness = 5 mm. The weighted images were obtained for 12 equispaced echo (TE) times ranging from 10 to 120 ms. Post reconstruction, the maps were obtained by fitting a mono-exponential model to each voxel.
We first study the effect of the filter size (dimensions of the block Toeplitz matrix) on the SNR of the images recovered from 8-fold undersampled multichannel Fourier data in Table. 1. We observe that the spatial dimensions of the filter have the largest influence on the SNR, as seen from Table. 1.(b). Specifically, we observe that filters with smaller support provides improved results, which demonstrates the benefit of exploiting spatial smoothness; a filter with size fails to exploit smoothness. We also observe from the Table.1.(a) the benefit of using temporal annihilation relations. Specifically, we obtain a 0.5 dB improvement over the filter with size 122x122x1, which only exploits joint sparsity, by increasing the length of the filter along time.
In Fig. 2, we compare the proposed approach with low rank algorithm on the recovery of single coil (coil compressed) weighted data from uniform random measurements. We chose a filter of size and a Schatten . We observe that the proposed scheme provides lower errors (see caption for details).
In Fig. 3, we compare the two methods on the recovery of multi-channel weighted data from 12-fold (3-fold variable density + 4-fold 22 Cartesian) undersampled data. We chose a filter of size and a Schatten for the proposed scheme (see caption for details).
We introduced a novel structured low rank algorithm to recover a 3-D volume consisting of a linear combination of exponentials, from undersampled Fourier measurements. A convolution relation was obtained between the Fourier samples and a 3-D FIR filter by exploiting the exponential structure of the time series at every pixel and the smoothness of the parameters. To speed up the computations, a hybrid approach was employed which consisted of performing circular and linear convolutions along the spatial and temporal dimensions respectively. The algorithm was applied in the context of MR parameter mapping and the reconstructed images and maps were sharper and had fewer errors than those obtained from state of the art methods.
-  P. Stoica and R. L. Moses, Introduction to spectral analysis. Prentice hall Upper Saddle River, 1997, vol. 1.
-  G. Ongie and M. Jacob, “GIRAF: A fast algorithm for structured low-rank matrix recovery.” [Online]. Available: http://arxiv.org/abs/1609.07429
-  D. Lee, K. H. Jin, E. Y. Kim, S.-H. Park, and J. C. Ye, “Acceleration of mr parameter mapping using annihilating filter-based low rank hankel matrix,” MRM, 2016.
-  X. Peng and et. al, “Accelerated exponential parameterization of T2 relaxation with model-driven low rank and sparsity priors (MORASA),” MRM, 2016.
-  J. Haldar, “Low-rank modeling of local k -space neighborhoods (loraks) for constrained mri,” IEEE Trans. Med. Imaging, no. 3, pp. 668–681, March 2014.
-  K. H. Jin, D. Lee, and J. C. Ye, “A general framework for compressed sensing and parallel MRI using annihilating filter based low-rank hankel matrix,” 2015. [Online]. Available: http://arxiv.org/abs/1504.00532
-  G. Ongie and M. Jacob, “Off-the-Grid Recovery of Piecewise Constant Images from Few Fourier Samples,” SIAM J Imaging Sciences, in press, 2015.