Extended Version of "New Theory and Faster Computations for Subspace-Based Sensitivity Map Estimation in Multichannel MRI''

TL;DR

This paper introduces a new theoretical perspective on subspace-based sensitivity map estimation in multichannel MRI, specifically ESPIRiT, and proposes computational acceleration methods called PISCO that significantly reduce computation time and memory usage.

Contribution

It provides a novel theoretical derivation of sensitivity map estimation based on linear predictability and low-rank modeling, and introduces PISCO for fast, memory-efficient computations.

Findings

Theoretical derivation offers an intuitive understanding of ESPIRiT.

PISCO accelerates computations by up to 100 times.

Memory usage is substantially reduced with PISCO.

Abstract

This is an unabridged version of a journal manuscript that has been submitted for publication [1]. (Due to length restrictions, we were forced to remove substantial amounts of content from the version that was submitted to the journal, including more detailed theoretical explanations, additional figures, and a more comprehensive bibliography. This content remains intact in this version of the document). Sensitivity map estimation is important in many multichannel MRI applications. Subspace-based sensitivity map estimation methods like ESPIRiT are popular and perform well, though can be computationally expensive and their theoretical principles can be nontrivial to understand. In the first part of this work, we present a novel theoretical derivation of subspace-based sensitivity map estimation based on a linear-predictability/structured low-rank modeling perspective. This results in an estimation approach that is equivalent to ESPIRiT, but with distinct theory that may be more intuitive for some readers. In the second part of this work, we propose and evaluate a set of computational acceleration approaches (collectively known as PISCO) that can enable substantial improvements in computation time (up to ~100x in the examples we show) and memory for subspace-based sensitivity map estimation.

Tables

Table 1. TABLE I: Normalized Projection Residual with 32 × \times 32 calibration data.

Dataset	ESPIRiT	ESPIRiT	Nullspace	Nullspace
		+ PISCO		+ PISCO
Brain MPRAGE	0.136	0.141	0.136	0.141
Brain TSE	0.081	0.087	0.081	0.087
Knee TSE	0.297	0.298	0.297	0.298

Table 2. TABLE II: Computation Time (secs.) with 32 × \times 32 calibration data.

Dataset	ESPIRiT	ESPIRiT	Nullspace	Nullspace
		+ PISCO		+ PISCO
Brain MPRAGE	41.1	0.8	n/a	0.8
Brain TSE	97.4	0.9	n/a	0.9
Knee TSE	64.0	0.3	n/a	0.3

Table 3. TABLE III: Memory Usage (GB) with 32 × \times 32 calibration data.

Dataset	ESPIRiT	ESPIRiT	Nullspace	Nullspace
		+ PISCO		+ PISCO
Brain MPRAGE	1.8	0.1	47.2∗	0.1
Brain TSE	3.2	0.1	127.0∗	0.1
Knee TSE	3.5	0.1	35.2∗	0.1

Equations

$$d_q({\mathbf{k}}_m)=\int c_q(\mathbf{x})\rho (\mathbf{x})e^{-i2\pi {\mathbf{k}}_m^T\mathbf{x}}d\mathbf{x}+{\eta }_{qm}$$

$$\rho (\mathbf{x})=\text{mathbbm}{1}_{\Gamma }(\mathbf{x})\sum _{\mathbf{n}\in {\mathbb{Z}}^D}s[\mathbf{n}]e^{i2\pi {\mathbf{n}}^T\mathbf{x}}$$

$${\rho }_q(\mathbf{x})=\text{mathbbm}{1}_{\Gamma }(\mathbf{x})\sum _{\mathbf{n}\in {\mathbb{Z}}^D}{s}_q[\mathbf{n}]e^{i2\pi {\mathbf{n}}^T\mathbf{x}}$$

$$\mathbbm{1}_{\Gamma}(\mathbf{x})=\left\{\begin{array}[]{ll}1,&\mathbf{x}\in\Gamma\\ 0,&\mathrm{else}.\end{array}\right.$$

$$\sum _{q=1}^Q\sum _{\mathbf{m}\in \Lambda }{h}_r[\mathbf{m},q]s_q[\mathbf{n}-\mathbf{m}]\approx 0\text{ for }\forall \mathbf{n}\in {\mathbb{Z}}^D$$

$$\sum _{\mathbf{m}\in {\mathbb{Z}}^D}\tilde{f}[\mathbf{m}]s[\mathbf{n}-\mathbf{m}]=0\text{ for }\forall \mathbf{n}\in {\mathbb{Z}}^D.$$

$$\sum _{\mathbf{m}\in \Lambda }\tilde{f}[\mathbf{m}]s[\mathbf{n}-\mathbf{m}]\approx 0\text{ for }\forall \mathbf{n}\in {\mathbb{Z}}^D.$$

$$h_{r}[\mathbf{n},q]=\left\{\begin{array}[]{cl}\alpha_{q}\tilde{f}[\mathbf{n}],&\mathbf{n}\in\Lambda\\ 0,&\mathrm{else}\end{array}\right.$$

$$c_q(\mathbf{x})\approx \text{mathbbm}{1}_{\Gamma }(\mathbf{x})\sum _{\mathbf{n}\in \Lambda }{\tilde{c}}_q[\mathbf{n}]e^{i2\pi {\mathbf{n}}^T\mathbf{x}}$$

$$\begin{array}{rl}\sum _{\mathbf{m}\in \Lambda }& {\tilde{c}}_i[\mathbf{m}]s_j[\mathbf{n}-\mathbf{m}]\\ & -\sum _{\mathbf{m}\in \Lambda }{\tilde{c}}_j[\mathbf{m}]s_i[\mathbf{n}-\mathbf{m}]\approx 0\text{ for }\forall \mathbf{n}\in {\mathbb{Z}}^D.\end{array}$$

$$h_{r}[\mathbf{n},q]=\left\{\begin{array}[]{cl}-\tilde{c}_{j}[\mathbf{n}],&\mathbf{n}\in\Lambda,q=i\\ \tilde{c}_{i}[\mathbf{n}],&\mathbf{n}\in\Lambda,q=j\\ 0,&\mathrm{else}\end{array}\right.$$

$$\sum _{q=1}^Q{w}_q{\rho }_q(\mathbf{x})\approx 0\text{ for }\forall \mathbf{x}\in \Gamma .$$

$$h_{r}[\mathbf{n},q]=\left\{\begin{array}[]{ll}w_{q},&\mathbf{n}=\mathbf{0}\\ 0,&\text{else}.\end{array}\right.$$

$$\mathbf{C}\triangleq \left[\begin{array}{cccc}{\mathbf{C}}_1& {\mathbf{C}}_2& \cdots & {\mathbf{C}}_Q\end{array}\right];$$

$$\sum _{q=1}^Q{h}_r(\mathbf{x},q){\rho }_q(\mathbf{x})\approx 0\text{ for }\forall \mathbf{x}\in \Gamma ,$$

$$h_r(\mathbf{x},q)\triangleq \sum _{\mathbf{n}\in \Lambda }{h}_r[\mathbf{n},q]e^{i2\pi {\mathbf{n}}^T\mathbf{x}}.$$

$$\rho (\mathbf{x})\sum _{q=1}^Q{h}_r(\mathbf{x},q)c_q(\mathbf{x})\approx 0\text{ for }\forall \mathbf{x}\in \Gamma .$$

$$\sum _{q=1}^Q{h}_r(\mathbf{x},q)c_q(\mathbf{x})\approx 0\text{ for }\forall \mathbf{x}\in \Omega ,$$

$$\mathbf{H}(\mathbf{x})\mathbf{c}(\mathbf{x})\approx 0\text{ for }\forall \mathbf{x}\in \Omega ,$$

$$\mathbf{B}(\mathbf{x})={\mathbf{I}}_Q-\frac{1}{|\Lambda |}(\mathbf{H}(\mathbf{x}){)}^H\mathbf{H}(\mathbf{x}),$$

$${\mathbf{C}}^H\mathbf{C}=\left[\begin{array}{ccccc}{\mathbf{C}}_1^H{\mathbf{C}}_1& \cdots & {\mathbf{C}}_1^H{\mathbf{C}}_q& \cdots & {\mathbf{C}}_1^H{\mathbf{C}}_Q\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ {\mathbf{C}}_p^H{\mathbf{C}}_1& \cdots & {\mathbf{C}}_p^H{\mathbf{C}}_q& \cdots & {\mathbf{C}}_p^H{\mathbf{C}}_Q\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ {\mathbf{C}}_Q^H{\mathbf{C}}_1& \cdots & {\mathbf{C}}_Q^H{\mathbf{C}}_q& \cdots & {\mathbf{C}}_Q^H{\mathbf{C}}_Q\end{array}\right].$$

$$b[\mathbf{n}]=v[\mathbf{n}]⊛s_q[\mathbf{n}]\triangleq \sum _{\mathbf{m}\in \Lambda }v[\mathbf{m}]s_q[\mathbf{n}-\mathbf{m}]$$

$$\mathbf{G}(\mathbf{x})\in {\mathbb{C}}^{Q\times Q}\triangleq (\mathbf{H}(\mathbf{x}){)}^H\mathbf{H}(\mathbf{x}).$$

$$\begin{array}{rl}[\mathbf{G}(\mathbf{x}){]}_{pq}& =\sum _{r=1}^R{h}_r^{\ast }(\mathbf{x},p)h_r(\mathbf{x},q)\\ & =\sum _{s=1}^{|\Lambda |}\sum _{t=1}^{|\Lambda |}{e}^{i2\pi ({\mathbf{n}}_s-{\mathbf{m}}_t{)}^T\mathbf{x}}\left(\sum _{r=1}^R{h}_r^{\ast }[{\mathbf{m}}_t,p]h_r[{\mathbf{n}}_s,q]\right).\end{array}$$

$$\mathbf{W}=\left[\begin{array}{ccccc}{\mathbf{W}}_{11}& \cdots & {\mathbf{W}}_{1p}& \cdots & {\mathbf{W}}_{1Q}\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ {\mathbf{W}}_{q1}& \cdots & {\mathbf{W}}_{qp}& \cdots & {\mathbf{W}}_{qQ}\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ {\mathbf{W}}_{Q1}& \cdots & {\mathbf{W}}_{Qp}& \cdots & {\mathbf{W}}_{QQ}\end{array}\right].$$

$$[\mathbf{G}(\mathbf{x}){]}_{pq}=\sum _{\mathbf{m}\in {\mathbb{Z}}^D}\left(\sum _{s=1}^{|\Lambda |}{w}_{qp}^s[\mathbf{m}+{\mathbf{n}}_s]\right)e^{-i2\pi {\mathbf{m}}^T\mathbf{x}}.$$

Full text

Sensitivity map estimation is important in many multichannel MRI applications. Subspace-based sensitivity map estimation methods like ESPIRiT are popular and perform well, though can be computationally expensive and their theoretical principles can be nontrivial to understand. In the first part of this work, we present a novel theoretical derivation of subspace-based sensitivity map estimation based on a linear-predictability/structured low-rank modeling perspective. This results in an estimation approach that is equivalent to ESPIRiT, but with distinct theory that may be more intuitive for some readers. In the second part of this work, we propose and evaluate a set of computational acceleration approaches (collectively known as PISCO) that can enable substantial improvements in computation time (up to $\sim$100$\times$ in the examples we show) and memory for subspace-based sensitivity map estimation.

I Introduction

Most modern magnetic resonance imaging (MRI) experiments are performed using multichannel phased-array receiver coils [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]. Assuming $Q$ receiver channels, the measurements are usually modeled using the Fourier transform as

[TABLE]

for $m=1,\ldots,M$ and $q=1,\ldots,Q$. In this expression, $M$ denotes the number of k-space sampling locations $\mathbf{k}_{m}\in\mathbb{R}^{D}$, where $D$ is the dimension of the image;111For simplicity, we have adopted the use of “spatial” notation throughout this paper, which is consistent with typical use cases in 2D or 3D imaging. In these cases, we would have either $D=2$ with $\mathbf{x}=[x,y]^{T}$ and $\mathbf{k}=[k_{x},k_{y}]^{T}$ or $D=3$ with $\mathbf{x}=[x,y,z]^{T}$ and $\mathbf{k}=[k_{x},k_{y},k_{z}]^{T}$. However, it should be noted that our approach is also compatible with other scenarios such as spatiotemporal imaging with time-varying sensitivities, which, for example, might use $D=3$ with $\mathbf{x}=[x,y,f]^{T}$ and $\mathbf{k}=[k_{x},k_{y},t]^{T}$. $d_{q}(\mathbf{k}_{m})$ denotes the complex-valued k-space data measured from the $q$th receiver channel at the $m$th k-space location; $\rho(\mathbf{x})$ is the complex-valued underlying MR image (reflecting the state of the excited magnetization at the time of data acquisition) as a function of the spatial position $\mathbf{x}\in\mathbb{R}^{D}$; $c_{q}(\mathbf{x})$ is a complex-valued sensitivity map describing the sensitivity of the $q$th receiver channel to excited magnetization at spatial position $\mathbf{x}$; and $\eta_{qm}$ represents complex-valued noise. We use $\rho_{q}(\mathbf{x})\triangleq c_{q}(\mathbf{x})\rho(\mathbf{x})$ to denote the $q$th sensitivity-weighted image.

This paper is focused on the problem of estimating the sensitivity maps $c_{q}(\mathbf{x})$ from a set of calibration measurements (i.e., a set of k-space data samples that are sampled at the Nyquist rate). While prior knowledge of the coil sensitivity maps is not always required, there are many tasks that can benefit from estimated sensitivity maps, including coil combination (in which the underlying MRI image $\rho(\mathbf{x})$ is estimated from the multichannel images $\rho_{q}(\mathbf{x})$ [2, 19, 20]) and accelerated image reconstruction (in which the underlying MRI image $\rho(\mathbf{x})$ is estimated from multichannel k-space data $d_{q}(\mathbf{k}_{m})$ that is sampled below the Nyquist rate [9, 10, 11, 12, 16, 17, 18]).

Over the years, many different sensitivity map estimation methods have been proposed [9, 10, 19, 21, 22, 23, 24, 25, 26, 27, 28, 29, 20, 30]. Among these, subspace-based approaches [25, 27, 28, 26] have proven to be particularly powerful and popular, with the subspace-based ESPIRiT method [26] rising to become the most widely-used sensitivity map estimation method in the modern literature because of its simplicity (with few tuning parameters) and excellent performance.

The mathematical principles underlying subspace-based sensitivity map estimation are nontrivial and are not always easy to understand. In the first part of this work, we present a novel theoretical derivation of subspace-based sensitivity map estimation that is based on a nullspace/linear predictability perspective [31], and makes use of concepts from the recent literature on structured low-rank matrix modeling [15, 32, 33, 34, 35, 36, 37, 31, 38]. This approach is distinct from and complementary to existing theoretical explanations, and we personally find it to be more intuitive than the alternatives.

In the second part of this work, we introduce several new computational approaches that can substantially reduce the time and memory required for subspace-based sensitivity map estimation. This can be important because existing subspace-based methods can be computationally demanding, with substantial memory requirements and slow computation times. While slow computations can already be problematic when considering individual datasets, they can be especially limiting in modern machine learning contexts where it may be of interest to calculate sensitivity maps for every element of a database containing hundreds or thousands of images.

Our computational contributions include the following:

• •

We introduce a fast FFT-based approach that allows us to rapidly calculate the nullspace vectors associated with a multichannel convolution-structured matrix without directly constructing the matrix. Our approach can be viewed as a multichannel generalization of an approach that was previously developed for single-channel MRI reconstruction [39]. In addition to being useful for subspace-based sensitivity map estimation, this type of approach may also be useful for the broad class of structured low-rank modeling methods for multichannel MRI image reconstruction [15, 31, 33, 32, 34, 35, 36, 37, 38]. Note that a similar FFT-based approach was proposed for multichannel MRI reconstruction in Ref. [40]. Compared to that approach, our implementation uses substantially fewer FFTs, which is enabled by making use of structural redundancies.

• •

We leverage the previous observation that using ellipsoidal convolution kernels instead of conventional rectangular kernels can lead to substantial reductions in memory usage and computation without sacrificing accuracy [41].

• •

Subspace-based methods often require constructing a small matrix for every spatial location in the image, where the subspaces of the matrix reflect sensitivity information [26, 19]. We propose an efficient approach that uses vector-space simplifications and FFTs to directly calculate dimension-reduced versions of these matrices.

• •

Noting that sensitivity maps arise from Maxwell’s equations and must be smooth, we propose to estimate the sensitivity maps on a low-resolution spatial grid and interpolate to the desired spatial resolution.

• •

A common step in subspace-based methods is to compute the singular value decomposition (SVD) of a large number of small matrices (i.e., one SVD for each spatial location) [26, 19]. We introduce an efficient iterative approach (based on classical power iteration [42]) to calculate many small-scale partial SVDs simultaneously. While this approach is directly useful for subspace-based sensitivity map estimation, we expect that it may also have some utility for a range of locally low-rank modeling and denoising methods that also require calculating partial SVDs over local image patches for every spatial location in the image [43, 44, 45, 46, 47].

Combined together, these algorithmic improvements — which we collectively call PISCO (Power iteration over simultaneous patches, Interpolation, ellipSoidal kernels, and FFT-based COnvolution) — enable a substantial computational acceleration for subspace-based sensitivity map estimation.

This paper is organized as follows. Section II describes fundamental ambiguities in sensitivity map estimation that will be relevant in the sequel. Section III describes our novel derivation of subspace-based sensitivity map estimation from a linear predictability/structured low-rank perspective. Section IV describes and evaluates the PISCO approaches for fast computations. Finally, discussion and conclusions are presented in Sec. V.

II Ambiguities in the Inverse Problem

Before describing our novel theory and methods, we first review inherent ambiguities in sensitivity map estimation.

As already mentioned, we are interested in the estimation of sensitivity maps $c_{q}(\mathbf{x})$ for $q=1,\ldots,Q$ from a set of Nyquist-sampled calibration measurements. We assume that these calibration measurements consist of a collection of $M$ k-space samples from each of the $Q$ channels, for a total of $MQ$ observations. Our goal will be to calculate the values of the sensitivity maps on a discrete voxel grid of $N$ spatial locations $\mathbf{x}_{n}$ for $n=1,\ldots,N$, resulting in a total of $NQ$ unknowns. Most practical scenarios will involve substantially more spatial locations than k-space samples (i.e., $NQ\gg MQ$), such that the inverse problem is underdetermined (with fewer measurements than unknowns) unless additional assumptions are made about the characteristics of $c_{q}(\mathbf{x})$. To overcome this issue, the literature typically makes the assumption that the sensitivity maps $c_{q}(\mathbf{x})$ are spatially smooth, which can be justified by physical principles (e.g., Maxwell’s equations).

A more pernicious issue is that, while our goal is to estimate the unknown sensitivity maps $c_{q}(\mathbf{x})$, the image $\rho(\mathbf{x})$ that appears in Eq. (1) is also generally an unknown. This means that Eq. (1) can be viewed as a bilinear inverse problem, which is problematic because bilinear inverse problems have well-known scaling ambiguities. In particular, given one set of estimated sensitivity maps $c_{q}(\mathbf{x})$ for $q=1,\ldots,Q$ and a corresponding estimated image $\rho(\mathbf{x})$ that together provide a good fit to the measured data, it is straightforward to observe that another equally-good fit to the measured data can be obtained by combining a rescaled version of the sensitivity maps $\alpha(\mathbf{x})c_{q}(\mathbf{x})$ for $q=1,\ldots,Q$ together with a rescaled version of the image $\frac{1}{\alpha(\mathbf{x})}\rho(\mathbf{x})$. The complex-valued scaling function $\alpha(\mathbf{x})$ must be nonzero for every point $\mathbf{x}$ within the field of view (FOV), but is otherwise completely unconstrained, and can have arbitrary spatial variations. As such, unless additional strong assumptions are made, there will be infinitely many ways of choosing $\alpha(\mathbf{x})$ that all lead to sensitivity map estimates that provide equally-good fits to the measured data. Practically, this means that both the magnitude and phase of the estimated sensitivity maps are inherently ambiguous, and we can only hope for sensitivity map estimates that are unique up to these scaling ambiguities.

Because of the inherent ambiguity, it is common for sensitivity map estimation methods to search for one set of estimated maps $c_{q}(\mathbf{x})$ that fit the data well, and subsequently rescale them with a heuristically-chosen scaling function $\alpha(\mathbf{x})$ that can be designed to endow the scaled sensitivity maps (and the corresponding scaled image) with properties that are deemed to be useful. Note also that the sensitivity map values $c_{q}(\mathbf{x})$ will have no impact on measured data at spatial locations where $\rho(\mathbf{x})=0$. As a result, the values of the sensitivity maps are completely unidentifiable at spatial locations outside the support of $\rho(\mathbf{x})$, leading to additional ambiguity.

III Nullspace/Linear Predictability Theory for Subspace-Based Sensitivity Map Estimation

In this section, we present a novel derivation of subspace-based sensitivity map estimation from a nullspace/linear predictability perspective [31]. We start by reviewing linear predictability for multichannel MRI before making novel links to sensitivity map estimation.

III-A Summary of Linear Predictability for Multichannel MRI

Define the spatial support of $\rho(\mathbf{x})$ as the set of points $\Omega=\{\mathbf{x}\in\mathbb{R}^{D}:|\rho(\mathbf{x})|>0\}$. In what follows, we assume (without loss of generality) that the spatial coordinate system has been normalized such that the support $\Omega$ is completely contained within a hypercube $\Gamma$ with sides of length one, defined as $\Gamma\triangleq\{\mathbf{x}\in\mathbb{R}^{D}:\|\mathbf{x}\|_{\infty}<\frac{1}{2}\}$. We will refer to $\Gamma$ as the FOV and use $\Omega^{\complement}\triangleq\Gamma\setminus\Omega$ to denote the complement of $\Omega$ in $\Gamma$. Our FOV assumption allows us to equivalently represent the continuous images $\rho(\mathbf{x})$ and $\rho_{q}(\mathbf{x})$ using infinite Fourier series representations as

[TABLE]

and

[TABLE]

for $q=1,\ldots,Q$, where $s[\mathbf{n}]$ and $s_{q}[\mathbf{n}]$ respectively represent samples of the Fourier transforms of $\rho(\mathbf{x})$ and $\rho_{q}(\mathbf{x})$ on the rectilinear Nyquist grid (taking $\Delta k=1$ along each dimension, as enabled by our assumptions about $\Gamma$), and $\mathbbm{1}_{\Gamma}(\mathbf{x})$ is the indicator function for the FOV:

[TABLE]

A basic assumption of the linear predictability formalism [31] is that $\{s_{q}[\mathbf{n}]\}_{q=1}^{Q}$ will be autoregressive, satisfying multiple multichannel shift-invariant linear predictability relationships [34, 31, 35, 32, 36, 38, 14, 26, 25, 28, 27, 13, 15, 33, 48]. Specifically, there should exist multiple distinct $Q$-channel finite impulse response (FIR) filters $h_{r}[\mathbf{n},q]$ that each satisfy

[TABLE]

for $r=1,\ldots,R$, where $R$ is the number of distinct filters and $\Lambda\subset\mathbb{Z}^{D}$ is a finite index set describing the support of each FIR filter.

There are many different situations for which the existence of such filters $h_{r}[\mathbf{n},q]$ has been theoretically proven [31, 38]. We review a few of these below:222Shift-invariant linear predictability relationships can also be proven for images that possess smooth phase [33, 31, 49] and sparsity characteristics [31, 50, 51, 36, 37, 38], although these relationships are not directly relevant for sensitivity map estimation and we do not discuss them further.

III-A1 Limited Image Support [33, 48, 31]

In many situations, the actual spatial support $\Omega$ of the original image $\rho(\mathbf{x})$ will be smaller than the FOV $\Gamma$, such that the complement of $\Omega$ in $\Gamma$ (denoted by $\Omega^{\complement}\triangleq\Gamma\setminus\Omega$), is a measurable nonempty set. In such cases, it is possible to define nonzero functions $f(\mathbf{x})$ that have supports fully-contained within $\Omega^{\complement}$. Since $\rho(\mathbf{x})$ and $f(\mathbf{x})$ have nonoverlapping support, we must have $f(\mathbf{x})\rho(\mathbf{x})=0$ for $\forall\mathbf{x}\in\Gamma$. If we let $\tilde{f}[\mathbf{n}]$ denote the Nyquist-sampled Fourier series representation of $f(\mathbf{x})$, then applying the Fourier convolution theorem to the Fourier transform of $f(\mathbf{x})\rho(\mathbf{x})$ yields that

[TABLE]

If $\Omega^{\complement}$ is large enough, then it becomes possible to choose relatively smooth functions $f(\mathbf{x})$ that satisfy the conditions specified above. If we assume that these functions are smooth enough that they satisfy $\tilde{f}[\mathbf{n}]\approx 0$ for $\forall\mathbf{n}\notin\Lambda$ (i.e., the functions $f(\mathbf{x})$ are approximately bandlimited to the frequency range defined by $\Lambda$),333While we can often obtain very accurate approximations of this form, this must always be an approximation, and will never be exact. Specifically, since $f(\mathbf{x})$ is assumed to have finite support, then $\tilde{f}[\mathbf{n}]$ cannot have exactly-limited support by the Fourier uncertainty principle [52]. While this paper describes this approximation informally to simplify the exposition, more formal treatments exist in the literature [31]. then we can approximate $\tilde{f}[\mathbf{n}]$ as an FIR filter, and obtain

[TABLE]

Note also that if we have $f(\mathbf{x})\rho(\mathbf{x})\approx 0$ for $\forall\mathbf{x}\in\Gamma$, then the relationship between $\rho(\mathbf{x})$ and $\rho_{q}(\mathbf{x})$ implies we must also have that $f(\mathbf{x})\rho_{q}(\mathbf{x})\approx 0$ for $\forall\mathbf{x}\in\Gamma$ for $q=1,\ldots,Q$. This implies that the multichannel FIR filter

[TABLE]

will satisfy Eq. (5) for any choices of $\alpha_{q}\in\mathbb{C}$ for $q=1,\ldots,Q$. It is easy to obtain many distinct $h_{r}[\mathbf{n},q]$ filters from this relationship because of the flexibility offered by choosing different functions $f(\mathbf{x})$ and different coefficients $\alpha_{q}$. Using similar arguments, additional channel-specific filters can also be obtained in common situations where some sensitivity-weighted images $\rho_{q}(\mathbf{x})$ have smaller support than the underlying image $\rho(\mathbf{x})$.

III-A2 Smooth Sensitivity Maps [31, 38, 25, 28]

As already mentioned, sensitivity maps are often spatially smooth. Linear prediction relationships can be derived if we assume the $c_{q}(\mathbf{x})$ are so smooth that they are approximately bandlimited to the index set $\Lambda$, such that

[TABLE]

for appropriate Fourier coefficients $\tilde{c}_{q}[\mathbf{n}]$. To see why, observe that $c_{i}(\mathbf{x})\rho_{j}(\mathbf{x})-c_{j}(\mathbf{x})\rho_{i}(\mathbf{x})=0$ for $\forall\mathbf{x}\in\Gamma$ for any choices of $i$ and $j$ (sometimes called the “cross-relation” in the multichannel deconvolution literature [53, 54, 55]). Applying bandlimitedness and the Fourier convolution theorem to the cross-relation leads to

[TABLE]

This implies that the multichannel FIR filter

[TABLE]

will satisfy Eq. (5) for any choices of $i$ and $j$ with $i\neq j$. It is easy to obtain many distinct $h_{r}[\mathbf{n},q]$ filters using different choices of $i$ and $j$.

III-A3 Compressible Receiver Arrays [56, 57, 58]

It is frequently observed that the different channels in a multichannel receiver array possess linear dependencies, such that there exist multiple distinct choices of linear combination weights $w_{q}$ that will each satisfy

[TABLE]

This observation is often used for dimensionality reduction (“coil compression”) [56, 57, 58]. However, this observation also implies that a relationship in the form of Eq. (5) will hold if $h_{r}[\mathbf{n},q]$ is chosen as

[TABLE]

Since there often exist multiple different choices of $w_{q}$ that satisfy Eq. (12), there can also be many distinct filters $h_{r}[\mathbf{n},q]$ that arise from this kind of relationship.

The preceding arguments all support the existence of multiple filters $h_{r}[\mathbf{n},q]$ that satisfy Eq. (5). If we were given full information about the support $\Omega$ and the sensitivity maps $c_{q}(\mathbf{x})$, then we would be able to calculate filters $h_{r}[\mathbf{n},q]$ directly from Eqs. (8), (11), and (13). However, in practice, we do not have information about $\Omega$ or $c_{q}(\mathbf{x})$ (these are quantities that we are interested in estimating!), and we must instead resort to other methods to obtain the FIR filters $h_{r}[\mathbf{n},q]$.

An important observation from the literature on structured low-rank matrix modeling [31, 38] is that the convolution appearing in Eq. (5) can be represented in matrix form as $\mathbf{C}\mathbf{h}_{r}\approx\mathbf{0}$, where the matrix $\mathbf{C}$ is given by

[TABLE]

the matrix $\mathbf{C}_{q}\in\mathbb{C}^{P\times|\Lambda|}$ is a convolution-structured matrix (i.e., a Hankel or Toeplitz matrix) formed from Fourier samples $s_{q}[\mathbf{n}]$, with each row of $\mathbf{C}_{q}$ corresponding to one value of $\mathbf{n}$ from Eq. (5); $P$ represents the number of different $\mathbf{n}$ values that are used in the construction of each $\mathbf{C}_{q}$ matrix; and $\mathbf{h}_{r}\in\mathbb{C}^{Q|\Lambda|}$ is the vector of $h_{r}[\mathbf{n},q]$ samples. See Ref. [34] for a detailed description of the matrix construction we use.

This matrix representation shows that all filters $h_{r}[\mathbf{n},q]$ that satisfy Eq. (5) must be approximate nullspace vectors of an appropriately-constructed $\mathbf{C}$ matrix. Importantly, it is possible to directly form such a $\mathbf{C}$ matrix from Nyquist-sampled calibration data, and its approximate nullspace can be identified using the SVD. It is therefore possible to learn a set of $R$ linearly independent filters $h_{r}[\mathbf{n},q]$ (with $R$ equal to the dimension of the approximate nullspace of $\mathbf{C}$) in a data-driven way [15, 35]. The full set of approximate nullspace vectors is expected to form a basis for the set of all possible $h_{r}[\mathbf{n},q]$ filters that are consistent with Eq. (5) and bandlimited to $\Lambda$.

III-B Linear Predictability and Sensitivity Map Estimation

In what follows, we assume that we have access to $R$ filters $h_{r}[\mathbf{n},q]$ that each satisfy Eq. (5), which can, e.g., be obtained from the nullspace of a calibration matrix as described above.

Observe that the linear prediction relationship from Eq. (5) can be rewritten in the image domain as

[TABLE]

where

[TABLE]

Using the relationship between $\rho(\mathbf{x})$ and $\rho_{q}(\mathbf{x})$, Eq. (15) can be further simplified to

[TABLE]

Since $|\rho(\mathbf{x})|>0$ for $\forall\mathbf{x}\in\Omega$ (by the definition of the image support set $\Omega$), this implies that the summation

[TABLE]

demonstrating that the sensitivity maps $c_{q}(\mathbf{x})$ and the spatial-domain representation of the filters $h_{r}(\mathbf{x},q)$ have strong dependencies that must be satisfied for all $\mathbf{x}\in\Omega$. On the other hand, the summation in Eq. (18) could result in arbitrary nonzero values for $\mathbf{x}\in\Omega^{\complement}$, since the fact that $\rho(\mathbf{x})=0$ at these spatial locations will still ensure that Eq. (17) is satisfied regardless of the properties of $h_{r}(\mathbf{x},q)$ and $c_{q}(\mathbf{x})$.

Since we have $R$ filters $h_{r}[\mathbf{n},q]$ that each satisfy Eq. (5), the relationship from Eq. (18) can be expressed simultaneously for all $R$ filters in matrix form as

[TABLE]

where $\mathbf{H}(\mathbf{x})\in\mathbb{C}^{R\times Q}$ has entries $[\mathbf{H}(\mathbf{x})]_{rq}=h_{r}(\mathbf{x},q)$, and $\mathbf{c}(\mathbf{x})\in\mathbb{C}^{Q}$ is the vector with the $q$th entry equal to $c_{q}(\mathbf{x})$. Importantly, this demonstrates that the sensitivity map values at location $\mathbf{x}\in\Omega$ will be one of the approximate nullspace vectors of $\mathbf{H}(\mathbf{x})$. Ideally, there is only one approximate nullspace vector at the the spatial location $\mathbf{x}\in\Omega$. In this case, because singular vectors associated with distinct singular values are always unique up to scaling ambiguities, we would able to use the unique approximate nullspace vector to obtain unique sensitivity map values $c_{q}(\mathbf{x})$ up to the same scaling ambiguities that are inherent to sensitivity map estimation (see previous discussion in Sec. II).

The ability to estimate unique sensitivity maps $c_{q}(\mathbf{x})$ (up to scaling ambiguities) for $\mathbf{x}\in\Omega$ requires that the matrix $\mathbf{H}(\mathbf{x})$ is approximately rank $Q-1$. A necessary condition for this to be true is that $R\geq Q-1$, although this condition is not sufficient (and it is easy to find counter examples where $R\geq Q-1$ but $\mathbf{H}(\mathbf{x})$ does not have the desired characteristics). In practice, we have observed empirically that if we use the entire collection of filters from the approximate nullspace of $\mathbf{C}$,444Our experience is based on applying a heuristic thresholding rule to the SVD to determine the approximate nullspace. Specifically, if $\sigma_{n}$ denotes the $n$th singular value of $\mathbf{C}$ (with singular values indexed in order from largest to smallest), then the $n$th right singular vector is included in the approximate nullspace whenever $\sigma_{n}<0.05\sigma_{1}$. The $0.05\sigma_{1}$ threshold was determined heuristically and is unlikely to be optimal, although worked well for the datasets considered in this work. then the $\mathbf{H}(\mathbf{x})$ matrices always have an approximate rank of $Q-1$ for $\mathbf{x}\in\Omega$ under the assumptions and scenarios considered in this paper.555Different rank characteristics may be observed under different assumptions. For example, if the image support $\Omega$ were larger than the FOV $\Gamma$ such that aliasing occurs, then we empirically observe multiple sets of nullspace vectors at the spatial locations where the image support aliases onto itself. This is expected since such aliased voxels will be a linear combination of multiple sets of sensitivity-weighted images, each with distinct sensitivity maps. This enables the computation of multiple sensitivity maps, similar to what has been described previously in work on ESPIRiT [26], although we do not belabor this point as this scenario is not our primary interest. Interestingly, when $\mathbf{x}\in\Omega^{\complement}$ (i.e., the situation where the sensitivity maps are unidentifiable) and we use the entire collection of filters from the approximate nullspace of $\mathbf{C}$, we empirically observe that the matrix $\mathbf{H}(\mathbf{x})$ always has rank $Q$ rather than $Q-1$, which provides a simple mechanism for identifying both $\Omega$ and $\Omega^{\complement}$.

Empirical demonstrations of the nullspace characteristics of $\mathbf{C}$ and $\mathbf{H}(\mathbf{x})$ are presented in the next subsection.

III-C Empirical Demonstration of Nullspace Characteristics

Throughout this paper, we will use the following three multichannel MRI datasets to illustrate theoretical principles and to validate our proposed methods:

• •

Brain MPRAGE. We used an MPRAGE sequence on a 3T scanner to acquire a T1-weighted human brain image using a 32-channel receiver array, as depicted in Fig. 1(a).

• •

Brain TSE. We used a turbo spin-echo sequence on a 3T scanner to acquire a T2-weighted human brain image using a 32-channel receiver array, as depicted in Fig. 1(b).

• •

Knee TSE. We used a 15-channel knee dataset from the fastMRI database [59], as depicted in Fig. 1(c). This image was acquired using a fat-suppressed turbo spin-echo sequence on a 3T scanner, producing a proton-density weighted image.

To illustrate the low-rank characteristics of the $\mathbf{C}$ matrix as described in Sec. III-A, Fig. 2 plots the singular values of $\mathbf{C}$ for all three datasets. In each case, the $\mathbf{C}$ matrix was calculated using the central 24$\times$24 k-space samples from the Nyquist grid and assuming a $7\times 7$ rectangular FIR filter support (i.e., $\Lambda=\{\mathbf{n}\in\mathbb{Z}^{2}:\|\mathbf{n}\|_{\infty}\leq 3\}$). As expected, all three $\mathbf{C}$ matrices are approximately low-rank with a substantial approximate nullspace.

To illustrate the characteristics of the $\mathbf{H}(\mathbf{x})$ matrix described in Sec. III-B, Fig. 3 shows spatial maps of the $Q=32$ squared singular values of $\mathbf{H}(\mathbf{x})$ for the Brain TSE dataset. As expected, the smallest singular value is approximately zero within the support of the image (i.e., $\mathbf{H}(\mathbf{x})$ is approximately rank $=Q-1$ for $\mathbf{x}\in\Omega$), while $\mathbf{H}(\mathbf{x})$ has full rank (i.e., rank $=Q$) outside the image support.

III-D Naive Algorithm for Nullspace-Based Sensitivity Map Estimation

Based on the theoretical framework described in the previous subsections, a naive nullspace-based algorithm for calulating sensitivity maps would proceed as follows:

Choose an FIR filter support set $\Lambda$, and construct the corresponding $\mathbf{C}$ matrix from calibration data. 2. 2.

Calculate a set of $R$ approximate nullspace vectors of the $\mathbf{C}$ matrix using the SVD. 3. 3.

Using the nullspace vectors from the previous step, calculate the $\mathbf{H}(\mathbf{x})$ matrices for each spatial position $\mathbf{x}$. 4. 4.

For each spatial location $\mathbf{x}$, calculate the SVD of $\mathbf{H}(\mathbf{x})$ and set $c_{q}(\mathbf{x})$ equal to the $q$th entry of the right singular vector associated with the smallest singular value. 5. 5.

(Optional) If desired, choose a scaling function $\alpha(\mathbf{x})$ that has useful properties and use it to normalize the sensitivity maps. 6. 6.

(Optional) If desired, obtain a binary mask for the image support $\Omega$ by applying thresholding to the spatial map of the smallest singular value of $\mathbf{H}(\mathbf{x})$.

As we will demonstrate in the sequel, this naive algorithm works well although can require substantial computational complexity and memory usage, while a modified algorithm that uses our proposed PISCO approach leads to substantial computational accelerations and memory savings.

III-E Comparisons of the Nullspace Algorithm and ESPIRiT

Before moving on to describe the rest of our novel PISCO computational acceleration techniques, we will first spend a moment to describe how the nullspace-based approach compares to the popular ESPIRiT approach for sensitivity map estimation [26], which is also subspace-based. While the two approaches are derived using very different arguments, the resulting approaches end up having quite similar structure. Specifically, in our nullspace-based approach, we use the nullspace vectors of the $\mathbf{C}$ matrix to identify matrices $\mathbf{H}(\mathbf{x})$ that satisfy $\mathbf{H}(\mathbf{x})\mathbf{c}(\mathbf{x})\approx\mathbf{0}$, where $\mathbf{c}(\mathbf{x})\in\mathbb{C}^{Q}$ is the vector of $c_{q}(\mathbf{x})$ sensitivity map values. This causes us to estimate sensitivity maps using the approximate nullspace vector associated with the smallest singular value of $\mathbf{H}(\mathbf{x})$. In contrast, ESPIRiT uses the rowspace of the $\mathbf{C}$ matrix to construct matrices $\mathbf{B}(\mathbf{x})\in\mathbb{C}^{Q\times Q}$ that are expected to satisfy $\mathbf{B}(\mathbf{x})\mathbf{c}(\mathbf{x})\approx\mathbf{c}(\mathbf{x})$. As a result, the sensitivity maps in ESPIRiT can be estimated using the eigenvector associated with the largest eigenvalue of $\mathbf{B}(\mathbf{x})$. Interestingly, it can be derived (with a substantial amount of tedious-but-straightforward derivations that are too involved to reproduce here) that the $\mathbf{B}(\mathbf{x})$ matrix used by ESPIRiT will satisfy

[TABLE]

where $\mathbf{I}_{Q}$ denotes the $Q\times Q$ identity matrix. This implies that the results of the nullspace-based approach and ESPIRiT will be formally equivalent if the same parameters are used for both methods.

Our empirical results confirm that ESPIRiT and the nullspace-based approach produce identical sensitivity map estimates (up to numerical errors resulting from finite precision arithmetic). To illustrate this, we used in-house MATLAB implementations of both approaches666Note that the naive nullspace algorithm described in Sec. III-D can require an excessive amount of memory that exceeds MATLAB’s limits on the maximum array size. As such, the results we present as the nullspace-based approach in this section actually use one component of the PISCO approach to bypass this issue, namely the efficient calculation of a dimension-reduced matrix $\mathbf{G}(\mathbf{x})$ instead of the matrix $\mathbf{H}(\mathbf{x})$ as will be described in Sec. IV-C. to estimate sensitivity maps for all three datasets from Fig. 1 using different amounts of calibration data (ranging from small 24$\times$24 calibration regions to larger 256$\times$256 regions, always centered at the k-space origin). As a precursor to our description of PISCO in the sequel, we also obtained results from the nullspace-based approach combined with PISCO and ESPIRiT combined with PISCO. Representative qualitative results are illustrated in Fig. 4, which demonstrate negligible visual differences between all four approaches.

To be more quantitative, we have also computed a “normalized projection residual” metric to quantify how well the estimated sensitivity maps match the fully-sampled data. Specifically, we first take the fully-sampled multichannel data and project it onto the subspace spanned by the sensitivity maps [26].777It was necessary to use a specialized error metric because the inherent scaling ambiguities described in Sec. II confound traditional error metrics. If sensitivity map estimation worked well, then the projection onto the sensitivity maps should capture all of the true signal, with only noise remaining in the residual after projection. Our metric is obtained as the $\ell_{2}$-norm of the projection residual divided by the $\ell_{2}$-norm of the original multichannel data, with smaller values interpreted as better.

Figure 5 plots this metric for the Knee TSE dataset (with a fixed 5$\times$5 rectangular FIR filter support $\Lambda$) and the Brain TSE dataset (with a fixed 7$\times$7 rectangular FIR filter support $\Lambda$) for different calibration sizes, while Table I shows numerical results for all three datasets using 32$\times$32 calibration data. As before, the differences between ESPIRiT and the nullspace-based approach are negligible. (This table also shows PISCO results which are also quite similar, although we defer a detailed discussion until after we have described PISCO completely).

Overall, these results confirm that our nullspace-based approach produces sensitivity maps that are numerically equivalent to the sensitivity maps produced by ESPIRiT. Although this is expected theoretically, it is notable given that the two approaches were derived from distinct perspectives. This concludes the first major contribution of this paper.

It should be noted that both ESPIRiT and the nullspace-based approach can require substantial computation times (on the order of minutes) and memory (on the order of gigabytes), as reflected in the computation times and memory usage reported in Tables II and III, respectively.888Note that, as described previously, the naive nullspace-based approach was not directly implemented because we would exceed MATLAB’s limits on the maximum array size, so we do not report computation times and only report projected values (denoted with an asterisk) for memory usage. The remainder of this paper will focus on another major contribution: the PISCO computational approaches that are applicable to a range of subspaced-based approaches, including the nullspace-based approach and ESPIRiT. As can already be seen from Tables II and III, the use of PISCO can substantially improve computation time and memory usage relative to baseline implementations, and both the nullspace-based approach and ESPIRiT have nearly identical computational complexity when PISCO is used.

IV Proposed Computational Improvements

PISCO is comprised of five separate computational approaches, which are independently described and individually evaluated in the first five subsections of this section (one subsection per approach). Subsequently, we discuss the combined use of all five approaches in Sec. IV-F.

The following subsections (as well as the results already reported in the previous section) use common infrastructure to quantify the complexity of different computational approaches. All methods were implemented in-house using MATLAB, and all computations were performed on a desktop computer with an Intel Xeon E5-1603 2.8GHz quad-core CPU and 32GB RAM. Computation times were determined by calculating the elapsed time (in seconds) between the start and end of the computation (i.e., the wall-clock time). To minimize the potential for confounding effects, these computations were performed during quiescent periods (without other major software applications running in the foreground, and away from times when major background processes were scheduled to run). To reduce the influence of random background activity, we repeated each computation at least five times and report the median. We also evaluated the peak memory required for each implementation. To illustrate performance in different scenarios, we performed computations for the Knee TSE dataset (which had the smallest number of channels with $Q=15$) using a smaller 5$\times$5 FIR filter support and the Brain TSE dataset (which had a larger number of channels with $Q=32$) using a larger 7$\times$7 FIR filter support.

In our empirical evaluations, we compare our acceleration techniques against what we call the “Nullspace Baseline” (NB) approach. The NB corresponds to the naive nullspace-based algorithm from Sec. III-D combined with the memory-saving approach that will be described in Sec. IV-C.999As described in the previous section, we do not directly implement the naive nullspace method as the most straightforward implementation would exceed MATLAB’s limits on the maximum array size.

The different components of PISCO are described below.

IV-A FFT-Based Calculations to Determine the Nullspace of $\mathbf{C}$

Many subspace-based methods require constructing the $\mathbf{C}$ matrix and then calculating $\mathbf{C}^{H}\mathbf{C}$ and its approximate nullspace (cf. steps 1 and 2 of the naive nullspace-based algorithm from Sec. III-D). This calculation can sometimes be expensive because the $\mathbf{C}$ matrix can be quite large in scenarios where $P\gg Q|\Lambda|$. Generally speaking, the value of $P$ will be larger whenever the FIR filter support $\Lambda$ is smaller and/or whenever the size of the calibration region is larger, and the acceleration approach described in this subsection will be particularly beneficial in those circumstances.

Rather than working directly with $\mathbf{C}$ (and similar to previous work on structured low-rank MR image reconstruction [39, 40]), our approach is based on directly calculating the matrix $\mathbf{C}^{H}\mathbf{C}$ without ever calculating $\mathbf{C}$. This can be beneficial because $\mathbf{C}$ and $\mathbf{C}^{H}\mathbf{C}$ share the same nullspace structure, while $\mathbf{C}^{H}\mathbf{C}$ will have fewer entries than $\mathbf{C}$ whenever $P>Q|\Lambda|$.

An efficient FFT-based approach for directly calculating $\mathbf{C}^{H}\mathbf{C}$ was originally introduced in Ref. [39] for the single-channel case, which was later extendeded to the multichannel case in Ref. [40]. Our proposed approach follows the same basic principles though includes some improvements over the approach of Ref. [40], as described below.

Based on the definition of $\mathbf{C}$ from Eq. (14), the matrix $\mathbf{C}^{H}\mathbf{C}$ has $Q\times Q$ block structure, where the $(p,q)$th block is given by $\mathbf{C}_{p}^{H}\mathbf{C}_{q}$ for $p,q\in\{1,\ldots,Q\}$. This is

[TABLE]

Construction of $\mathbf{C}^{H}\mathbf{C}$ can thus be performed efficiently if we have efficient methods for computing $\mathbf{C}_{p}^{H}\mathbf{C}_{q}$.

Note that efficient methods for computing $\mathbf{C}_{p}^{H}\mathbf{C}_{q}$ can potentially be obtained if we have efficient methods for computing matrix-vector multiplications of the form $\mathbf{C}_{p}^{H}\mathbf{C}_{q}\mathbf{v}$ for arbitrary vectors $\mathbf{v}\in\mathbb{C}^{|\Lambda|}$. Specifically, the $i$th column of the matrix $\mathbf{C}_{p}^{H}\mathbf{C}_{q}$ can be computed using the matrix-vector multiplication $\mathbf{C}_{p}^{H}\mathbf{C}_{q}\mathbf{e}_{i}$, where $\mathbf{e}_{i}$ is the $i$th column of the identity matrix $\mathbf{I}_{|\Lambda|}$. This allows the entire $\mathbf{C}^{H}\mathbf{C}$ matrix to be calculated by calculating $\mathbf{C}_{p}^{H}\mathbf{C}_{q}\mathbf{e}_{i}$ for each $(p,q,i)$ combination with $p,q\in\{1,\ldots,Q\}$, and $i\in\{1,\ldots,|\Lambda|\}$.

The convolution structure of the $\mathbf{C}_{q}$ matrices implies that, for any vector $\mathbf{v}\in\mathbb{C}^{|\Lambda|}$, the matrix-vector multiplication $\mathbf{C}_{q}\mathbf{v}$ can be expressed as first performing the simple convolution

[TABLE]

(where $v[\mathbf{n}]$ is a simple reformatting of $\mathbf{v}$ and $\circledast$ denotes convolution) and then extracting the $P$ samples of $b[\mathbf{n}]$ corresponding to the $P$ values of $\mathbf{n}$ that were used in the construction of $\mathbf{C}$ (which can be viewed as masking the convolution output). Similarly, for any vector $\mathbf{u}\in\mathbb{C}^{P}$, the matrix-vector multiplication $\mathbf{C}_{p}^{H}\mathbf{u}$ can be expressed as reformatting the vector $\mathbf{u}$ into a sequence $u[\mathbf{n}]$, computing $u[\mathbf{n}]\circledast s_{p}^{*}[-\mathbf{n}]$, and finally extracting relevant samples from the convolution output, where ∗ denotes complex conjugation.

As described by Ongie and Jacob [39], the separate convolution operations that can be used to calculate $\mathbf{C}_{q}\mathbf{v}$ and $\mathbf{C}^{H}_{p}\mathbf{u}$ can each be efficiently computed (without approximation) using standard FFT-based methods. Moreover, if the masking operations are neglected (which typically only has a minor effect on the final results [39]), then the composite matrix-vector multiplication $\mathbf{C}_{p}^{H}\mathbf{C}_{q}\mathbf{v}$ can be accurately approximated by calculating $v[\mathbf{n}]\circledast(s_{q}[\mathbf{n}]\circledast s_{p}^{*}[-\mathbf{n}])$ and extracting relevant samples from the convolution output. If zero-padded FFTs of $s_{q}[\mathbf{n}]$ and $s_{p}[\mathbf{n}]$ have been precomputed, then this convolution can be easily computed by taking the zero-padded FFT of $v[\mathbf{n}]$, multiplying the result with the precomputed FFTs of $s_{q}[\mathbf{n}]$ and $s_{p}^{*}[-\mathbf{n}]$,101010Note that the FFT of $s_{p}^{*}[-\mathbf{n}]$ can be easily retrieved from the FFT of $s_{p}[\mathbf{n}]$ using the conjugation property of the Fourier transform. and then taking an inverse FFT.

The steps described above closely resemble the multichannel approach of Zhao et al. [40]. In our implementation (and similar to the original single-channel work of Ongie and Jacob [39]), we additionally use the fact that if $\mathbf{v}=\mathbf{e}_{i}$, then the corresponding sequence representation $v[\mathbf{n}]$ will be equal to $\delta[\mathbf{n}-\mathbf{n}_{i}]$, where $\delta[\cdot]$ is the Kronecker delta function and $\mathbf{n}_{i}\in\mathbb{Z}^{D}$ is the $i$th element of $\Lambda$. As a result, the matrix-vector multiplications $\mathbf{C}_{p}^{H}\mathbf{C}_{q}\mathbf{e}_{i}$ for different choices of $i$ can be simply approximated as convolutions of the form $\delta[\mathbf{n}-\mathbf{n}_{i}]\circledast(s_{q}[\mathbf{n}]\circledast s_{p}^{*}[-\mathbf{n}])$ with different choices of $\mathbf{n}_{i}$. Importantly, the sifting property of the Kronecker delta allows us to retrieve the results of $\delta[\mathbf{n}-\mathbf{n}_{i}]\circledast(s_{q}[\mathbf{n}]\circledast s_{p}^{*}[-\mathbf{n}])$ for every choice of $\mathbf{n}_{i}$ by simply shifting the output of $\delta[\mathbf{n}]\circledast(s_{q}[\mathbf{n}]\circledast s_{p}^{*}[-\mathbf{n}])=s_{q}[\mathbf{n}]\circledast s_{p}^{*}[-\mathbf{n}]$ by a number of samples according to $\mathbf{n}_{i}$. This means that it is not necessary to separately calculate distinct convolution results for each $i\in\{1,\ldots,|\Lambda|\}$. Instead, it suffices to calculate a single convolution $s_{q}[\mathbf{n}]\circledast s_{p}^{*}[-\mathbf{n}]$ per $(p,q)$ pair and then construct approximations of the columns of $\mathbf{C}_{p}^{H}\mathbf{C}_{q}$ by simple shifts of the convolution output!

A close approximation of $\mathbf{C}^{H}\mathbf{C}$ can thus be obtained by precomputing zero-padded FFTs of each $s_{q}[\mathbf{n}]$ for $q=1,\ldots,Q$ (a total of $Q$ FFTs) and then multiplying the results together and computing inverse FFTs for each $(p,q)$ pair (a total of $Q^{2}$ additional FFTs). Once this approximation of $\mathbf{C}^{H}\mathbf{C}$ is obtained, it is straightforward to calculate its approximate nullspace using standard SVD algorithms.

The computation associated with this acceleration approach scales with $Q^{2}$, so is expected to be more advantageous in situations where the number of channels is small, and may not be beneficial in situations where $Q$ is large. And as already mentioned, this approach is expected to be more beneficial in situations where $P>Q|\Lambda|$, and is less likely to be beneficial in situations where $P<Q|\Lambda|$.

Our empirical results for this acceleration approach are shown in Fig. 6, and are consistent with our theoretical expectations. While our FFT-based approach does involve a small amount of approximation, the effects of this approximation are observed to be relatively minor, and the normalized projection residual is not substantially changed. On the other hand, we see substantial (though somewhat nuanced) changes in computation time. For large calibration regions (corresponding to situations where $P\gg Q|\Lambda|$), we consistently observe that the convolutional acceleration outperforms the Nullspace Baseline. However, when the calibration region is small, the convolutional acceleration is only better than the Nullspace Baseline for the Knee TSE data, but not for the Brain TSE data. This is consistent with expectations, as the Knee TSE data has both larger $P$ than the Brain TSE data (due to the use of a smaller FIR filter support $\Lambda$) and smaller $Q$. While there are subtle differences in memory usage, these differences were not substantial and are not shown.

IV-B *Ellipsoidal FIR Filter Support *

When constructing the $\mathbf{C}$ matrix, methods like ESPIRiT have traditionally chosen the FIR filter support $\Lambda$ to have a rectangular shape (e.g., in the 2D case, a $(2\tau+1)$ ​$\times$ ​$(2\tau+1)$ rectangular filter support is obtained when using $\Lambda=\{\mathbf{n}\in\mathbb{Z}^{2}:\|\mathbf{n}\|_{\infty}\leq\tau\}$). However, previous work has shown that the corners of the rectangular filter have a negligible impact for methods that rely on the $\mathbf{C}$ matrix [41], and that removing these unimportant corners (i.e., simplifying the rectangle into an ellipsoid using $\Lambda=\{\mathbf{n}\in\mathbb{Z}^{2}:\|\mathbf{n}\|_{2}\leq\tau\}$) offers substantial computational advantages with negligible impact on the final results. This same approach is adopted in PISCO.

Our empirical results for this acceleration approach are shown in Fig. 7. Although this acceleration approach does slightly change the final sensitivity maps compared to the Nullspace Baseline, the differences are not substantial as can be observed from the normalized projection residual. However, the acceleration approach does provide a slight but noticeable benefit to the total computation time when calculating sensitivity maps, particularly as the size of the calibration region grows. Although this acceleration approach also results in memory benefits as previously described [41], these benefits were not very noticeable (and are not shown) since the steps involving the $\mathbf{C}$ matrix are not the most memory-consuming steps of the Nullspace Baseline.

IV-C Replacing $\mathbf{H}(\mathbf{x})$ with $\mathbf{G}(\mathbf{x})$

As already noted, constructing $\mathbf{H}(\mathbf{x})$ (i.e., constructing one $R\times Q$ matrix for every one of the $N$ spatial locations in the image) can require a substantial amount of memory, and can impede implementation of the naive nullspace-based algorithm. While there are simple ways to bypass this problem (i.e., only computing $\mathbf{H}(\mathbf{x})$ for one voxel at a time instead of storing the matrix values for all voxels simultaneously), this can result in a loop over all voxels that can be very slow in the common case where the number of voxels is large.

Our proposed approach is based on replacing $\mathbf{H}(\mathbf{x})$ with

[TABLE]

Similar to the arguments we made for preferring $\mathbf{C}^{H}\mathbf{C}$ over $\mathbf{C}$ in Sec. III-A, the matrix $\mathbf{G}(\mathbf{x})\in\mathbb{C}^{Q\times Q}$ shares the same nullspace structure with $\mathbf{H}(\mathbf{x})\in\mathbb{C}^{R\times Q}$, while $\mathbf{G}(\mathbf{x})$ will be substantially smaller when the number of nullspace filters $R$ is much larger than the number of channels $Q$ (which is a common occurence as can be seen from the large nullspaces observed in Fig. 2). Directly calculating $\mathbf{G}(\mathbf{x})$ instead of $\mathbf{H}(\mathbf{x})$ is thus expected to reduce memory utilization by at least a factor of $R/Q$, which is often substantial (e.g., the value of $R/Q$ ranges from between 20 to 50 for the datasets and parameter settings considered in this paper!).

A further observation that underlies our proposed approach is that $\mathbf{G}(\mathbf{x})$ is structured, with entries that can be written as

[TABLE]

A naive calculation of $\mathbf{G}(\mathbf{x})$ at every spatial location using this formula would use $Q^{2}|\Lambda|^{2}(R+1)N$ multiplications, which could be quite expensive.

To proceed further, let $\mathbf{W}\in\mathbb{C}^{Q|\Lambda|\times Q|\Lambda|}$ be defined as $\mathbf{W}=\sum_{r=1}^{R}\mathbf{h}_{r}\mathbf{h}_{r}^{H}$, where $\mathbf{h}_{r}$ is the $r$th approximate nullspace vector of $\mathbf{C}$ as previously defined in Sec. III-A. Note that $\mathbf{W}$ is easily calculated using $Q^{2}|\Lambda|^{2}R$ multiplications. It is straightforward to observe that $\mathbf{W}$ has a $Q\times Q$ block structure, where we denote the $(q,p)$th block by $\mathbf{W}_{qp}\in\mathbb{C}^{|\Lambda|\times\Lambda|}$ for $p,q\in\{1,\ldots,Q\}$. Specifically,

[TABLE]

The entry in the $s$th row and $t$th column of $\mathbf{W}_{qp}$ is given by $[\mathbf{W}_{qp}]_{st}=\sum_{r=1}^{R}h_{r}^{*}[\mathbf{m}_{t},p]h_{r}[\mathbf{n}_{s},q]$, which corresponds to one of the terms appearing in Eq. (24). Using $w_{qp}^{s}[\mathbf{m}_{t}]$ to denote the $s$th row of $\mathbf{W}_{qp}$ reformatted as a sequence and leveraging the shift property of the discrete Fourier transform, Eq. (24) simplifies to

[TABLE]

Importantly, the inner summation over $s$ does not require any additional multiplications beyond the $Q^{2}|\Lambda|^{2}R$ multiplications that were already used in the formation of $\mathbf{W}$, while the outer summation can be easily evaluated for every spatial location using a zero-padded FFT ($O(N\log_{2}N)$ multiplications). This calculation needs to be repeated for every $(p,q)$ combination, such that the total number of complex multiplications required for our proposed approach is thus $O(Q^{2}(N\log_{2}N+|\Lambda|^{2}R))$. This can be substantially smaller than the $Q^{2}|\Lambda|^{2}(R+1)N$ multiplications needed for the naive calculation of $\mathbf{G}(\mathbf{x})$.

Since we could not implement the naive nullspace-based algorithm using $\mathbf{H}(\mathbf{x})$, our new ability to calculate $\mathbf{G}(\mathbf{x})$ is an enabling technology for the nullspace-based algorithm. Figure 8 shows a comparison between the memory utilized by our proposed $\mathbf{G}(\mathbf{x})$ approach compared against our projections of the memory required for $\mathbf{H}(\mathbf{x})$. As expected, there are orders-of-magnitude improvements in required memory.

IV-D Smoothness-Based Interpolation

The fact that subspace-based methods require manipulating matrices like $\mathbf{H}(\mathbf{x})$ or $\mathbf{G}(\mathbf{x})$ at every spatial location $\mathbf{x}$ can also be burdensome if the number of spatial locations $N$ is large. Using the well-known smoothness characteristics of sensitivity maps, we observe that instead of estimating sensitivity maps directly for the nominal high-resolution voxel grid, it should be sufficient to compute $\mathbf{G}(\mathbf{x})$ and estimate sensitivity maps on a low-resolution voxel grid, and then interpolate back to the desired spatial resolution. Our implementation relies on periodic sinc interpolation, which can be implemented efficiently using zero-padded FFTs. For the implementation reported in this work, we used a low-resolution voxel grid that was heuristically chosen to be 24 samples bigger than the calibration region along each dimension (i.e., for a 24$\times$24 calibration region, we used a 48$\times$48 low-resolution voxel grid). This choice worked well for the datasets considered in this paper, although alternative heuristics may lead to even better results, and the optimal selection of interpolation factors is likely to be application-dependent.

For interpolation to be successful, it is necessary for the low-resolution sensitivity map estimates to be smooth, which is potentially confounded by the ambiguities described in Sec. II, including the inherent phase-ambiguities associated with the SVD of $\mathbf{G}(\mathbf{x})$ which can result in quasi-random non-smooth phase characteristics. To avoid this issue, we choose a scaling function $\alpha(\mathbf{x})$ to provide smooth phase characteristics. Specifically, we use the initial sensitivity map estimates (which often have arbitrary nonsmooth phase) to perform coil-combination of low-resolution (Gaussian-apodized) Fourier reconstructions of the calibration data, and then choose $\alpha(\mathbf{x})$ to cancel the phase of the coil-combined reconstruction (so that the coil-combined reconstruction becomes real-valued). This approach is based on the assumption that the underlying image $\rho(\mathbf{x})$ does not have substantial phase variations such that the observed phase can be largely be ascribed to the sensitivity maps, which are presumed smooth. This approach worked well for the datasets considered in this paper, although more care may be needed in cases where the image phase is expected to be more complicated.

Our empirical results for interpolation are reported in Fig. 9, which demonstrates that interpolation leads to substantial reductions in the total computation time and memory utilization with only minor effects on the normalized projection residual.

IV-E Power Iteration

One of the most expensive steps of a subspace-based approach is the calculation of SVDs of $\mathbf{H}(\mathbf{x})$ or $\mathbf{G}(\mathbf{x})$ at every spatial location. Rather than performing full SVDs, we propose to use iterative methods to estimate only the right singular vector associated with the smallest singular value. For example, this can be done using classical power iteration [42], and it is also straightforward to modify the iterative procedure to perform this estimation process for all spatial locations simultaneously in parallel using vectorization techniques, thereby avoiding a loop over spatial locations.

Fast convergence of power iteration is expected when there is a big difference in magnitude between the smallest two singular values of $\mathbf{G}(\mathbf{x})$ [42], which we have observed is usually the case (cf. Fig. 3). As a result, we only need to use a few iterations to achieve suitable convergence. However, in cases where multiple sensitivity maps are required to describe a specific spatial location (e.g., when aliasing occurs within the FOV or in the presence of substantial partial voluming), then the smallest two singular values may have similar magnitudes and power iteration may be slower to converge. In such cases, the use of more advanced iterative algorithms (e.g., the Lanczos algorithm) may be preferable.

Our empirical results for power iteration (which used a fixed number of 10 iterations) are shown in Fig. 10, and demonstrate substantial improvements in the total computation time relative to the Nullspace Baseline with only minor changes in the normalized projection residual. The use of power iteration did not have a major impact on memory utilization (not shown).

IV-F All Components of PISCO Combined

The previous subsections demonstrated that each of the individual components of PISCO can separately provide computational advantages. In this section, we demonstrate that the combination of all five PISCO approaches leads to even further improvements. Note that while some acceleration techniques did not provide improvements in every situation (e.g., the FFT-based calculation of the nullspace of $\mathbf{C}$ from Sec. III-A), we still always used all five techniques together for simplicity. Even better performance would have been obtained if we selectively turned off acceleration techniques in regimes where they were not useful.

Figure 11 compares the results of the Nullspace Baseline against the nullspace-based method with PISCO. In all scenarios, we observe consistently substantial improvements in computation time (up to $\sim$100$\times$ acceleration in some cases) and substantial improvements in memory usage (even though the Nullspace Baseline already makes use of our $\mathbf{G}(\mathbf{x})$ computation, which is one of the most memory-saving components of PISCO). Figure 12 shows representative sensitivity maps, which qualitatively demonstrates the strong similarity between the sensitivity maps obtained with both approaches. Detailed numbers for 32$\times$32 calibration data were already presented in Tables I-III and support the same conclusions.

V Discussion and Conclusions

This work presented a new theoretical framework for subspace-based sensitivity map estimation based on principles from the literature on linear predictability and structured low-rank modeling. This led to a nullspace-based sensitivity map estimation method that is mathematically equivalent to ESPIRiT, but which has a distinct and complementary theoretical explanation that will be potentially more intuitive to some readers. In addition, we introduced the set of PISCO computational accelerations for subspace-based sensitivity map estimation that enable substantial reductions in computation time (up to $\sim$100$\times$ in some of our examples) without a substantial loss in estimation quality. The PISCO class of techniques is not only useful for our nullspace-based algorithm, but also offers substantial improvements for other subspace-based sensitivity map estimation methods like ESPIRiT. While we have not demonstrated further extensions to applications beyond sensitivity map estimation, we expect that some of the PISCO techniques may also enable computational improvements for other types of similar computations, including subspace-based adaptive coil combination [19], structured low-rank image reconstruction [15, 32, 33, 34, 35, 36, 37, 31, 38], and other patch-based low-rank modeling methods that require the simultaneous SVDs of a large number of small-scale matrices [43, 46, 47].

One of the optional steps given in the nullspace-based algorithm in Sect. III-D, corresponds to choosing a scaling function $\alpha(\mathbf{x})$ with useful properties to normalize the estimated sensitivity maps. In the experiments shown in this work we have chosen this scaling function such that the magnitude of the estimated sensitivity maps is normalized following a sum-of-squares approach, and the phase is referenced with respect to the phase of one of the estimated sensitivity maps [26]. However, this type of scaling function is not necessarily optimal, and different applications can be benefited from other normalization approaches.

It should be noted that the results presented in this paper were all based on in-house MATLAB implementations. We do not expect MATLAB computation times to be competitive with computation times achieved using software written in programming languages like C, and expect that computational improvements will easily be achieved by implementing PISCO in another programming language. Interestingly, we have anecdotally compared a C-based implementation of ESPIRiT (i.e., the implementation available from the BART software package [60]) against our PISCO-based MATLAB implementation of ESPIRiT using the Brain TSE dataset with 24$\times$24 calibration data. While the MATLAB implementation is naturally expected to have a major disadvantage in this comparison, it is notable that the PISCO-based MATLAB computation completed in $\sim$1 second, which was substantially faster than the $\sim$19 seconds used by the C-based implementation. This underscores the promise of PISCO, and strongly suggests that a C-based implementation of PISCO might be worthwhile.

The authors have made an open-source implementation of the nullspace-based algorithm with PISCO available at http://mr.usc.edu/download/pisco/.