# Novel Structured Low-rank algorithm to recover spatially smooth exponential image time series

We propose a structured low rank matrix completion algorithm to recover a time series of images consisting of linear combination of exponential parameters at every pixel, from under-sampled Fourier measurements. The spatial smoothness of these parameters is exploited along with the exponential structure of the time series at every pixel, to derive an annihilation relation in the k-t domain. This annihilation relation translates into a structured low rank matrix formed from the k-t samples. We demonstrate the algorithm in the parameter mapping setting and show significant improvement over state of the art methods.

## Authors

• 2 publications
• 24 publications
• ### Calibration-free B0 correction of EPI data using structured low rank matrix recovery

We introduce a structured low rank algorithm for the calibration-free co...
04/20/2018 ∙ by Arvind Balachandrasekaran, et al. ∙ 0

• ### Tight Risk Bound for High Dimensional Time Series Completion

Initially designed for independent datas, low-rank matrix completion was...
02/15/2021 ∙ by Pierre Alquier, et al. ∙ 0

• ### Recovery of Piecewise Smooth Images from Few Fourier Samples

We introduce a Prony-like method to recover a continuous domain 2-D piec...
02/03/2015 ∙ by Greg Ongie, et al. ∙ 0

• ### Bayesian Learning for Low-Rank matrix reconstruction

We develop latent variable models for Bayesian learning based low-rank m...
01/23/2015 ∙ by Martin Sundin, et al. ∙ 0

• ### Image Restoration: Structured Low Rank Matrix Framework for Piecewise Smooth Functions and Beyond

Recently, mapping a signal/image into a low rank Hankel/Toeplitz matrix ...
12/12/2020 ∙ by Jian-Feng Cai, et al. ∙ 0

• ### Low-rank Matrix Completion in a General Non-orthogonal Basis

This paper considers theoretical analysis of recovering a low rank matri...
12/14/2018 ∙ by Abiy Tasissa, et al. ∙ 0

• ### Time-Series Analysis via Low-Rank Matrix Factorization: Applied to Infant-Sleep Data

We propose a nonparametric model for time series with missing data based...
04/09/2019 ∙ by Sheng Liu, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 introduction

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 [1]

. 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 [2]

, 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 in

space, 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, [3] 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 [3]. This decoupled approach is less constrained and cannot handle efficient 2-D sampling patterns; Cartesian patterns with fully sampled lines were considered in [3]. In [4], 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 [5][6],[7] 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 :

 ρ[r,n]=L∑i=1αi(r) βi(r)n. (1)

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 .

It is well known that the exponential signal in (1), corresponding to a specific spatial location , can be annihilated by a filter [1]:

 ∑k∈θρ[r,n−k] h[r,k]=0=μ[r,n],  ∀r. (2)

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

 ^ρ[k,n]⊗c[k,n]=0. (3)

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

 T(^ρ) c=Qc=0 (4)

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

 ^ρ[k,n]⊗d[k,n]=0, (5)

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

 b=A(^ρ)+η, (6)

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:

 ^ρ⋆=argmin^ρ∥T(^ρ)∥p+λ2∥A(^ρ)−b∥22 (7)

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):

 Hn =[T(^ρn−1)T(^ρn−1)∗R+ϵnI]p2−1 (8) ^ρn =argmin^ρ12∥H12nT(^ρ)∥2F+λ2∥A(^ρ)−b∥22 (9)

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

 ^ρ∗=argmin^ρ12M∑i=1∥h(i)T(^ρ)∥22+λ2∥A(^ρ)−b∥22 (10)

Note that the term is the linear convolution between the 3-D sequences and . In the GIRAF algorithm [2], 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,

 h(i)T(^ρ)=(h(i)Nt…h(i)1)⎛⎜ ⎜⎝T(^ρk)…T(^ρ1)⋮⋮⋮T(^ρk+Nt−1)⋯T(^ρNt)⎞⎟ ⎟⎠ (11)

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.

### 3.2 Weight-Update

We consider the Gram matrix in (8) that has the following structure:

 ⎛⎜ ⎜⎝R1,1R1,2…R1,Nt⋮⋮⋯⋮RNt,1RNt,2⋯RNt,Nt⎞⎟ ⎟⎠ (12)

with column and row partitions and is a matrix block of dimension . We obtain the general expression for the matrix block as.

 Rp,q=k∑i=1T(^ρp+i−1)T(^ρq+i−1)∗ (13)

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

 Pi,j=P∗ΛCirc(^ρi)Circ(^ρj)∗Circ(g)PΛ (14)

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

 H=[U(Λ+ϵI)U∗]p2−1=U(Λ+ϵI)p2−1U∗,

where is the eigen decomposition of . Hence, one choice of the matrix square root is .

## 4 Experiments

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).

## 5 Conclusion

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.