Cover Tree Compressed Sensing for Fast MR Fingerprint Recovery

TL;DR

This paper introduces a cover tree-based compressed sensing method for fast MRI fingerprint recovery, significantly reducing computation time while maintaining accuracy by leveraging approximate nearest neighbor searches on smooth manifolds.

Contribution

It presents a novel integration of cover trees with the IPG algorithm for MRI compressed sensing, achieving logarithmic projection costs and substantial computational speedups.

Findings

2-3 orders of magnitude reduction in computation time

Maintains similar or better reconstruction accuracy

Effective on Magnetic Resonance Fingerprinting data

Abstract

We adopt data structure in the form of cover trees and iteratively apply approximate nearest neighbour (ANN) searches for fast compressed sensing reconstruction of signals living on discrete smooth manifolds. Levering on the recent stability results for the inexact Iterative Projected Gradient (IPG) algorithm and by using the cover tree's ANN searches, we decrease the projection cost of the IPG algorithm to be logarithmically growing with data population for low dimensional smooth manifolds. We apply our results to quantitative MRI compressed sensing and in particular within the Magnetic Resonance Fingerprinting (MRF) framework. For a similar (or sometimes better) reconstruction accuracy, we report 2-3 orders of magnitude reduction in computations compared to the standard iterative method which uses brute-force searches.

Figures

Equations

$$y\approx A{x}_0,$$

$$\hat{x}={\mathrm{argmin}}_{x\in \mathcal{C}}||y-Ax|{|}_2^2.$$

$$x^{k+1}={\mathcal{P}}_{\mathcal{C}}\left(x^k-\mu A^H(A{x}^k-y)\right),$$

$$\alpha ||x-x^{\prime }|{|}^2⩽||A(x-x^{\prime })|{|}^2⩽\beta ||x-x^{\prime }|{|}^2$$

$$x^k={\mathcal{P}}_{\mathcal{C}}^{\epsilon }\left(x^{k-1}-\mu A^H(A{x}^k-y)\right).$$

$$\mathcal{P}_{\mathcal{C}}^{\varepsilon}(x)\in\Big{\{}u\in\mathcal{C}:\,|\!|u-x|\!|\leqslant(1+\varepsilon)|\!|\mathcal{P}_{\mathcal{C}}(x)-x|\!|\Big{\}}.$$

$$\sqrt{\epsilon +{\epsilon }^2}⩽\delta \sqrt{\alpha }/\left|\left|\left|A\right|\right|\right|\text{and}\beta <(2-2\delta +{\delta }^2)\alpha$$

$$||x^k-x_0||⩽{\rho }^k||x_0||+\frac{{\kappa }_w}{1-\rho }w$$

$$\rho =\sqrt{\frac{1}{\mu \alpha }-1}+\delta ,{\kappa }_w=2\frac{\sqrt{\beta }}{\alpha }+\sqrt{\mu }\delta ,w=||y-A{x}_0||.$$

$$\tilde{\mathcal{C}}:=\bigcup _{i=1}^d\{{\psi }_i\}\in {\mathbb{C}}^{\tilde{n}}$$

$$X_j\in \text{cone}(\tilde{\mathcal{C}}),\forall j=1,\dots ,J,$$

$$\text{cone}(\tilde{\mathcal{C}}):=\{x\in {\mathbb{C}}^{\tilde{n}}:x/\gamma \in \tilde{\mathcal{C}}\text{for some}\gamma >0\}.$$

$$\mathcal{C}:=\prod _{j=1}^J\text{cone}(\tilde{\mathcal{C}})$$

$$i^{\ast }={\mathrm{argmin}}_i||Z_j-{\psi }_i/||{\psi }_i||||$$

$$X_j^{k+1}={\mathcal{P}}_{\text{cone}(\tilde{\mathcal{C}})}(Z_j)={\gamma }_j{\psi }_{i^{\ast }}$$

$$2^{O({\dim }_D(\tilde{\mathcal{C}}))}\log \Delta +(1/\epsilon {)}^{O({\dim }_D(\tilde{\mathcal{C}}))}$$

$${\psi }_{i^{\ast }}=\text{(1+ε)-ANN}\left(\mathcal{T},\frac{Z_j}{||Z_j||},\frac{X_j^k}{{\gamma }_j}\right),\forall j=1,\dots ,J.$$

$$Y^l\approx {\mathcal{F}}_{{\Omega }_l}(X^l),\forall l=1,\dots \tilde{n},$$

$$(X^l{)}^{k+1}={\mathcal{P}}_{\text{cone}(\tilde{\mathcal{C}})}^{\epsilon }\left((X^l{)}^k-\mu {\mathcal{F}}_{{\Omega }_l}^H\left({\mathcal{F}}_{{\Omega }_i}((X^l{)}^k)-Y^l\right)\right)$$

Full text

Index Terms—  Data driven compressed sensing, iterative projected gradient, magnetic resonance fingerprinting, cover tree, approximate nearest neighbour search.

1 Introduction

Compressed sensing (CS) algorithms adopt efficient signal models to achieve accurate reconstruction given small number of non-adaptive measurements [1, 2]. CS consists of a linear sampling model:

[TABLE]

where $y$ is a noisy $m$-dimensional vector of measurements taken from the ground truth signal $x_{0}$ in a high dimension $n\gg m$ through using a linear mapping $A$. The success of CS for stably inverting such ill-posed systems relies on the fact that natural signals often follow low-dimensional models $\mathcal{C}$ with limited degrees of freedom compared to their original ambient dimension, and that certain random sampling schemes can compactly (in terms of the measurements) preserve the key signal information (see e.g. [3] for an overview on different CS models and their corresponding sample complexities). Variational approach for CS reconstruction consists of solving e.g. the following constrained least-square problem:

[TABLE]

where the constraint takes into account the signal model.

First order methods based on iterative local gradient and projection (or also proximal) updates are popular for CS recovery due to their scalability for big data problems [4] and their flexibility in handling various complicated (convex/nonconvex and sometimes combinatorial) models [5]. Starting from an initial point e.g. $x^{k=0}=\mathbf{0}$ and for a given step size $\mu$, the solution of Iterative Projected Gradient (IPG) iteratively updates as follows:

[TABLE]

where the projection is defined as $\mathcal{P}_{\mathcal{C}}(x)\in\operatorname{argmin}_{u\in\mathcal{C}}|\!|x-u|\!|$ and throughout, $|\!|.|\!|$ stands for the Euclidean norm.

In this work we consider compressed sensing of signals living on a low-dimensional smooth manifold. In this case, the model projection might not be computationally easy for certain complex manifolds. A popular approach consists of collecting discrete samples of the manifold in a dictionary [6, 7], i.e. a data driven CS approach, and approximate the projection step by searching the nearest atom in this dictionary. Depending on the manifold complexity and the desired approximation accuracy, the size of this dictionary can be very large in practice. In this case performing an exhaustive nearest neighbour (NN) search could become computationally challenging.

We address this shortcoming by preprocessing the dictionary and build a tree structure suitable for fast approximate nearest neighbour (ANN) searches. Recently, it as been shown that provided enough measurements and for moderate approximations the inexact IPG achieves similar solution accuracy as for the exact algorithm [8]. In this regard, we use the cover tree structure [9, 10] which is proposed for fast ANN searches on large datasets with low intrinsic dimensionality and with a search complexity $O(\log(d))$ in time, where $d$ is the dataset population (in contrast with the linear complexity of an exact exhaustive search $O(d)$).

We apply our results to accelerate the Magnetic Resonance Fingerprint (MRF) recovery from under-sampled k-space measurements; an emerging field of study in computational medical imaging [11]. Since the MR fingerprints belong to the low-dimensional manifold of Bloch dynamic equations, we show that using the cover tree ANN search significantly shortcuts the computations compared to an exact IPG with exhaustive searches as in [12].

2 Preliminaries

The following embedding assumption plays a critical role in stability analysis of the (convex/nonconvex) IPG algorithm [13, 5, 8]:

Definition 1**.**

$A$* is bi-Lipschitz w.r.t. $\mathcal{C}$, if $\forall x,x^{\prime}\in\mathcal{C}$ there exists constants $0<\alpha\leqslant\beta$ such that*

[TABLE]

Dealing with complicated signal models, the projection step can bring a huge computational burden at each iteration, and thus a natural line of thought is to perform this step with cheaper (available) approximations. The inexact IPG takes the following form:

[TABLE]

We consider the class of $(1+\varepsilon)$-approximate projections: for any $\varepsilon\geqslant 0$:

[TABLE]

The following result states that if the embedding holds (with a better conditioning compared to the exact IPG case $\varepsilon=0$) then the inexact algorithm still achieves a solution accuracy comparable to the exact IPG [8]:

Theorem 1**.**

Assume $(A,\mathcal{C})$ is bi-Lipschitz, ${x_{0}}\in\mathcal{C}$ and that

[TABLE]

for $\varepsilon\geqslant 0$ and some constant $\delta\in[0,1)$. Set the step size $\left((2-2\delta+\delta^{2}){\alpha}\right)^{-1}<\mu\leqslant{\beta}^{-1}$. The sequence generated by Algorithm (5) obeys the following bound:

[TABLE]

where

[TABLE]

Remark 1*.*

Increasing $\varepsilon$ slows down the rate $\rho$ of linear convergence, however the iterations are potentially cheaper. Too large approximations $\varepsilon$ may result in divergence.

Remark 2*.*

After a finite $K=O(\log(\tau^{-1}))$ number of iterations Algorithm (5) achieves the solution accuracy $|\!|x^{K}-x_{0}|\!|=O(w)+\tau$, for any $\tau>0$.

3 Data driven CS in product space

Data driven CS corresponds to cases where in the absence of an algebraic (semi-algebraic) physical model one resorts to collecting a large number of data samples in a dictionary and use it e.g. as a point cloud model for CS reconstruction [3]. Data driven CS finds numerous applications e.g. in hyperspectral imaging (HSI) [14], mass spectroscopy (e.g. MALDI imaging) [15] and MRF recovery [11, 12] just to name a few. For instance the USGS Hyperspectral library111http://speclab.cr.usgs.gov contains the spectral signatures of thousands of substances. This side-information is shown to enhance CS recovery and classification [14, 16]. Here the model is discrete and nonconvex:

[TABLE]

and it corresponds to atoms of a (possibly complex-valued) dictionary $\Psi$ customized for e.g. HSI, MALDI or MRF data.

A multichannel image such as HSI, MALDI or MRF can be represented by a ${\tilde{n}}\times J$ matrix $X$, where $n={\tilde{n}}J$ is the total number of pixels, $J$ is the spatial resolution and ${\tilde{n}}$ is the number of channels. Under the pixel purity assumption each spatial pixel corresponds to a certain material with specific signature. Considering a notion of signal intensity in this model (e.g. due to variable illumination conditions in HSI, or a non-uniform proton density in MRF), the columns of $X$ belong to the cone associated with the signatures (7). Denoting by $X_{j}$ the $j$th column of $X$ (i.e. a multi-dimensional spatial pixel), we have

[TABLE]

where the cone is also a discrete set defined as

[TABLE]

Here, $\gamma$ corresponds to the signal intensity (per spatial pixel).

The CS sampling and reconstruction follows (1) and (2) by setting $x^{*}:=X_{\text{vec}}$ (by $X_{\text{vec}}$ we denote the vector-rearranged form of the matrix $X$) and the only update is that now the solution lives in a product space of the same model i.e.

[TABLE]

In a single dimensional setting i.e. $\mathcal{C}=\text{cone}(\widetilde{\mathcal{C}})$, there will be no need for an IPG algorithm to solve problem (2). If the embedding holds with a ${\alpha}>0$ (i.e. an identifiability condition) then the atom with a maximum coherence $\langle A\psi_{i}/|\!|A\psi_{i}|\!|,y\rangle$ with the measurements identifies the correct signature in a single iteration. However such a direct approach in a multi-dimensional setting $J\gg 1$ (e.g. when the sampling model non-trivially combines columns of $X$) generally has an exponential $O(d^{J})$ complexity because of the combinatorial nature of the product space constraints. A tractable scheme which has been frequently considered for this case e.g. in [12] would be an IPG type algorithm (3).

After the gradient update $Z_{\text{vec}}$ = $X_{\text{vec}}^{k}+\mu A^{H}(y-AX_{\text{vec}}^{k})$, the projection onto the product model (9) decouples into separate cone projections for each spatial pixel $j$. The cone projection $\mathcal{P}_{\text{cone}(\widetilde{\mathcal{C}})}(Z_{j})$ follows the steps below:

[TABLE]

where, $\gamma_{j}=\max\left(\text{real}(\langle Z_{j},\psi_{i^{*}}\rangle)/{|\!|\psi_{i^{*}}|\!|^{2}},0\right)$ is the per-pixel signal intensity.

As a result, the IPG algorithm breaks down the computations into local updates namely, the gradient and projection steps for which an exact projection step e.g. by using a brute-force nearest neighbour (NN) search, has complexity $O(Jd)$ in time.

4 Cover tree for fast nearest neighbour search

With the data driven CS formalism and discretization of the model the projection step of IPG reduces to a search for the nearest atom in each of the product spaces, however in a potentially very large $d$ size dictionary. Thus search strategies with linear complexity in $d$, e.g. an exhaustive search, can be a serious bottleneck for solving such problems. A very well-established approach to overcome the complexity of an exhaustive NN search on a large dataset consists of hierarchically partitioning the solution space and forming a tree structure whose nodes representing those partitions, and then using a branch-and-bound method on the resulting tree for a fast Approximate NN (ANN) search with $o(d)$ complexity. In this regard, we address the computational shortcoming of the projection step in the exact IPG by preprocessing $\widetilde{\mathcal{C}}$ and form a cover tree structure suitable for fast ANN searches [9, 10].

A cover tree is a levelled tree whose nodes (i.e. associated data points) at different scales form covering nets for data points at multiple resolutions. Denote by $\mathcal{S}_{l}$ the set of nodes appearing at scale $l=1,\ldots,L_{\max}$ and by $\sigma:=\max_{\psi\in\widetilde{\mathcal{C}}}|\!|\psi_{\text{root}}-\psi|\!|$ the maximal coverage by the root, then a cover tree structure must have the following three properties:

Nesting: $\mathcal{S}_{l}\subseteq\mathcal{S}_{l+1}$, once a point $p$ appears as a node in $\mathcal{S}_{l}$, then every lower level in the tree has that node. 2. 2.

Covering: every node $q\in\mathcal{S}_{l+1}$ has a parent node $p\in\mathcal{S}_{l}$, where $|\!|p-q|\!|\leqslant\sigma 2^{-l}$. As a result, covering becomes finer at higher scales in a dyadic fashion. 3. 3.

Separation: nodes belonging to the same scale are separated by a minimal distance which dyadically shrinks at higher scales i.e. $\forall q,q^{\prime}\in\mathcal{S}_{l}$ we have $|\!|q-q^{\prime}|\!|>\sigma 2^{-l}$.

Each node $p$ also keeps the maximum distance to its descendants denoted by $\mathrm{maxdist}(q):=\max_{q^{\prime}\in\text{descendant}(q)}|\!|q-q^{\prime}|\!|\leqslant\sigma 2^{-l+1}$. which will be useful for fast NN/ANN search.222Note that any node $q\in\mathcal{S}_{l}$ due to the covering property satisfies $maxdist(q)\leqslant\sigma\left(2^{-l}+2^{-l-1}+2^{-l-2}+\dots\right)<\sigma 2^{-l+1}$ and thus, one might avoid saving $maxdist$ values and use this upper bound instead. The explicit tree representation requires the storage space $O(d)$ [9].

Definition 2**.**

Given a query $p$, $q\in\widetilde{\mathcal{C}}$ is a $(1+\varepsilon)$-ANN of $p$ for an $\varepsilon\geqslant 0$ if it holds: $|\!|p-q|\!|\leqslant(1+\varepsilon)\inf_{u\in\widetilde{\mathcal{C}}}|\!|p-u|\!|.$

Algorithm 1 details the branch-and-bound procedure for $(1+\varepsilon)$-ANN search on a given cover tree (proof of correctness is available in [9]). In short, we iteratively traverse down the cover tree and at each scale we populate the set of candidates $Q_{l}$ with nodes which could be the ancestors of the solution and discard others (this refinement uses the triangular inequality and the lower bound on the distance of the grandchildren of $Q$ to $p$, based on $\mathrm{maxdist}(q)$). At the finest scale (before stopping) we search the whole set of final candidates and report an $(1+\varepsilon)$-ANN point. Note that at each scale we only compute distances for non self-parent nodes (we pass, without computation, distance information of the self-parent children to finer scales). The case $\varepsilon=0$ refers to the exact tree NN search where one has to continue Algorithm 1 until the finest level of the tree. One should distinguish between this strategy and performing a brute-force search. Although they both perform an exact NN search, the complexity of Algorithm 1 is empirically shown to be way less in practical datasets.

Although the cover tree construction is blind to the explicit structure of data, several key growth properties such as the tree’s explicit depth, the number of children per node, and importantly the overall search complexity are characterized by the intrinsic doubling dimension of the model [10]. Practical datasets are often assumed to have small doubling dimensions e.g. when $\widetilde{\mathcal{C}}\subseteq\mathcal{M}$ samples a $K$-dimensional manifold $\mathcal{M}$ with certain regularity, one has $\dim_{D}(\widetilde{\mathcal{C}})\leqslant\dim_{D}(\mathcal{M})=O(K)$ [17]. The following theorem bounds the complexity of $(1+\varepsilon)$-ANN cover tree search [10, 9]:

Theorem 2**.**

Given a query which might not belong to $\widetilde{\mathcal{C}}$, the approximate $(1+\varepsilon)$-ANN search on a cover tree takes at most

[TABLE]

computations in time with $O(\#\widetilde{\mathcal{C}})$ memory requirement, where $\Delta$ is the aspect ratio333The aspect ratio of a set $\mathcal{C}$ is $\Delta:=\frac{\max|\!|q-q^{\prime}|\!|}{\min|\!|q-q^{\prime}|\!|}$, $\forall q\neq q^{\prime}\in\mathcal{C}$. of $\widetilde{\mathcal{C}}$ and $\#$ stands for set cardinality.

For most applications $\log(\Delta)=O(\log(d))$ [10] and thus for datasets with low dimensional structures i.e. $\dim_{D}=O(1)$ and by using moderate approximations one achieves a search complexity logarithmic in $d$, as opposed to the linear complexity of a brute-force search.

Note that the complexity of an exact cover tree search could be arbitrary high i.e. linear in $d$, at least in theory (unless the query belongs to the dataset [9] which does not generally apply e.g. for the intermediate steps of the IPG). However in our experiments we empirically observe that the complexity of an exact cover tree NN is still much lower than performing an exhaustive search.

4.1 Inexact data driven IPG algorithm

We use cover tree’s $(1+\varepsilon)$-ANN approximate search to efficiently implement the cone projection step of the data driven IPG algorithm. We only replace the exact search step (10) with the following approximation:

[TABLE]

We denote by $\mathcal{T}$ the cover tree structure built for the normalized dictionary with atoms $\psi_{i}/|\!|\psi_{i}|\!|$. For the current estimate we use the atom previously selected i.e. $X_{j}^{k}/\gamma_{j}$ for all $j$. We additionally normalize the gradient update $Z_{j}$ for the search input. This step may sound redundant however since all dictionary atoms are normalized and live on the hypersphere, we empirically observed that this trick leads to a better acceleration.

As a result by this update we obtain a $(1+\varepsilon)$-approximate cone projection denoted by $\mathcal{P}^{\varepsilon}_{\text{cone}(\widetilde{\mathcal{C}})}(Z_{j})$ satisfying definition (6). Therefore supported with an embedding assumption w.r.t. to the product model (9) and according to Theorem 1 we can guarantee stable CS recovery using this algorithm with a small computational effort.

5 Application in Quantitative MRI

In quantitative MRI rather than simply forming a single MR image with only contrast information (i.e. qualitative MRI), physicists are interested in measuring the NMR properties of tissues namely, the $T1,T2$ relaxation times, $\delta f$ magnetic resonance and proton density (PD), and use these parameters to differentiate different biological tissues [18]. The standard approach for parameter estimation is to acquire a large sequence of images in different times from a simple excitation pulse (i.e. excitations in terms of rotating the magnetic field) and then use an exponential model curve fitting to recover the decay exponents $T1,T2$ for each voxel. This procedure runs separately to estimate each relaxation time. The long process of acquiring multiple fully sampled images, brings serious limitation to standard approaches in quantitative imaging to apply within a reasonable time and with an acceptable signal-to-noise ratio (SNR) and resolution.

5.1 Magnetic Resonance Fingerprinting

Recently, a novel parameter estimation process coined as the Magnetic Resonance Fingerprinting (MRF) has been proposed to address this shortcoming. For the acquisition, MRF applies a shorter sequence of random excitations (i.e. $\Gamma\in[0,\pi/2]^{{\tilde{n}}}$ random flip angels with typically ${\tilde{n}}\approx 300-1000$ and the time interval $TR\sim 10-40$ msec) for estimating all parameters $T1,T2,\delta f$, and proton density at once. After each excitation pulse the response is recorded through the measurements taken from a small portion of k-space. As a result the acquisition time can significantly decreases to a couple of minutes [11].

The exponential model as for the traditional acquisition sequences will no longer hold for the MRF model. Due to using random excitations, the signatures casts a non-trivial implicit expression readable by solving the Bloch dynamic equations $\mathcal{B}(\Gamma,\Theta)\in\mathbb{C}^{\tilde{n}}$ customized for the set of parameters $\Theta=\{T1,T2,\delta f\}$. The MRF problem is an example of data driven CS model with an algebraic signal model (i.e. the Bloch equations) but a non-trivial projection. For this problem one constructs off-line a dictionary of fingerprints (i.e. magnetization responses) for all possible $d=\#\{\Theta_{i}\}_{i}$ parameters presented in normal tissues where the atoms $\psi_{i}:=\mathcal{B}(\Gamma,\Theta_{i})$ uniquely represent the underlying set of parameters $\Theta_{i}$ . This corresponds to sampling a low-dimensional manifold associated with the solutions of the Bloch dynamic equations [12], which for complicated excitation patterns neither the response nor the projection has an analytic solution.

5.2 MRF image model and parameter estimation

An MRF image can be represented by a matrix $X\in\mathbb{C}^{{\tilde{n}}\times J}$. Each column $X_{j}$ represents a ${\tilde{n}}$-dimensional spatial pixel $j$ and each row $X^{l}$ represents a temporal image slice (with the spatial resolution equal to $J$) corresponding to the measurements taken after the excitation pulse $\Gamma_{l}$, $\forall l=1,\ldots,{\tilde{n}}$. The total spatio-temporal resolution of the MRF image is $n:={\tilde{n}}J$. By the pixel purity assumption, the image follows the product space data driven model (7) and (8) where each column/pixel $j$ correspond to a single fingerprint $\psi$ with a unique parameter $\Theta$. The cone($\widetilde{\mathcal{C}}$) corresponds to the fingerprints dictionary and $\rho>0$ represents the proton density.444The proton density can be generally a complex number for which only the rescaling step of the cone projection slightly changes and the rest of the analysis similarly holds, for more details see [12]. The MRF acquisition model is linear as follows:

[TABLE]

where $\mathcal{F}$ denotes the two-dimensional Fourier transform, and $\mathcal{F}_{\Omega_{l}}:\mathbb{C}^{J}\rightarrow\mathbb{C}^{m/{\tilde{n}}}$ subsamples the k-space coefficients of the image slice $X^{l}$ according a pattern/set $\Omega_{l}$. The k-space subsampling could correspond to a spiral pattern [11] or a random horizontal/vertical line subselection [12] e.g. in the Echo Planar Imaging (EPI). By vectorization of $X,Y$, we have a product space CS model similar to (1) for a sampling protocol $A$ acting spatially however with different k-space sampling pattern per image slice.

The CS reconstruction proposed by [12] is based on a data driven IPG algorithm (with the exact cone projection described in Section 3) where one by identifying the correct fingerprint per pixel $X_{j}$ can recover the underlying parameters $\Theta$ stored in a look-up table (the last rescaling step also recovers the proton density $\gamma_{j}$). It has been shown in [12, Theorem 1] that if $\Omega_{l}$ independently at random subselects the rows (or columns) of the k-space, which corresponds to a random EPI sampling protocol, the resulting forward map $A$ is bi-Lipschitz w.r.t. the product model $\mathcal{C}=\prod_{j=1}^{J}\text{cone}(\widetilde{\mathcal{C}})$, provided sufficient discriminations between the fingerprints.555One achieves a suitable discrimination between the MR fingerprints by e.g. choosing the random excitation sequence long enough. This result guarantees stability of the parameter estimation using this algorithm with an step size $\mu=n/m$.

Since the Bloch manifold is parametrized by only three parameters in $\Theta$ and thus has a low intrinsic dimension, we use the cover tree ANN searches to accelerate the projection step. We perform the following procedure iteratively until convergence:

[TABLE]

The approximate cone projection follows the update in (13).

Provided with the same embedding result (for the random EPI acquisition protocol) and according to Theorem 1 we can deduce the linear convergence of this algorithm to the true solution for mildly-chosen approximations $\varepsilon$ (and the same step size $\mu=n/m$) where the projection step has now the complexity $O(J\log(d))$. For large-sized dictionaries the gap between the computation costs of the exact and inexact IPG algorithm can be significantly large.

6 Numerical Experiments

We investigate application of our proposed scheme for accelerating the MRF reconstruction. Parameter estimation consists of two off-line steps: i) forming the fingerprints dictionary by solving the Bloch dynamic equations for a wide range of parameters. ii) building a cover tree structure on the resulting dictionary which enables using fast ANN searches within the iterations of the inexact IPG algorithm.

We construct a cover tree on a dictionary composed of $d=48682$ fingerprints which are sampled from the Bloch response of $d$ pairs of $(T1,T2)$ relaxation times ranging between 100-5000 (ms) and 20-1800 (ms), respectively.666In our experiments we set the off-resonance frequency parameter to zero, the repetition time $TR=37$ (ms) and the echo time $TE=TR/2$. Fingerprints $\psi_{i}\in\mathbb{C}^{1024}$ are normalized magnetization responses of the corresponding parameters to a balanced steady-state free-precession (bSSFP or TrueFISP) sequence of 1024 slowly varying flip angles (excitations) ranging from $0^{\circ}$ to $60^{\circ}$, shown in Figure 1. From Figure 2 we can observe the MRF low dimensional manifold across its first three principal components. The hierarchical data partitions resulted by the cover tree is also depicted for four scales (the tree has total 14 scales with the finest coverage resolution $\sim 10^{-4}$). The ground truth image is a synthetic $256\times 256$ brain phantom with six $(T1,T2)$ segments previously used for evaluations in [12]. Figure 3 shows the corresponding segments, $T1$, $T2$ and proton density maps. This phantom is synthesized with the corresponding Bloch responses and proton density according to model (8) to build the magnetization image $X$. The k-space subsampling is based on a lattice-based non random EPI protocol (i.e. uniform row subselection with a shifted pattern for different image slices). We assume no noise on the measurements and we set the step size $\mu=n/m$ in all experiments.

We apply the $(1+\varepsilon)$-ANN tree search for different values of $\varepsilon=\{0,0.2,0.4,0.6,0.8\}$. We compare the performances of the exact and inexact data driven IPG on this synthetic phantom. We only report the cost of projections i.e. the total number of pairwise distances calculated for performing the NN or ANN searches.777Since the MRF’s forward/adjoint sampling operator is based on the Fast Fourier Transform, the gradient updates cost a tiny fraction of the search step (particularly when dealing with large-size datasets) and thus we exclude them from our evaluations. Results are depicted in Figure 4 and indicate that the inexact algorithm (for a suitable parameter $\varepsilon$) achieves a similar level of accuracy as for the exact IPG (with an exhaustive search) by saving 2-3 orders of magnitude in computations. Using an exact tree search (i.e. $\varepsilon=0$) also saves computations by an order of magnitude compared to using a brute-force search. Remarkably, for subsampling ratio $8:1$ the solution accuracy of the approximate scheme is about an order of magnitude better than the exact IPG. Although our theoretical results do not cover such observation, we shall relate it to a common practical knowledge that using soft decisions (e.g. here approximations) generally improves the performance of nonconvex algorithms compared to making hard decisions, and introduces a notion of robustness against undesirable local minima in such settings.

7 Conclusions

We proposed an inexact IPG algorithm for data driven CS problem with an application to fast MRF parameter estimation. The projection step is implemented approximately by using cover trees ANN searches. Thanks to the IPG’s robustness against approximations, by using this algorithm one achieves a comparable (and sometimes empirically enhanced) solution precision to the exact algorithm however with significantly less amount of computations regarding iteratively searching through large datasets.