# Tuning-free multi-coil compressed sensing MRI with Parallel Variable Density
Approximate Message Passing (P-VDAMP)

Charles Millard  Mark Chiew  Jared Tanner  Aaron T Hess and Boris Mailhé We
are grateful to Thijs de Buck for sharing his angiogram data. This work was
supported in part by an EPSRC Industrial CASE studentship with Siemens
Healthineers under voucher 17000051, in part by The Alan Turing Institute
under EPSRC grant EP/N510129/1, in part by EPSRC under grant EP/T013133/1, in
part by the Royal Academy of Engineering under grant RF201617/16/23, and in
part by the Wellcome Trust under grant 203139/Z/16/Z. The concepts and
information presented in this paper are based on research results that are not
commercially available.Charles Millard, Mark Chiew and Aaron T. Hess are with
the Wellcome Centre for Integrative Neuroimaging, FMRIB, Nuffield Department
of Clinical Neurosciences, University of Oxford, Level 0, John Radcliffe
Hospital, Oxford, OX3 9DU, UK. (email<EMAIL_ADDRESS><EMAIL_ADDRESS>and<EMAIL_ADDRESS>Jared Tanner is with
the Mathematical Institute at the University of Oxford, Oxford, OX2 6GG, UK
(email: tanner@maths.ox.ac.uk)Boris Mailhé is with Siemens Healthineers, 755
College Rd E, Princeton, NJ 08540, USA (email: boris.mailhe@siemens-
healthineers.com)

###### Abstract

Magnetic Resonance Imaging (MRI) has excellent soft tissue contrast but is
hindered by an inherently slow data acquisition process. Compressed sensing,
which reconstructs sparse signals from incoherently sampled data, has been
widely applied to accelerate MRI acquisitions. Compressed sensing MRI requires
one or more model parameters to be tuned, which is usually done by hand,
giving sub-optimal tuning in general. To address this issue, we build on
previous work by the authors on the single-coil Variable Density Approximate
Message Passing (VDAMP) algorithm, extending the framework to multiple
receiver coils to propose the Parallel VDAMP (P-VDAMP) algorithm. For
Bernoulli random variable density sampling, P-VDAMP obeys a “state evolution”,
where the intermediate per-iteration image estimate is distributed according
to the ground truth corrupted by a zero-mean Gaussian vector with
approximately known covariance. To our knowledge, P-VDAMP is the first
algorithm for multi-coil MRI data that obeys a state evolution with accurately
tracked parameters. We leverage state evolution to automatically tune sparse
parameters on-the-fly with Stein’s Unbiased Risk Estimate (SURE). P-VDAMP is
evaluated on brain, knee and angiogram datasets and compared with four
variants of the Fast Iterative Shrinkage-Thresholding algorithm (FISTA),
including two tuning-free variants from the literature. The proposed method is
found to have a similar reconstruction quality and time to convergence as
FISTA with an optimally tuned sparse weighting and offers substantial
robustness and reconstruction quality improvements over competing tuning-free
methods.

Magnetic Resonance Imaging (MRI), Compressed Sensing, Parallel Imaging,
Approximate Message Passing, Stein’s Unbiased Risk Estimate (SURE)

##  1 Introduction

Reducing acquisition time is one of the main challenges for rapid, dynamic
Magnetic Resonance Imaging (MRI). Accordingly, there is considerable research
attention on techniques that reconstruct sub-sampled k-space data without
substantially sacrificing image quality. Parallel imaging exploits the non-
uniformity of coil sensitivites to reconstruct data that have been
simultaneously acquired over a multichannel array of receiver coils [1, 2, 3,
4]. More recently, compressed sensing [5, 6] has been applied to MRI to
accommodate further acceleration [7, 8], using presumed redundancy in the
signal and regularized optimization techniques to reconstruct incoherently
sampled data.

A frequently noted challenge for compressed sensing MRI (CS-MRI) is how to
tune model parameters [9, 10, 11], which is typically done by hand [12]. The
ideal parameter choice depends on a multitude of case-specific factors
including the anatomy, underlying MR physics and sampling strategy, so hand-
tuning generally yields sub-optimal reconstruction. For reliably high quality
CS-MRI reconstruction, a robust, accurate and computationally efficient method
for automatic parameter tuning is desirable.

Many existing methods for automatic parameter tuning depend on a testing
procedure, where the reconstruction algorithm is run multiple times [13, 14,
15, 16]. In general this is a highly computationally intensive process,
especially for models with multiple parameters. The focus of this paper is on
algorithms that tunes model parameters on-the-fly, so that the algorithm only
needs to be run once [17, 18].

A natural candidate for automatic on-the-fly parameter tuning is Stein’s
Unbiased Risk Estimate (SURE) [19]. SURE efficiently estimates the squared
error of a denoised vector, so can be used as a proxy for the actual error to
tune parameters of a denoising function [20]. In CS-MRI, random sampling
ensures that the aliasing is “noise-like”, suggesting that the application of
SURE may be valid. However, SURE only gives near-optimal parameter selection
when the noise is drawn from a zero-mean Gaussian distribution, and the per-
iteration aliasing of standard CS-MRI algorithms deviates considerably from a
zero-mean Gaussian in practice. Despite no guarantee that the tuning is near-
optimal, there has been empirical success employing SURE for tuning-free CS-
MRI [21, 17], suggesting that SURE is somewhat robust to deviations from the
appropriate noise distribution.

This paper presents a solution to the problem of aliasing model mismatch via
the Approximate Message Passing (AMP) framework[22, 23, 24]. The AMP algorithm
was originally developed for sparse signals sampled with i.i.d. Gaussian
sensing matrices[22]. For such sensing matrices, it is known that AMP obeys a
“state evolution”, where the intermediate signal estimate is distributed
according to the reference signal corrupted by white Gaussian aliasing. State
evolution ensures that SURE’s assumed aliasing distribution matches the actual
aliasing, so automatic parameter tuning with SURE is truly near-optimal [25,
26]. However, when the sensing matrix does not satisfy its i.i.d. Gaussian
assumption, such as for CS-MRI, AMP’s state evolution may fail. The naive
application of AMP to CS-MRI does not obey state evolution and does not
perform well in practice. There has been a number of attempts to adapt AMP-
based algorithms to CS-MRI including: Location Constrained AMP [27], Block-
Matching 3D AMP-MRI [28, 29] Learned Denoising AMP [30] and, for single coil
measurements, Vector AMP (VAMP) for Image Reconstruction [24, 31] and Damped
Denoising VAMP [32]. Although these modifications were found to empirically
improve the convergence behaviour of AMP, they do not obey state evolution so
cannot be leveraged for near-optimal parameter tuning with SURE.

The authors’ recently proposed single-coil VDAMP algorithm [33, 34] is, to our
knowledge, the first algorithm for Fourier sampled images that obeys a state
evolution. Specifically, single-coil VDAMP obeys a “colored state evolution”,
where the per-iteration aliasing is accurately characterized by one real
number per wavelet subband [33]. This paper proposes the Parallel Variable
Density Approximate Message Passing (P-VDAMP) algorithm, which extends single-
coil VDAMP to multiple coils. For Bernoulli sampled k-space, where sampling is
random and independent, P-VDAMP’s intermediate aliasing is distributed
according to a zero-mean Gaussian with a covariance matrix represented with a
wavelet-domain vector with $N$ unique entries, where $N$ is the dimension of
the reference image. We refer to this as as “parallel colored state
evolution”. To our knowledge, P-VDAMP is the first algorithm for multi-coil
MRI that obeys a state evolution with accurately tracked parameters, and
therefore the first method where parameter tuning with SURE is truly near-
optimal.

##  2 Theory

Denote the $N$-dimensional data on coil $c$ as

$\displaystyle\bm{y}_{c}$
$\displaystyle=\bm{M}_{\Omega}(\bm{F}\bm{S}_{c}\bm{x}_{0}+\bm{\varepsilon}_{c}),\hskip
28.45274ptc=\\{1,2,\ldots,N_{c}\\}$ (1)

where $\bm{M}_{\Omega}$ is a diagonal undersampling mask with $1$ on the $j$th
diagonal entry if $j\in\Omega$ and $0$ otherwise, where $\Omega$ is a sampling
set with $|\Omega|=n$ for $n<N$. The matrix $\bm{F}$ is the discrete Fourier
transform and $\bm{S}_{c}$ is the diagonal sensitivity profile of coil $c$,
which we assume have been normalized so that they satisfy
$\sum_{c}\bm{S}_{c}^{H}\bm{S}_{c}=\mathds{1}_{N}$, and $\bm{x}_{0}$ is the
image of interest. The vector $\bm{\varepsilon}_{c}$ represents the
measurement noise, which is modelled as complex Gaussian distributed with zero
mean.

The mathematical results in this section assume that the entries of the
sampling set $\Omega$ are generated independently with arbitrary variable
density $\mathrm{Pr}(j\in\Omega)=p_{j}$, referred to as “Bernoulli sampling”.
In other words, the sampling points correspond to Bernoulli distributed phase
encodes with a fully sampled readout dimension in a Cartesian sampling
framework, and (1) is a 2D subproblem that follows an inverse Fourier
transform along the readout dimension.

The experimental work in this paper can be reproduced with publicly available
MATLAB code111https://github.com/charlesmillard/P-VDAMP.

###  2.1 Density compensation and error propagation

In this section, we demonstrate that a Bernoulli sampled zero-filled image
estimate has aliasing with a fully characterisable distribution. In other
words, we show that state evolution holds at the first iteration. The
mathematical tools developed in this section are preliminary to the full
P-VDAMP algorithm presented in Section 2.3.

Consider the density compensated coil combined estimate from zero-filled data:

$\hat{\bm{x}}_{0}=\sum_{c}\bm{S}_{c}^{H}\bm{F}^{H}\bm{P}^{-1}\bm{y}_{c},$ (2)

where $\bm{P}=\mathds{E}_{\Omega}\\{\bm{M}_{\Omega}\\}$ is the expected value
of the random sampling mask, so that $\bm{P}$ is diagonal with $j$th diagonal
$p_{j}$. The $\bm{P}^{-1}$ matrix ensures that $\hat{\bm{x}}_{0}$ is an
unbiased estimate of $\bm{x}_{0}$:

$\displaystyle\mathds{E}\\{\hat{\bm{x}}_{0}\\}$
$\displaystyle=\mathds{E}\\{\sum_{c}\bm{S}_{c}^{H}\bm{F}^{H}\bm{P}^{-1}\bm{M}_{\Omega}(\bm{F}\bm{S}_{c}\bm{x}_{0}+\bm{\varepsilon}_{c})\\}$
$\displaystyle=\sum_{c}\bm{S}_{c}^{H}\bm{F}^{H}\bm{P}^{-1}\bm{P}\bm{F}\bm{S}_{c}\bm{x}_{0}=\bm{x}_{0},$

where $\mathds{E}\\{\bm{M}_{\Omega}\\}=\bm{P}$ and
$\mathds{E}\\{\bm{\varepsilon}_{c}\\}=\bm{0}$ has been used.

It is often stated that the role of random sampling in CS-MRI is to ensure
that the undersampling artifacts of the zero-filled estimate are “noise-like”.
A primary result of this paper is to show that when $\bm{M}_{\Omega}$ is
Bernoulli random the undersampling artifacts of $\hat{\bm{x}}_{0}$ behave
exactly as Gaussian noise:
$\hat{\bm{x}}_{0}=\bm{x}_{0}+\mathcal{CN}(\bm{0},\bm{\Sigma}^{2}_{0})$, where
$\mathcal{CN}(\bm{0},\bm{\Sigma}^{2})$ denotes the complex Gaussian
distribution with zero mean and independent real and imaginary parts with
covariance matrix $\bm{\Sigma}^{2}/2$. Moreover, we show that the covariance
matrix $\bm{\Sigma}^{2}_{0}$ can be accurately and efficiently estimated with
an $N$-dimensional vector, which we denote as $\bm{\tau}_{0}$.

Image-domain aliasing is colored by the non-uniform spectral density of the
image and the variable density sampling scheme [35], so has a covariance
matrix with strong off-diagonals. We employ an orthogonal wavelet transform to
approximately decorrelate the aliasing, so that $\tilde{\bm{\Sigma}}^{2}_{0}$
is modelled as diagonal in the wavelet domain [33]. Wavelets are commonly used
in CS-MRI for sparsely representing an image, however, we emphasize that
wavelets are used here as a tool for modelling the aliasing, and that P-VDAMP
is not constrained to wavelet-domain regularisation. The competitive
performance of more sophisticated image models, including via neural networks,
was demonstrated for single-coil VDAMP in [36], although is not the focus of
this paper.

Figure 1:  The density compensated, zero-filled estimate $\hat{\bm{x}}_{0}$
for a $R=10$ Bernoulli undersampled brain dataset. The image labelled
“Aliasing model” shows $\sqrt{\bm{\tau}_{0}}$ computed via (3), and the image
labelled “True aliasing/model” shows
$\bm{\Psi}(\bm{x}_{0}-\hat{\bm{x}}_{0})\oslash\sqrt{\bm{\tau}_{0}}$, where
$\oslash$ is entry-wise division.

The appendix derives an approximation of the wavelet-domain covariance
$\bm{\tau}_{0}\approx\mathrm{diag}(\bm{\Psi}\tilde{\bm{\Sigma}}^{2}_{0}\bm{\Psi}^{H})$,
where $\bm{\Psi}$ is an orthogonal wavelet transform. The $j$th entry is

$\tau_{0,j}=\bm{\xi}_{j}^{H}\left(\sum_{i}|\hat{\Psi}_{ji}|^{2}\left[\left(\frac{1-p_{i}}{p_{i}}\right)\bm{y}_{i}\bm{y}_{i}^{H}+\bm{\Sigma}_{\varepsilon,i}^{2}\right]\frac{m_{i}}{p_{i}}\right){\bm{\xi}_{j}}.$
(3)

Here, $\hat{\bm{\Psi}}=\bm{\Psi F}^{H}$, $p_{i}$ and $m_{i}$ are the $i$th
diagonal of $\bm{P}$ and $\bm{M}_{\Omega}$ respectively, and
$\bm{\Sigma}_{\varepsilon,i}^{2}=\mathds{E}\\{\bm{\varepsilon}_{i}\bm{\varepsilon}_{i}^{H}\\}$.
The vectors $\bm{\xi}_{j},\bm{y}_{i}$ and $\bm{\varepsilon}_{i}$ have
dimensionality $N_{c}$; for instance, $\bm{y}_{i}$ is the vector of k-space
location $i$ across all coils. The notation $|\cdot|^{2}$ refers to the entry-
wise magnitude squared. As described in the appendix, the vector
$\bm{\xi}_{j}$ approximates the coil-weighted spectral densities
$\hat{\bm{\Psi}}^{c}=\bm{\Psi S}_{c}^{H}\bm{F}^{H}$ with
$\hat{\bm{\Psi}}_{ji}^{c}\approx\xi_{c,j}\hat{\bm{\Psi}}_{ji}$. We suggest
using $\xi_{c,j}=\braket{}{\bm{\Psi}_{j}}{{}^{2},\mathrm{Diag}(\bm{S}_{c})}$,
where $\bm{\Psi}_{j}$ is the $j$th row of $\bm{\Psi}$.

A visualisation of the error propagation formula is shown in Fig. 1 for a
$256\times 256$ brain dataset. Here, $N_{c}=16$ coils sensitivities were
estimated with ESPIRiT [4] via a $24\times 24$ central autocalibration region
and compressed to $8$ virtual coils with Principle Component Analysis (PCA)
[37], and an order 4 Daubechies wavelet $\bm{\Psi}$ at $4$ decomposition
scales was employed. The histograms verify that the aliasing of
$\hat{\bm{x}}_{0}$ is accurately approximated by $\bm{\tau}_{0}$, so that
$\hat{\bm{x}}_{0}\approx\bm{x}_{0}+\mathcal{CN}(\bm{0},\bm{\Psi}\mathrm{Diag}(\bm{\tau}_{0})\bm{\Psi}^{H})$.

The P-VDAMP algorithm described in Section 2.3 has such aliasing not only for
the zero-filled initial estimate $\hat{\bm{x}}_{0}$ but also for all further
iterations $\hat{\bm{x}}_{k}$:

$\hat{\bm{x}}_{k}\approx\bm{x}_{0}+\mathcal{CN}(\bm{0},\bm{\Sigma}_{k}^{2}),$
(4)

where
$\bm{\Sigma}_{k}^{2}\approx\bm{\Psi}\mathrm{Diag}(\bm{\tau}_{k})\bm{\Psi}^{H}$
and $\bm{\tau}_{k}$ is estimated with an error propagation formula based on
(3). We refer to (4) as P-VDAMP’s “parallel coloured state evolution”.

###  2.2 Automatic parameter tuning with SURE

This section details how P-VDAMP’s state evolution can be leveraged by
applying SURE [19], building on previous work on tuning-free AMP in [25, 26,
38]. Consider a random vector
$\bm{r}=\bm{w}_{0}+\mathcal{CN}(\bm{0},\mathrm{Diag}(\bm{\tau}))$, where
$\bm{\tau}$ is known, and let an estimate of $\bm{w}_{0}$ be given by
$\bm{f}(\bm{r};\bm{\Theta})$, where $\bm{\Theta}$ is a parameter vector.
Complex SURE takes the form [19, 33]

$cSURE(\bm{f}(\bm{r};\bm{\Theta});\bm{\tau})\\\
=\|\bm{f}(\bm{r};\bm{\Theta})-\bm{r}\|^{2}_{2}+\bm{\tau}^{T}(2\bm{\partial}(\bm{f}(\bm{r};\bm{\Theta})))-\bm{1}),$
(5)

where

$\partial_{i}(\bm{f}(\bm{r};\bm{\Theta}))=\frac{1}{2}\left(\frac{\partial\Re[f_{i}(\bm{r};\bm{\Theta})]}{\partial\Re[r_{i}]}+\frac{\partial\Im[f_{i}(\bm{r};\bm{\Theta})]}{\partial\Im[r_{i}]}\right).$
(6)

Here, $\Re[\cdot]$ and $\Im[\cdot]$ are the real and imaginary parts
respectively. cSURE is an unbiased estimate of the squared error of the
denoised image,

$\displaystyle\mathds{E}\\{\|\bm{f}(\bm{r};\bm{\Theta})-\bm{w}_{0}\|^{2}_{2}\\}=\mathds{E}\\{cSURE(\bm{f}(\bm{r};\bm{\Theta});\bm{\tau})\\},$
(7)

where $\mathds{E}\\{\cdot\\}$ is the expectation over the random Gaussian
noise. Therefore cSURE can be minimized as a proxy for minimizing the true
squared error without requiring the ground truth $\bm{w}_{0}$, so a near-
optimal denoiser satisfying

$\min_{\bm{\Theta}}cSURE(\bm{f}(\bm{r};\bm{\Theta});\bm{\tau})\approx\min_{\bm{\Theta}}\|\bm{f}(\bm{r};\bm{\Theta})-\bm{w}_{0}\|^{2}_{2}$
(8)

can be employed.

In this paper, as in [33], we use the aliasing model (3) in conjunction with
cSURE to employ soft thresholding with a per-subband threshold. To tune such a
model by hand would be highly impractical; for instance, per-subband
thresholding of a wavelet transform with $4$ decomposition scales would
require 13 parameters to be tuned by hand.

###  2.3 Statement of algorithm

Algorithm 1 P-VDAMP

Require: Sampling set $\Omega$, coil sensitivities $\bm{S}_{c}$, wavelet
transform $\bm{\Psi}$, probability matrix $\bm{P}$, measurements $\bm{y}_{c}$,
measurement noise covariance $\bm{\Sigma}_{\varepsilon}^{2}$, denoiser
$\bm{g}(\bm{r};\bm{\tau})$, number of iterations $K_{it}$.

1: Compute $|\hat{\bm{\Psi}}|^{2}=|\bm{\Psi}\bm{F}^{H}|^{2}$ and
$\xi_{c,j}=\braket{}{\bm{\Psi}_{j}}{{}^{2},\mathrm{Diag}(\bm{S}_{c}^{H})}$ and
set $\widetilde{\bm{r}}_{0}=\bm{0}_{N}$

2: for $k=0,1,\ldots,K_{it}-1$ do

3:
$\bm{z}_{k,c}=\bm{y}_{c}-\bm{M}_{\Omega}\bm{FS}_{c}\bm{\Psi}^{H}\widetilde{\bm{r}}_{k}$

4:
$\bm{r}_{k}=\widetilde{\bm{r}}_{k}+\bm{\Psi}\sum_{c}\bm{S}_{c}^{H}\bm{F}^{H}\bm{P}^{-1}\bm{z}_{k,c}$

5:
$\tau_{k,j}=\bm{\xi}_{j}^{H}\left(\sum_{i}|\hat{\Psi}_{ji}|^{2}\left[\left(\frac{1-p_{i}}{p_{i}}\right)\bm{z}_{k,i}\bm{z}_{k,i}^{H}+\bm{\Sigma}_{\varepsilon,i}^{2}\right]\frac{m_{i}}{p_{i}}\right)\bm{\xi}_{j}$

6: $\hat{\bm{w}}_{k}=\bm{g}(\bm{r}_{k};\bm{\tau}_{k})$

7:
$\bm{\alpha}_{k}=\braket{\bm{\partial}(\bm{g}(\bm{r}_{k};\bm{\tau}_{k}))}_{\mathrm{sband}}$

8: $\bm{c}_{k}=\bm{1}_{N}\oslash(\bm{1}_{N}-\bm{\alpha}_{k})$

9:
$\widetilde{\bm{r}}_{k+1}=\bm{c}_{k}\odot(\hat{\bm{w}}_{k}-\bm{\alpha}_{k}\odot\bm{r}_{k})$

10: end for

11: return
$\hat{\bm{x}}=\bm{\Psi}^{H}\hat{\bm{w}}_{k}+\sum_{c}\bm{S}_{c}^{H}\bm{F}^{H}(\bm{y}_{c}-\bm{M}_{\Omega}\bm{FS}_{c}\bm{\Psi}^{H}\hat{\bm{w}}_{k})$
or $\hat{\bm{x}}=\bm{\Psi}^{H}\bm{r}_{k}$

The P-VDAMP algorithm is shown in Algorithm 1. Here, vectors with subscript
$c$ ($\bm{y}_{c}$ and $\bm{z}_{k,c}$) refer to an $N$-dimensional vector on
the $c$th coil, and vectors with subscript $i$ or $j$ ($\bm{\xi}_{j}$ and
$\bm{z}_{k,i}$) refer to the $N_{c}$-dimensional vector at the $i$th or $j$th
k-space coefficient. In lines 3 and 5, the computation is for all coils $c$
and k-space locations $j$ respectively.

P-VDAMP is closely related to the single-coil VDAMP algorithm [33]. The
crucial difference between VDAMP and P-VDAMP is the aliasing model update in
line 5. While VDAMP represents the aliasing with one real number per wavelet
subband, P-VDAMP uses (3) to represent the aliasing with an $N$-dimensional
vector, modelling the intra-subband modulation induced by the coil sensitivity
profiles.

Lines 3-4, referred to as “density compensated gradient descent”, is the
multi-coil version of lines 3-4 of VDAMP [33]. This can be understood as
gradient descent with a modified forward model, so that (1) is left-multiplied
by $\bm{P}^{-1/2}$,

$\displaystyle\bm{P}^{-1/2}\bm{y}_{c}$
$\displaystyle=\bm{P}^{-1/2}\bm{M}_{\Omega}(\bm{F}\bm{S}_{c}\bm{x}_{0}+\bm{\varepsilon}_{c}),$
(9)

for all coils. This weighting ensures that the columns of the multi-coil
sensing matrix are normalized in expectation over the random mask, a necessary
component of state evolution. At the first iteration $k=0$, where
$\widetilde{\bm{r}}_{0}=\bm{0}$, lines 3-4 correspond to the density-
compensated zero-filled estimate stated in (2). Using density compensated
gradient descent implies that the data consistency term in P-VDAMP’s cost
function is a weighted $\ell_{2}$ norm,
$\sum_{c}\|\bm{P}^{-1/2}(\bm{y}_{c}-\bm{M}_{\Omega}\bm{F}\bm{S}_{c}\bm{x})\|_{2}^{2}$,
the consequences of which are discussed in Section 4.1.

Line 5 applies the aliasing model update stated in (3) to the vector
$\bm{z}_{k}$, which estimates the wavelet-domain covariance matrix of the
aliasing of $\bm{r}_{k}$. Lines 6-9 are identical to lines 6-9 of the single-
coil VDAMP algorithm [33]. In line 6, a function
$\bm{g}(\bm{r}_{k};\bm{\tau}_{k})$ is applied to the wavelet-domain estimate
$\bm{r}_{k}$. Since, by state evolution, $\bm{r}_{k}$ is the Gaussian
corrupted ground truth, $\bm{g}(\bm{r}_{k};\bm{\tau}_{k})$ is referred to as a
denoiser. In this work, the denoiser was soft thresholding with per-subband
soft thresholds tuned with cSURE via the aliasing model $\bm{\tau}_{k}$, as
described in Section 2.2. Lines 7-9 correspond to a correction of the denoiser
referred to in the AMP literature as “Onsager correction”, which removes the
correlation that causes state evolution to fail for standard CS-MRI algorithms
[22, 23]. In line 7, the notation $\braket{\cdot}_{\mathrm{sband}}$ averages
subbands, so that $\bm{\alpha}_{k}$ has the structure

$\bm{\alpha}_{k}=\begin{bmatrix}\alpha_{k,1}\bm{1}_{N_{1}}\\\
\alpha_{k,2}\bm{1}_{N_{2}}\\\ \vdots\\\
\alpha_{k,1+3s}\bm{1}_{N_{1+3s}}\end{bmatrix},$ (10)

with

$\alpha_{k,b}=\frac{1}{N_{b}}\sum_{j\in
J_{b}}\partial_{j}(\bm{g}(\bm{r}_{k};\bm{\tau}_{k})),$

where $J_{b}$ is the set of indices associated with subband $b$,
$N_{b}=|J_{b}|$ and $\bm{\partial}(\cdot)$ is defined in (6). The notation
$\oslash$ in line 8 refers to entry-wise division and $\odot$ in line 9 refers
to entry-wise multiplication.

Line 11 of P-VDAMP gives two options for the algorithm output. The first uses
gradient descent without density compensation. We also suggest
$\bm{\Psi}^{H}\bm{r}_{k}$, which is unbiased when state evolution holds. The
qualitative and quanitative differences between these outputs are discussed in
Section 4.

###  2.4 Comments on P-VDAMP in practice

We have found that some damping aids P-VDAMP’s convergence, so that the
denoiser is damped by a constant factor $0<\rho\leq 1$, as in [39, 24]. For
$k>0$, we suggest replacing lines 6 and 7 of P-VDAMP, Algorithm 1, with

$\displaystyle\hat{\bm{w}}_{k}$
$\displaystyle=\rho\bm{g}(\bm{r}_{k};\bm{\tau}_{k})+(1-\rho)\bm{\hat{w}}_{k-1}$
$\displaystyle\bm{\alpha}_{k}$
$\displaystyle=\rho\braket{\bm{\partial}(\bm{g}(\bm{r}_{k};\bm{\tau}_{k}))}_{\text{sband}}.$

Note that $\rho=1$ corresponds to zero damping. In all experiments in this
paper, the damping factor was set to $\rho=0.75$.

Rather than using a fixed number of iterations, we used a more practical
stopping criterion based on the aliasing model $\bm{\tau}_{k}$. [36] In the
experiments in this paper, P-VDAMP stops iterating when the mean-squared of
the aliasing is predicted to have increased,
$\braket{\bm{\tau}_{k}}>\braket{\bm{\tau}_{k-1}}$, or when
$\braket{\bm{\tau}_{k}}$ changes by a small amount:
$|\braket{\bm{\tau}_{k}}-\braket{\bm{\tau}_{k-1}}|/\braket{\bm{\tau}_{k-1}}<\epsilon$.
We used $\epsilon=10^{-3}$ throughout this paper.

##  3 Methods

###  3.1 Description of data

P-VDAMP was evaluated on three types of MRI data:

  * •

Brain. $256\times 256$ acquisitions acquired on $N_{c}=16$ coils, from an on-
site Siemens 3T scanner.

  * •

Knee. $368\times 640$ acquisitions acquired on $N_{c}=15$ coils, from the
publicly available fastMRI dataset [40], which includes a mixture of field
strengths.

  * •

Time-of-flight Angiography. $640\times 506\times 56$ acquisitions acquired on
$N_{c}=32$ coils, from an on-site Siemens 7T scanner, decomposed in to a
series of 2D problems by taking an inverse Fourier transform along the first
dimension.

For each of the three categories, we used two fully sampled datasets. The two
knee datasets were randomly selected from the fastMRI dataset.

###  3.2 Experimental method

All datasets were retrospectively undersampled at factors $R=5,10$. with
variable density Bernoulli masks. The masks had $p_{j}=1$ in a $24\times 24$
central autocalibration region used to estimate the coil sensitivities
$\bm{S}_{c}$ with ESPIRiT [4] via the BART toolbox [41]. An example of the
sampling probability and instance of a corresponding mask for a $256\times
256$ k-space is shown in Fig. 2.

To reduce the reconstruction time, the coils were compressed to $8$ virtual
coils with PCA [37]: equation (3) is linear in $N$ but quadratic in $N_{c}$,
so P-VDAMP especially benefits from coil compression. For all algorithms we
used the order 4 Daubechies mother wavelet at $4$ decompositions scales for
$\bm{\Psi}$. P-VDAMP’s sensitivity to measurement noise is evaluated in
Section 4.3; otherwise, the covariance $\bm{\Sigma}_{\varepsilon}^{2}$ was set
to a matrix of zeros.

(a) Sampling density $\bm{P}$

(b) Random mask $\bm{M}_{\Omega}$

Figure 2: The sampling density and Bernoulli random mask employed in this
paper for a $256\times 256$ k-space at $R=10$.

###  3.3 Comparative algorithms

We evaluated the reconstruction quality of both output options in line 11 of
P-VDAMP, Algorithm 1, referring to the $\bm{\Psi}^{H}\bm{r}_{k}$ output as
“Unbiased P-VDAMP”, and the alternative as simply “P-VDAMP”.

For comparative hand-tuned methods, we implemented the Fast Iterative
Shrinkage-Thresholding algorithm (FISTA) [42] with a $\lambda$ tuned using the
reference to be approximately MSE optimal, referred to as “Optimal FISTA”,
using an exhaustive search within a reasonable range. Secondly, to reflect the
more realistic case where the sparse parameter is not optimally tuned, we also
ran FISTA with a sparse weighting shared across all reconstructions, which we
refer to as “Shared FISTA”. For the shared sparse weighting we used the mean
of the optimally tuned sparse weightings.

For comparative tuning-free methods, we ran FISTA with a per-subband threshold
tuned automatically with SURE using a white alising model, referred to as
“SURE-IT”[17]. We also evaluated a recently proposed parameter-free variant of
FISTA, which automatically selects a per-scale sparse weighting based on
k-means clustering of the zero-filled estimate, referred to in this work as
“Automatic-FISTA” (A-FISTA)[18].

The stopping criteria employed for Optimal FISTA, Shared FISTA and SURE-IT was
analogous to the P-VDAMP criteria described Section 2.4. For A-FISTA, we found
that such a stopping criterion typically stopped the algorithm prematurely, so
to give A-FISTA the best chance of a high-quality reconstruction we used a
fixed, large number of iterations. As in [18], we chose 200 iterations.

Figure 3: Reconstruction quality for all algorithms, where the result for each
data type has been averaged over the two datasets.

###  3.4 Performance metrics

Prior to computing performance metrics, we masked the images to mitigate
background effects. For the brain and knee datasets, all pixels with absolute
value less than 5% of the maximum of the reference were masked. We also skull-
stripped the angiograms with a hand-selected mask [43] and computed the
maximum projection along the readout dimension.

To quantify reconstruction quality we used the Normalized Mean-Squared Error
(NMSE), Structural Similarity Index Metric (SSIM) [44] and High-Frequency
Error Norm (HFEN) [45], which quantifies the reconstruction quality of the
image’s fine features and is computed by applying a $15\times 15$ Laplacian of
Gaussian filter with standard deviation $1.5$ to $\hat{\bm{x}}-\bm{x}_{0}$ and
computing the $\ell_{2}$ norm. We normalized the HFEN by dividing by the
$\ell_{2}$ norm of the reference image. For the angiograms, we computed the
SSIM after applying a vessel mask computed by thresholding with Otsu’s method
[46, 47].

Figure 4: Reconstructions of a brain dataset sampled at $R=5$ with error maps.
Shared FISTA’s sparse weighting is around a factor of 10 larger than for
Optimal FISTA, which led to significant reconstruction artifacts. P-VDAMP
performs competitively with Optimal FISTA and substantially better than the
competing parameter-free methods SURE-IT and A-FISTA. The scale on the
colorbar indicates that its maximum is 0.08 times the maximum of the reference
image.

Figure 5: Reconstructions of a knee dataset sampled at $R=10$, with error
maps, cropped to a central $320\times 320$ region. Unbiased P-VDAMP’s NMSE is
2.5dB higher than the next worst, SURE-IT, but has competitive HFEN and
performs well qualitatively. The competing parameter-free algorithms SURE-IT
and A-FISTA smooth over details fine details in the bone.

Figure 6: An angiogram sampled at $R=5$, with a magnified $100\times 100$
region. The green arrow shows a challenging narrow blood vessel that is
retained well for Unbiased P-VDAMP but mostly lost in competing methods.
Figure 7: NMSE vs time for a brain and knee dataset at $R=5$ and $R=10$
respectively.

##  4 Results and Discussion

###  4.1 Reconstruction quality

The NMSE, HFEN and SSIM of each is algorithm is shown in Fig. 3. P-VDAMP
consistently outperformed the competing parameter-free methods SURE-IT and
A-FISTA. The improvement over SURE-IT emphasizes that P-VDAMP’s state
evolution is key for high-quality threshold selection with SURE. A-FISTA’s
performace is quite varied, which indicates that the quality of its parameter
tuning is not as robust as P-VDAMP. The highly mixed performance of Shared
FISTA shows that that the choice of sparse weighting is key, and that sharing
its value for different types of data cannot be expected to yield a
consistently high quality reconstruction. P-VDAMP peformed competitively with
Optimal FISTA for all reconstructions, indicating that it is robust to the
forward model and anatomy. Since the fully sampled reference is not available
for parameter tuning in prospectively undersampled scans, FISTA’s parameter
tuning cannot be expected to be near-optimal, so P-VDAMP is likely to
outperform FISTA in practice. Unbiased P-VDAMP’s NMSE and SSIM was generally
poorer than P-VDAMP, but had comparable HFEN.

Since P-VDAMP’s sparse model has one parameter per subband while Optimal FISTA
has a single, global parameter, one might expect that P-VDAMP would
consistently outperform Optimal FISTA. There are two reasons why P-VDAMP does
not always outperform Optimal FISTA. Firstly, P-VDAMP’s $\bm{P}^{-1}$ in lines
3-4 of Algorithm 1, which is essential for state evolution, is not
advantageous from a reconstruction error perpective as it boosts errors at
high frequencies, where the sampling probability is low. P-VDAMP outperforms
Optimal FISTA when the benefits of its richer sparse model “wins” over the
disadvantage of needing to include $\bm{P}^{-1}$. The second reason is that
damping causes P-VDAMP’s aliasing model $\bm{\tau}_{k}$ to be imperfect,
especially for high $k$, as discussed in Section 4.2, so parameter selection
with SURE is not always near-optimal in practice.

Figure 8:  P-VDAMP’s state evolution for the same brain dataset and sampling
pattern as shown in Fig. 1 at iterations $k=1,2,3$, showing that the Gaussian
aliasing illustrated at $k=0$ approximately holds for subsequent iterations.

Figure 9:  P-VDAMP’s state evolution for a knee dataset sampled at $R=10$ at
iterations $k=0,5,10$.

Fig. 6 shows one of the brain datasets at $R=5$. Optimal FISTA and both
P-VDAMP estimates retain considerably more of the fine detail than competing
methods. Shared FISTA’s $\lambda$ is around $10$ larger than Optimal FISTA in
this case, and compression artifacts are clearly visible. The competing
parameter-free methods SURE-IT and A-FISTA peform better than Shared FISTA in
this case, but still have substantial reconstruction artifacts. The underlying
anatomy is not clearly visible in Unbiased P-VDAMP’s error maps, with a
similar contribution to the error at tissue boundaries and bulk tissue, which
is a consequence of P-VDAMP’s state evolution.

One of the knee datasets undersampled at $R=10$ is shown in Fig. 6. P-VDAMP
outperforms both competing parameter-free methods. Shared FISTA’s $\lambda$ is
6 times smaller than Optimal FISTA here, and performs reasonably well, but has
a higher NMSE than P-VDAMP.

Fig. 6 shows the reconstructions of one of the Angiograms at $R=5$, where
P-VDAMP outperforms Optimal FISTA. The green arrow in the magnified region
highlights a narrow region of a blood vessel, which is visible for Unbiased
P-VDAMP and partially visible for P-VDAMP but mostly lost in the competing
methods.

The $\bm{\tau}_{k}$ update in line 5 of P-VDAMP, Algorithm 1, is reasonably
computationally expensive; P-VDAMP took around 4 times longer per iteration
than FISTA in our MATLAB implementation. However, we found that P-VDAMP
typically requires substantially fewer iterations, so required only a little
longer than Optimal FISTA to reach convergence. This is illustrated in two
NMSE vs time plots shown in Fig. 7.

###  4.2 P-VDAMP’s state evolution

Fig. 9 shows the wavelet-domain aliasing for one of the brain datasets at
$k=1,2,3$ and one of the knee datasets at iterations $k=0,5,10$, both at
$R=10$, which shows that P-VDAMP’s parallel coloured state evolution (4)
approximately holds. There is a slight deviation from exactly unbiased
aliasing apparent at larger $k$. This is due to P-VDAMP’s damping factor
$\rho$, which causes P-VDAMP’s Onsager correction to be inexact. We found that
using less damping caused state evolution to be more exact, but that the
algorithm converged to a poorer reconstruction. We therefore recommend using
damping at the cost of a slight breakdown in state evolution.

###  4.3 Robustness to measurement noise

To evaluate the effect of measurement noise on P-VDAMP’s performance, we
simulated Gaussian additive noise. We constructed an artificial measurement
noise covariance matrix by generating a $N_{c}$-dimensional random vector
$\bm{v}$ with complex entries with real and imaginary parts drawn from a
uniform distribution between $-1$ and $1$, and scaled the vector to
$|\bm{v}|=c\bm{1}_{N_{c}}$, where $c$ is a constant that determines the
signal-to-noise ratio (SNR). For each measurement noise SNR we averaged over
10 random instances of $\bm{v}$. We compared P-VDAMP with FISTA with a fixed
sparse weighting, tuned to be optimal in the absence of simulated measurement
noise.

Fig. 10 shows the dependence of the reconstruction NMSE on measurement noise
SNR for each algorithm for one of the knee datasets sampled at $R=5,10$. While
FISTA’s sparse weighting is fixed, P-VDAMP’s sparse parameters adapt to the
measurement noise on-the-fly and converges to a better solution for low SNR in
these examples. Since the measurement noise is amplified by $\bm{P}^{-1}$ in
P-VDAMP, the unbiased P-VDAMP output is more sensitive to measurement noise,
and did not perform well for low SNR, indicating that it may not be
appropriate for especially low SNR CS-MRI applications.

##  5 Conclusions

P-VDAMP is an efficient, robust and principled tuning-free algorithm for CS-
MRI that is competitive with Optimal FISTA and offers substantial robustness
and reconstruction quality improvements over competing parameter-free methods.

State evolution essentially recasts multi-coil CS-MRI reconstruction as a
sequence of Gaussian denoising tasks, where the Gaussian “effective noise” has
known distribution parameters. In other words, roughly speaking, how well one
can reconstruct an image depends on how well one can denoise a Gaussian-
corrupted version of it [29]. This suggests a number of potential advantages
of P-VDAMP’s state evolution beyond parameter tuning that could be explored in
future work. For instance, benefits for deep-learning based reconstruction
[48, 49, 36] and uncertainty quantification [50] have already been explored in
the context of single-coil VDAMP [36]. Initial work by the authors on
Denoising-P-VDAMP (D-P-VDAMP) [51, 52] has shown that the deep Plug-and-Play
framework [53, 54] can improve P-VDAMP’s performance.

We have found that P-VDAMP naively applied to Poisson disc sampled data also
performs well, despite a violation of the independent sampling assumption.
However, we have found that it typically diverges when applied to problems
that contain smooth sampling trajectories, so cannot yet be applied to
sampling schemes such as 3D radial [55] or 3D cones [56], nor to purely 2D
acquisitions. P-VDAMP could be extended to such sampling schemes by an
appropriate modification of the aliasing model update in line 5 of Algorithm
1.

Very recently an AMP-based algorithm was proposed [57] for single and multi-
coil MRI acquisitions that uses similar principles to single-coil VDAMP but
with an alternative procedure for data consistency at each iteration, which
may have better fixed points. Extending the richer aliasing model presented in
this paper to the multi-coil version of this alternative algorithm is a
desirable avenue for future work.

This appendix derives an approximate expression for the diagonal of the
covariance of the aliasing of $\bm{r}_{0}=\bm{\Psi}\hat{\bm{x}}_{0}$, where
$\hat{\bm{x}}_{0}=\sum_{c}\bm{S}_{c}^{H}\bm{F}^{H}\bm{P}^{-1}\bm{y}_{c}$. Let
$\bm{y}_{0,c}=\bm{F}\bm{S}_{c}\bm{x}_{0}$ and $\bm{w}_{0}=\bm{\Psi}\bm{x}_{0}$
be the signal of interest’s coil-weighted Fourier coefficients and wavelet
coefficients respectively. Using
$\bm{x}_{0}=\sum_{c}\bm{S}_{c}^{H}\bm{F}^{H}\bm{y}_{0,c}$, the aliasing of
$\bm{r}_{0}$ is

$\displaystyle|{\bm{r}}_{0}-\bm{w}_{0}|^{2}=$
$\displaystyle|\bm{\Psi}(\hat{\bm{x}}_{0}-\bm{x}_{0})|^{2}$ $\displaystyle=$
$\displaystyle|\bm{\Psi}\sum_{c}\bm{S}_{c}^{H}\bm{F}^{H}[\bm{P}^{-1}\bm{y}_{c}-\bm{y}_{0,c}]|^{2}$
$\displaystyle=$
$\displaystyle|\sum_{c}\hat{\bm{\Psi}}^{c}[\bm{P}^{-1}\bm{y}_{c}-\bm{y}_{0,c}]|^{2},$

where $\hat{\bm{\Psi}}^{c}=\bm{\Psi}\bm{S}_{c}^{H}\bm{F}^{H}$. Since
$\bm{y}_{c}=\bm{M}_{\Omega}(\bm{y}_{0,c}+\bm{\varepsilon}_{c})$,

$\displaystyle|\sum_{c}\hat{\bm{\Psi}}^{c}[\bm{P}^{-1}\bm{y}_{c}-\bm{y}_{0,c}]|^{2}$
$\displaystyle=|\sum_{c}\hat{\bm{\Psi}}^{c}[(\bm{P}^{-1}\bm{M}_{\Omega}-\mathds{1})\bm{y}_{0,c}-\bm{P}^{-1}\bm{M}_{\Omega}\bm{\varepsilon}_{c}]|^{2}$
$\displaystyle=|\sum_{c}\hat{\bm{\Psi}}^{c}\bm{q}_{c}|^{2}$ (11)

where we define
$\bm{q}_{c}=(\bm{P}^{-1}\bm{M}_{\Omega}-\mathds{1})\bm{y}_{0,c}-\bm{P}^{-1}\bm{M}_{\Omega}\bm{\varepsilon}_{c}.$
Computing the expectation of (11) with a method analogous to the single-coil
approach in [33] would involve computing and storing in memory the $N\times N$
matrices $|\bm{\Psi}^{c}|^{2}$ for all $N_{c}$ coils, which is not practical
for large $N$. We suggest approximating the $\hat{\bm{\Psi}}^{c}$ in terms of
the unweighted spectra $\bm{\hat{\Psi}}=\bm{\Psi F}^{H}$, so that the spectral
density of the $j$th row of $\hat{\bm{\Psi}}^{c}$ is expressed as

$\hat{\bm{\Psi}}^{c}_{j}\approx\xi_{c,j}\hat{\bm{\Psi}}_{j},$ (12)

where $\xi_{c,j}\in\mathds{C}$ is a constant. Since $\hat{\bm{\Psi}}^{c}_{j}$
is the $j$th coil-weighted wavelet spectrum, (12) is interpretable as
approximating the coil sensitivity as flat over the support of the $j$th
wavelet filter, so that its spectrum is the unweighted spectrum multiplied by
a filter- and coil-dependent scalar. We suggest using the $\xi_{c,j}$ that
minimizes $\|\hat{\bm{\Psi}}_{j}^{c}-\xi\hat{\bm{\Psi}}_{j}\|^{2}_{2}$:

$\displaystyle\xi_{c,j}$
$\displaystyle=\operatorname*{arg\,min}_{\xi}\|\hat{\bm{\Psi}}_{j}^{c}-\xi\hat{\bm{\Psi}}_{j}\|^{2}_{2}$
$\displaystyle=\operatorname*{arg\,min}_{\xi}\|\bm{\Psi}_{j}\bm{S}_{c}^{H}\bm{F}^{H}-\xi\bm{\Psi}_{j}\bm{F}^{H}\|^{2}_{2}$
$\displaystyle=\frac{\bm{\Psi}_{j}\bm{S}_{c}^{H}\bm{\Psi}_{j}^{H}}{\bm{\Psi}_{j}\bm{\Psi}_{j}^{H}}$
$\displaystyle=\braket{}{\bm{\Psi}_{j}}{{}^{2},\mathrm{Diag}(\bm{S}_{c}^{H})},$

where $\bm{\Psi}_{j}$ is the $j$th row of $\bm{\Psi}$, and we have used the
orthonormality of the rows of $\bm{\Psi}_{j}$ to set
$\bm{\Psi}_{j}\bm{\Psi}_{j}^{H}=1$.

Figure 10: Reconstruction NMSE’s dependence on simulated measurement noise SNR
for one of the knee datasets.

Claim 1. Let ${\xi}_{c,j}\in\mathds{C}$ be a constant that approximates
$\hat{\bm{\Psi}}_{j}^{c}\approx\xi_{c,j}\hat{\bm{\Psi}}_{j}$, for all $j$ rows
of $\bm{\hat{\Psi}}^{c}$. Then the expectation of (11) is

$\mathds{E}\\{|r_{0,j}-w_{0,j}|^{2}\\}\\\
\approx\bm{\xi}_{j}^{H}\left(\sum_{i}|\hat{\Psi}_{ji}|^{2}\left[\left(\frac{1-p_{i}}{p_{i}}\right)\bm{y}_{0,i}\bm{y}_{0,i}^{H}+\frac{1}{p_{i}}\bm{\Sigma}_{\varepsilon,i}^{2}\right]\right)\bm{\xi}_{j},$
(13)

where the $\bm{y}_{0,i}\in\mathds{C}^{N_{c}}$ and
$\bm{\xi}_{0,j}\in\mathds{C}^{N_{c}}$ are comprised of the $i$th k-space
location over all coils, having $c$th entry
$[y_{0,i}]_{c}=[\bm{FS}_{c}\bm{x}_{0}]_{i}$ and
$[\bm{\xi}_{0,j}]_{c}=\xi_{c,j}$, and
$\bm{\Sigma}_{\varepsilon,i}^{2}=\mathds{E}\\{\bm{\varepsilon}_{i}\bm{\varepsilon}_{i}^{H}\\}$
where $\bm{\varepsilon}_{i}\in\mathds{C}^{N_{c}}$ similarly has $c$th entry
$\varepsilon_{i,c}=[\bm{\varepsilon}_{c}]_{i}$.

######  Proof .1.

Using (12), the $j$th entry of (11) is approximately

$\displaystyle|{{r}}_{0,j}-{w}_{0,j}|^{2}$
$\displaystyle\approx|\sum_{c}\hat{\bm{\Psi}}_{j}\bm{q}_{c}\xi_{c,j}|^{2}$
$\displaystyle=|\sum_{c,i}\hat{{\Psi}}_{ji}{q}_{c,i}\xi_{c,j}|^{2}.$ (14)

Using the result from Appendix C of [33], the expectation of (14) over the
sampling mask and measurement noise is

$\displaystyle\mathds{E}\\{|\sum_{c,i}\hat{{\Psi}}_{ji}{q}_{c,i}\xi_{c,j}|^{2}\\}$
$\displaystyle=\sum_{i}|\hat{{\Psi}}_{ji}|^{2}\mathds{E}\\{|\sum_{c}\xi_{c,j}q_{c,i}|^{2}\\}$
$\displaystyle=\sum_{i}|\hat{{\Psi}}_{ji}|^{2}\mathds{E}\\{|{\bm{\xi}_{j}^{H}}\bm{q}_{i}|^{2}\\},$
(15)

where $\bm{\xi}_{j}\in\mathds{C}^{N_{c}}$ and
$\bm{q}_{i}\in\mathds{C}^{N_{c}}$. Resolving the expectation,

$\displaystyle\mathds{E}$
$\displaystyle\\{|{\bm{\xi}_{j}^{H}}\bm{q}_{i}|^{2}\\}=\mathds{E}\\{{\bm{\xi}_{j}^{H}}\bm{q}_{i}\bm{q}_{i}^{H}\bm{\xi}_{j}\\}$
$\displaystyle=$
$\displaystyle\mathds{E}\left\\{\bm{\xi}_{j}^{H}\left[\left(\frac{m_{i}}{p_{i}}-1\right)^{2}\bm{y}_{0,i}\bm{y}_{0,i}^{H}+\left(\frac{m_{i}}{p_{i}}\right)^{2}\bm{\varepsilon}_{i}\bm{\varepsilon}_{i}^{H}\right]\bm{\xi}_{j}\right\\},$
(16)

where $\mathds{E}\\{\bm{\varepsilon_{i}}\\}=\bm{0}_{N_{c}}$ has been applied
to evaluate the cross-terms to zero. For Bernoulli sampling,

$\mathds{E}\left\\{\left(\frac{m_{i}}{p_{i}}-1\right)^{2}\right\\}=\frac{1-p_{i}}{p_{i}},$

Therefore, defining
$\mathds{E}\\{\bm{\varepsilon}_{i}\bm{\varepsilon}_{i}^{H}\\}=\bm{\Sigma}_{\varepsilon,i}^{2}$,
(16) is

$\displaystyle\mathds{E}\\{|\bm{\xi}_{j}^{H}\bm{q}_{i}|^{2}\\}=\bm{\xi}_{j}^{H}\left[\left(\frac{1-p_{i}}{p_{i}}\right)\bm{y}_{0,i}\bm{y}_{0,i}^{H}+\frac{1}{p_{i}}\bm{\Sigma}_{\varepsilon,i}^{2}\right]\bm{\xi}_{j}.$
(17)

Substituting (17) in to (15) gives (13) as required.

Eqn. (13) requires knowledge of $\bm{y}_{0}$ so is of limited practical use.
However, we can approximate
$\mathds{E}_{\Omega,\varepsilon}\\{|r_{0,j}-w_{0,j}|^{2}\\}$ without the
ground truth with

$\tau_{j}=\bm{\xi}_{j}^{H}\left(\sum_{i}|\hat{\Psi}_{ji}|^{2}\left[\left(\frac{1-p_{i}}{p_{i}}\right)\bm{y}_{i}\bm{y}_{i}^{H}+\bm{\Sigma}_{\varepsilon,i}^{2}\right]\frac{m_{i}}{p_{i}}\right){\bm{\xi}_{j}}.$
(18)

Claim 2. The estimate of the multi-coil wavelet-domain error stated in (18) is
unbiased:

$\mathds{E}\\{\tau_{j}\\}=\bm{\xi}_{j}^{H}\left(\sum_{i}|\hat{\Psi}_{ji}|^{2}\left[\left(\frac{1-p_{i}}{p_{i}}\right)\bm{y}_{0,i}\bm{y}_{0,i}^{H}+\frac{1}{p_{i}}\bm{\Sigma}_{\varepsilon,i}^{2}\right]\right){\bm{\xi}_{j}}.$
(19)

######  Proof .2.

Since $\bm{y}_{i}=m_{i}(\bm{y}_{0,i}+\bm{\varepsilon}_{i})$, and
$\mathds{E}\\{\bm{\varepsilon}_{i}\\}=\bm{0}_{N_{c}}$,

$\displaystyle\mathds{E}\left\\{\left[\left(\frac{1-p_{i}}{p_{i}}\right)\bm{y}_{i}\bm{y}_{i}^{H}+\bm{\Sigma}_{\varepsilon,i}^{2}\right]\frac{m_{i}}{p_{i}}\right\\}$
$\displaystyle=\mathds{E}\left\\{\left[m_{i}\left(\frac{1-p_{i}}{p_{i}}\right)(\bm{y}_{0,i}+\bm{\varepsilon}_{i})(\bm{y}_{0,i}+\bm{\varepsilon}_{i})^{H}+\bm{\Sigma}_{\varepsilon,i}^{2}\right]\frac{m_{i}}{p_{i}}\right\\}$
$\displaystyle=\left(\frac{1-p_{i}}{p_{i}}\right)\bm{y}_{0,i}\bm{y}_{0,i}^{H}+\frac{1}{p_{i}}\bm{\Sigma}_{\varepsilon,i}^{2}$
(20)

Substituting (20) in to the expectation of (18) gives the right-hand-side of
(19).

## References

  * [1] J. B. Ra and C. Y. Rim, “Fast imaging using subencoding data sets from multiple detectors,” Magnetic resonance in medicine, vol. 30, no. 1, pp. 142–145, 1993.
  * [2] K. P. Pruessmann, M. Weiger, M. B. Scheidegger, and P. Boesiger, “SENSE: sensitivity encoding for fast MRI.,” Magnetic resonance in medicine, vol. 42, pp. 952–62, nov 1999.
  * [3] M. A. Griswold, P. M. Jakob, R. M. Heidemann, M. Nittka, V. Jellus, J. Wang, B. Kiefer, and A. Haase, “Generalized autocalibrating partially parallel acquisitions (GRAPPA),” Magnetic Resonance in Medicine, vol. 47, pp. 1202–1210, jun 2002.
  * [4] M. Uecker, P. Lai, M. J. Murphy, P. Virtue, M. Elad, J. M. Pauly, S. S. Vasanawala, and M. Lustig, “ESPIRiT-an eigenvalue approach to autocalibrating parallel MRI: Where SENSE meets GRAPPA,” Magnetic Resonance in Medicine, vol. 71, pp. 990–1001, mar 2014.
  * [5] E. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Transactions on Information Theory, vol. 52, pp. 489–509, feb 2006.
  * [6] D. L. Donoho, “Compressed sensing,” IEEE Transactions on Information Theory, vol. 52, pp. 1289–1306, apr 2006.
  * [7] M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The application of compressed sensing for rapid MR imaging,” Magnetic Resonance in Medicine, vol. 58, pp. 1182–1195, dec 2007.
  * [8] D. Liang, B. Liu, J. Wang, and L. Ying, “Accelerating SENSE using compressed sensing,” Magnetic Resonance in Medicine, vol. 62, pp. 1574–1584, dec 2009.
  * [9] A. Hsiao, M. Lustig, M. T. Alley, M. Murphy, F. P. Chan, R. J. Herfkens, and S. S. Vasanawala, “Rapid pediatric cardiac assessment of flow and ventricular volume with compressed sensing parallel imaging volumetric cine phase-contrast MRI,” American Journal of Roentgenology, vol. 198, pp. W250–W259, mar 2012.
  * [10] K. G. Hollingsworth, D. M. Higgins, M. McCallum, L. Ward, A. Coombs, and V. Straub, “Investigating the quantitative fidelity of prospectively undersampled chemical shift imaging in muscular dystrophy with compressed sensing and parallel imaging reconstruction,” Magnetic Resonance in Medicine, vol. 72, pp. 1610–1619, dec 2014.
  * [11] T. Akasaka, K. Fujimoto, T. Yamamoto, T. Okada, Y. Fushumi, A. Yamamoto, T. Tanaka, and K. Togashi, “Optimization of regularization parameters in compressed sensing of magnetic resonance angiography: Can statistical image metrics mimic radiologists’ perception?,” PLOS ONE, vol. 11, pp. 1–14, 01 2016.
  * [12] O. N. Jaspan, R. Fleysher, and M. L. Lipton, “Compressed sensing MRI: a review of the clinical literature,” The British Journal of Radiology, vol. 88, p. 20150487, dec 2015.
  * [13] P. C. Hansen, “The L-curve and its use in the numerical treatment of inverse problems,” 1999.
  * [14] S. Ramani, Z. Liu, J. Rosen, J.-F. Nielsen, and J. A. Fessler, “Regularization parameter selection for nonlinear iterative image restoration and MRI reconstruction using GCV and SURE-based methods,” IEEE Transactions on Image Processing, vol. 21, no. 8, pp. 3659–3672, 2012.
  * [15] R. Giryes, M. Elad, and Y. C. Eldar, “The projected GSURE for automatic parameter tuning in iterative shrinkage methods,” Applied and Computational Harmonic Analysis, vol. 30, no. 3, pp. 407–422, 2011.
  * [16] D. S. Weller, S. Ramani, J.-F. Nielsen, and J. A. Fessler, “Monte Carlo SURE-based parameter selection for parallel magnetic resonance imaging reconstruction,” Magnetic Resonance in Medicine, vol. 71, no. 5, pp. 1760–1770, 2014.
  * [17] K. Khare, C. J. Hardy, K. F. King, P. A. Turski, and L. Marinelli, “Accelerated MR imaging using compressive sensing with no free parameters,” Magnetic Resonance in Medicine, vol. 68, pp. 1450–1457, nov 2012.
  * [18] G. Varela-Mattatall, C. A. Baron, and R. S. Menon, “Automatic determination of the regularization weighting for wavelet-based compressed sensing MRI reconstructions,” Magnetic Resonance in Medicine, 2021.
  * [19] C. M. Stein, “Estimation of the Mean of a Multivariate Normal Distribution,” The Annals of Statistics, vol. 9, pp. 1135–1151, nov 1981\.
  * [20] D. L. Donoho and I. M. Johnstone, “Adapting to Unknown Smoothness via Wavelet Shrinkage,” Journal of the American Statistical Association, vol. 90, p. 1200, dec 1995.
  * [21] F. Ong, M. Uecker, U. Tariq, A. Hsiao, M. T. Alley, S. S. Vasanawala, and M. Lustig, “Robust 4D Flow Denoising Using Divergence-Free Wavelet Transform,” Magnetic Resonance in Medicine, vol. 73, p. 828, feb 2015\.
  * [22] D. L. Donoho, A. Maleki, and A. Montanari, “Message-passing algorithms for compressed sensing,” Proceedings of the National Academy of Sciences, vol. 106, no. 45, pp. 18914–18919, 2009.
  * [23] J. Ma and L. Ping, “Orthogonal AMP,” IEEE Access, vol. 5, pp. 2020–2033, 2017.
  * [24] S. Rangan, P. Schniter, and A. K. Fletcher, “Vector Approximate Message Passing,” IEEE Transactions on Information Theory, pp. 1–1, 2019.
  * [25] A. Mousavi, A. Maleki, and R. G. Baraniuk, “Parameterless Optimal Approximate Message Passing,” arXiv preprint arXiv:1311.0035, oct 2013.
  * [26] C. Guo and M. E. Davies, “Near Optimal Compressed Sensing Without Priors: Parametric SURE Approximate Message Passing,” IEEE Transactions on Signal Processing, vol. 63, pp. 2130–2141, apr 2015.
  * [27] K. Sung, B. L. Daniel, and B. A. Hargreaves, “Location Constrained Approximate Message Passing for Compressed Sensing MRI,” Magnetic Resonance in Medicine, vol. 70, pp. 370–381, 2013.
  * [28] E. M. Eksioglu and A. K. Tanc, “Denoising AMP for MRI Reconstruction: BM3D-AMP-MRI,” SIAM Journal on Imaging Sciences, vol. 11, pp. 2090–2109, jan 2018.
  * [29] C. A. Metzler, A. Maleki, and R. G. Baraniuk, “From Denoising to Compressed Sensing,” IEEE Transactions on Information Theory, vol. 62, no. 9, pp. 5117–5144, 2016.
  * [30] H. Kim, D. U. Kang, and S. Y. Chun, “Preliminary Studies On Training And Fine-Tuning Deep Denoiser Neural Networks In Learned D-Amp For Undersampled Real Mr Measurements,” in 2020 IEEE 17th International Symposium on Biomedical Imaging Workshops (ISBI Workshops), pp. 1–4, 2020.
  * [31] P. Schniter, S. Rangan, and A. Fletcher, “Plug-and-play Image Recovery using Vector AMP,” BASP Frontiers Workshop 2017, 2017.
  * [32] S. Sarkar, R. Ahmad, and P. Schniter, “MRI Image Recovery using Damped Denoising Vector AMP,” no. 2, pp. 1–5, 2020.
  * [33] C. Millard, A. T. Hess, B. Mailhe, and J. Tanner, “Approximate Message Passing with a Colored Aliasing Model for Variable Density Fourier Sampled Images,” IEEE Open Journal of Signal Processing, p. 1, 2020.
  * [34] C. Millard, A. T. Hess, B. Mailhe, and J. Tanner, “An Approximate Message Passing Algorithm For Rapid Parameter-Free Compressed Sensing MRI,” in 2020 IEEE International Conference on Image Processing (ICIP), pp. 91–95, 2020\.
  * [35] P. Virtue and M. Lustig, “The Empirical Effect of Gaussian Noise in Undersampled MRI Reconstruction,” Tomography, vol. 3, pp. 211–221, dec 2017.
  * [36] C. A. Metzler and G. Wetzstein, “D-VDAMP: Denoising-based Approximate Message Passing for Compressive MRI,” arXiv:2010.13211, 2020.
  * [37] F. Huang, S. Vijayakumar, Y. Li, S. Hertel, and G. R. Duensing, “A software channel compression technique for faster reconstruction with many channels,” Magnetic Resonance Imaging, vol. 26, no. 1, pp. 133–141, 2008\.
  * [38] M. Bayati, M. A. Erdogdu, and A. Montanari, “Estimating LASSO risk and noise level,” Advances in Neural Information Processing Systems, pp. 1–9, 2013\.
  * [39] S. Rangan, P. Schniter, and A. Fletcher, “On the convergence of approximate message passing with arbitrary matrices,” in 2014 IEEE International Symposium on Information Theory, pp. 236–240, IEEE, jun 2014.
  * [40] J. Zbontar, F. Knoll, A. Sriram, T. Murrell, Z. Huang, M. J. Muckley, A. Defazio, R. Stern, P. Johnson, M. Bruno, et al., “fastMRI: An open dataset and benchmarks for accelerated MRI,” arXiv preprint arXiv:1811.08839, 2018.
  * [41] M. Uecker, J. I. Tamir, F. Ong, and M. Lustig, “The BART Toolbox for Computational Magnetic Resonance Imaging,” in ISMRM, 2016.
  * [42] A. Beck and M. Teboulle, “A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems,” SIAM Journal on Imaging Sciences, vol. 2, pp. 183–202, jan 2009.
  * [43] P. Kalavathi and V. S. Prasath, “Methods on skull stripping of MRI head scan images—a review,” Journal of digital imaging, vol. 29, no. 3, pp. 365–379, 2016.
  * [44] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: From error visibility to structural similarity,” IEEE Transactions on Image Processing, vol. 13, no. 4, pp. 600–612, 2004.
  * [45] S. Ravishankar and Y. Bresler, “MR image reconstruction from highly undersampled k-space data by dictionary learning,” IEEE transactions on medical imaging, vol. 30, no. 5, pp. 1028–1041, 2010.
  * [46] N. Otsu, “A threshold selection method from gray-level histograms,” IEEE transactions on systems, man, and cybernetics, vol. 9, no. 1, pp. 62–66, 1979\.
  * [47] T. Akasaka, K. Fujimoto, T. Yamamoto, T. Okada, Y. Fushumi, A. Yamamoto, T. Tanaka, and K. Togashi, “Optimization of regularization parameters in compressed sensing of magnetic resonance angiography: can statistical image metrics mimic radiologists’ perception?,” PloS one, vol. 11, no. 1, p. 146548, 2016.
  * [48] C. Metzler, A. Mousavi, and R. Baraniuk, “Learned D-AMP: Principled Neural Network based Compressive Image Recovery,” in Advances in Neural Information Processing Systems 30 (I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, eds.), pp. 1772–1783, Curran Associates, Inc., 2017.
  * [49] M. Borgerding, P. Schniter, and S. Rangan, “AMP-inspired deep networks for sparse linear inverse problems,” IEEE Transactions on Signal Processing, vol. 65, no. 16, pp. 4293–4308, 2017.
  * [50] R. Kitichotkul, C. A. Metzler, F. Ong, and G. Wetzstein, “SUREMap: Predicting Uncertainty in CNN-based Image Reconstruction Using Stein’s Unbiased Risk Estimate,” arXiv preprint arXiv:2010.13214, 2020.
  * [51] C. Millard, A. Hess, J. Tanner, and B. Mailhe, “Deep Plug-and-Play multi-coil compressed sensing MRI with matched aliasing: the Denoising-P-VDAMP algorithm,” in 2022 ISMRM annual meeting, 2022.
  * [52] C. Millard, Approximate message passing for compressed sensing magnetic resonance imaging.  PhD thesis, University of Oxford, 2021.
  * [53] S. V. Venkatakrishnan, C. A. Bouman, and B. Wohlberg, “Plug-and-Play priors for model based reconstruction,” in 2013 IEEE Global Conference on Signal and Information Processing, pp. 945–948, IEEE, dec 2013.
  * [54] R. Ahmad, C. A. Bouman, G. T. Buzzard, S. Chan, S. Liu, E. T. Reehorst, and P. Schniter, “Plug-and-Play Methods for Magnetic Resonance Imaging: Using Denoisers for Image Recovery,” IEEE Signal Processing Magazine, vol. 37, pp. 105–116, jan 2020.
  * [55] K. M. Johnson, S. B. Fain, M. L. Schiebler, and S. Nagle, “Optimized 3D ultrashort echo time pulmonary MRI.,” Magnetic resonance in medicine, vol. 70, pp. 1241–1250, nov 2013.
  * [56] P. T. Gurney, B. A. Hargreaves, and D. G. Nishimura, “Design and analysis of a practical 3D cones trajectory,” Magnetic Resonance in Medicine, vol. 55, no. 3, pp. 575–582, 2006.
  * [57] S. K. Shastri, R. Ahmad, C. A. Metzler, and P. Schniter, “Expectation Consistent Plug-and-Play for MRI,” arXiv preprint arXiv:2202.05820, 2022\.


††thanks: X. Chen and Z. Byambadorj contributed equally to this work.††thanks:
X. Chen and Z. Byambadorj contributed equally to this work.

# CMOS-based area-and-power-efficient neuron and synapse circuits for time-
domain analog spiking neural networks

Xiangyu Chen  Systems Design Lab., School of Engineering, the University of
Tokyo, Tokyo, Japan  Zolboo Byambadorj  Systems Design Lab., School of
Engineering, the University of Tokyo, Tokyo, Japan  Takeaki Yajima  Department
of Electrical and Electronic Engineering, Kyushu University, Fukuoka, Japan
Hisashi Inoue  National Institute of Advanced Industrial Science and
Technology (AIST), Ibaraki, Japan  Isao H. Inoue  National Institute of
Advanced Industrial Science and Technology (AIST), Ibaraki, Japan  Tetsuya
Iizuka<EMAIL_ADDRESS>Systems Design Lab., School of Engineering,
the University of Tokyo, Tokyo, Japan

###### Abstract

Conventional neural structures tend to communicate through analog quantities
such as currents or voltages, however, as CMOS devices shrink and supply
voltages decrease, the dynamic range of voltage/current-domain analog circuits
becomes narrower, the available margin becomes smaller, and noise immunity
decreases. More than that, the use of operational amplifiers (op-amps) and
continuous-time or clocked comparators in conventional designs leads to high
energy consumption and large chip area, which would be detrimental to building
spiking neural networks. In view of this, we propose a neural structure for
generating and transmitting time-domain signals, including a neuron module, a
synapse module, and two weight modules. The proposed neural structure is
driven by a leakage current of MOS transistors and uses an inverter-based
comparator to realize a firing function, thus providing higher energy and area
efficiency compared to conventional designs. The proposed neural structure is
fabricated using TSMC 65 nm CMOS technology. The proposed neuron and synapse
occupy the area of 127 $\bm{\mu m^{2}}$ and 231 $\bm{\mu m^{2}}$,
respectively, while achieving millisecond time constants. Actual chip
measurements show that the proposed structure implements the temporal signal
communication function with millisecond time constants, which is a critical
step toward hardware reservoir computing for human-computer interaction.
Simulation results of the spiking-neural network for reservoir computing with
the behavioral model of the proposed neural structure demonstrate the learning
function.

††preprint: AIP/123-QED

Deep neural networks (DNNs), which are the second generation of artificial
neural networks (ANNs), have extensively explored in recent years for growing
number of applications. However, their huge energy consumption especially for
the memory access in conventional von-Neumann architecture has forced people
to find an alternative way to achieve more power-efficient solutions Zhang
_et al._ (2020); Shin and Yoo (2020); LeCun, Bengio, and Hinton (2015); Kohno
and Aihara (2008); Chicca and Indiveri (2020); Bo  _et al._ (2020). Spiking
neural network (SNN) is one of the attractive solutions as the third
generation of ANNs that can realize learning function with low power by
mimicking biological neurons. SNNs consist of neurons and synapses, and are
usually built using a bottom-up approach, which means that each component of
the SNNs needs to be designed firstBo  _et al._ (2020); Yang and Sengupta
(2020); Chen  _et al._ (2022); Maass (1997); Radhakrishnan  _et al._ (2021);
Chen  _et al._ (2021); Jeong (2018).

Many hardware implementations of pulsed neurons or synapses have been
reportedIndiveri, Chicca, and Douglas (2006); Wu  _et al._ (2015); Joubert
_et al._ (2012); Aamir  _et al._ (2018a); Basu and Hasler (2010); Dutta  _et
al._ (2017); Rubino, Payvand, and Indiveri (2019); Aamir  _et al._ (2018b);
Saber Moradi (2017). To implement the leaky integrate function of neurons,
conventional designs usually build integrators with operational amplifiers
(op-amps)Wu  _et al._ (2015) and often use large on-chip capacitors and
resistors to mimic the millisecond time constants of biological neuronsBasu
and Hasler (2010); Aamir  _et al._ (2018a). Moreover, to implement the neuron
“fire” function, a dedicated circuit structure of a continuous-time or clocked
comparator is usually used to set the threshold for neuron excitationAamir
_et al._ (2018b); Wu  _et al._ (2015); Joubert  _et al._ (2012); Indiveri,
Chicca, and Douglas (2006); Aamir  _et al._ (2018a). The bias current of the
continuous-time comparator undoubtedly increases the power consumption of the
neuron, while the clocked comparator requires additional clock signal
distribution and the complex comparator structure occupies a large chip area.
While more advanced processes can achieve low power consumption by reducing
supply voltage and static leakage currentSaber Moradi (2017), this also leads
to a narrower dynamic range, smaller available margin, and degraded noise
immunity of the voltage/current-domain analog circuitsAsada  _et al._ (2018).
This is detrimental to conventional neural networks that use analog quantities
such as voltage and current to communicate with each other. On the other hand,
thanks to the scaled transistors that have an improved operation speed with
sharp signal transitions, the analog information can be represented more
efficiently in time domain, i.e. a time interval of two signal transitions.
This so-called time-domain circuit have another advantage in its power
efficiency as it often consists of inverters or logic gates that ideally
consume no DC power Staszewski  _et al._ (2004); Asada  _et al._ (2018). Thus,
time-domain circuits are ideal for future implementations of low-power SNNs.

In this paper we propose an original neural structure for generating and
transmitting time-domain signals to compose a time-domain neural network. The
integrated structure includes neuron and synapse modules that respectively
generate and transmit time-domain signals, as well as weight modules for
learning functions. One of our main target applications is reservoir
computing, which processes information related to human activity. Our
application targets simpler and less data-intensive processing such as bio-
signals. In reservoir computing, learning functions such as ECG and speaker
recognition as well as handwriting recognition can be implemented using only a
few hundreds of neurons. Ref. Gallicchio and Micheli, 2011 shows that learning
performance improves when the time constants of the input effects are matched
between the target function and reservoir dynamics, and we use millisecond
time constants as a design target for a neural structure that will be used to
process time-series information of human activities. We use the behavioral
model of the proposed neural structure to construct the SNN for reservoir
computing and implement the learning function, which proves that our proposed
neural structure can be used for reservoir computing.

(a)

(b)

Figure 1: (a) The proposed structure and (b) a micrograph of the chip.

The designed and fabricated neural structure is shown in Fig. 1(a), which is
based on the proposed neuron, synapse, and weight modules, which will be
described in detail below. In this structure, the input of the neuron module
is connected to two weight modules, one for tuning the inhibitory signal and
the other for the excitatory signal. We fabricated the proposed neural
structure shown in Fig. 1(a) with TSMC 65 nm standard CMOS technology. The
micrograph of the chip is shown in Fig. 1(b), where the die area of neuron,
synapse and weight modules are 127 $\rm{\mu m^{2}}$, 231 $\mu m^{2}$ and 525
$\mu m^{2}$, respectively.

The LIF neuron model consists mainly of a membrane capacitor, a leaky
resistor, and a voltage comparator. Neurons receive signals from other neurons
via synapses and the soma generates action potentials in response to these
external signals. If a neuron receives a sufficient number of spikes through
the synapse, its membrane potential will reach a threshold value, causing the
neuron to “fire”Abbott and Dayan (2005); Gerstner and Kistler (2012); Chen
_et al._ (2022).

(a)

(b)

(c)

(d)

Figure 2: (a) Circuit diagram of the proposed neuron module, (b) behaviours of
proposed LIF neuron and synapse modules, (c) circuit diagram of the proposed
synapse module, and (d) circuit diagram of the proposed weight module.

The use of inverters to implement the “fire” function is already known as an
alternative to comparators. Ref. Yajima, 2022 has proposed an inverter-based
neuron, which is well suited for use in the proposed neural structure, and
therefore the neuron used in this study was designed based on Ref. Yajima,
2022, which is shown in Fig. 2(a). It consists of input device, leaky
integrator device, fire device and delay device. Originally in Ref. Yajima,
2022, the circuit is not assumed to be designed as an element to build a
neural network, thus does not have a structure to receive excitatory and
inhibitory signals. In the proposed circuit, on the other hand, the input
device that consists of $\mathrm{M}_{\mathrm{1}}$ and
$\mathrm{M}_{\mathrm{2}}$ receives an excitatory input and an inhibitory
input, respectively. The inputs to $\mathrm{M}_{\mathrm{1}}$ and
$\mathrm{M}_{\mathrm{2}}$ are narrow pulse signals as shown in Fig. 2(a),
which is generated from a pre-stage synapse. The activity of the pre-stage
synapse is represented by the pulse frequency, and the coupling weight is
represented by the pulse width. When more than one pre-stage synapses are
connected to compose a network, the multiple pulses can be applied through OR
logic, or by adding input devices connected in parallel. With the parallel
input devices, the neuron circuit can accept multiple pulses even at the same
time.

In leaky-integrator device, $C_{\mathrm{mem}}$ represents the cell membrane of
the neuron, and $\mathrm{M}_{\mathrm{5}}$ can be regarded as a leaky resistor
in the resting state. When there is no external input to the input device, the
capacitor is charged by the leakage current of $\mathrm{M}_{\mathrm{3}}$ and
$\mathrm{M}_{\mathrm{4}}$, and the membrane potential $V_{\mathrm{mem}}$
increases continuously with the inflow of the leakage current (the current is
integrated as shown in Fig. 2(b)(i)). At this point, since
$\mathrm{M}_{\mathrm{5}}$ is in off state, it can be considered as a resistor
in parallel with the capacitor, i.e. leaky resistor, capable of achieving a
long time constant.

Once $V_{\mathrm{mem}}$ rises to the threshold voltage
$V_{\mathrm{th(Fire)}}$, the fire device is activated (Fig. 2(b)(ii)). In
conventional designs, LIF neurons mostly use dedicated circuit structures of
continuous-time or clocked comparator to set the threshold voltage. This is
not friendly for building SNNs that are as energy efficient and bio-scale as
the brain. In this study, the fire device is implemented by an inverter-based
comparator that can set the threshold voltage by two transistors instead of a
continuous-time or clocked comparator. To realize an accurate threshold
voltage for an inverter-based comparator, we may use an auto-zeroing technique
that periodically sense, store and cancel the offset with switches and
capacitors Razavi (1995). However, it requires multiple-phase clocks to
control the switches, thus is not suitable for area-and-power-efficient
reservoir implementations. Though with a simple inverter-based comparator
there may be a threshold variation due to process, voltage and temperature
fluctuations, it can be seen as mimicking the difference between individuals
of real neurons. In addition, the learning function is able to compensate for
threshold differences and process variations Wunderlich  _et al._ (2019). When
there is an excitatory pulse input, $\mathrm{M}_{\mathrm{1}}$ will be turned
on instantaneously, which causes more current to charge $C_{\mathrm{mem}}$ and
$V_{\mathrm{mem}}$ to rise rapidly. Conversely, an inhibitory pulse input
signal will cause $\mathrm{M}_{\mathrm{2}}$ to turn on momentarily, causing
$C_{\mathrm{mem}}$ to charge slower or even discharge through
$\mathrm{M}_{\mathrm{2}}$, which in turn slows down the rate of
$V_{\mathrm{mem}}$ rise or makes it fall.

(a)

(b)

(c)

(d)

Figure 3: (a) A photo of the experimental setup, (b) the measured firing rate
of the neuron for 4 chips, (c) the measured waveforms of the neuron output,
and (d) the measured waveforms of the synapse output.

(a)

(b)

Figure 4: (a) Another combined structure fabricated to evaluate the synapse
and (b) the measured waveforms of $V_{\mathrm{Ring}}$ and $V_{\mathrm{SYN}}$.

(a)

(b)

Figure 5: (a) The behavioral model of the SNN for reservoir computing based on
the proposed neural structure. (b) The system-level behavioral simulation
results: (i) based on a model with $15$ Hz$-200$ Hz frequency tuning range,
zoomed-in view of the (ii) excitatory and (iii) inhibitory input signals
converted from the output, (iv) based on $15$ Hz$-2$ kHz and (v) $15$ Hz$-20$
kHz frequency tuning ranges.

When the fire device is activated, it generates a low level of
$V_{\mathrm{Fire}}$ to be connected to $\mathrm{M}_{\mathrm{4}}$, which will
increase the current to charge the membrane capacitor $C_{\mathrm{mem}}$,
resulting in an instantaneous increase of the membrane potential
$V_{\mathrm{mem}}$, which promotes the triggering of the fire device. This
mimics the influx of Na+ into the cell membrane prompting a rapid increase in
membrane voltage, i.e., a positive feedback effect. Finally, the low level of
$V_{\mathrm{Fire}}$ generated from the fire device is converted to a high
level of $V_{\mathrm{Spike}}$ (Fig. 2(b)(iii)) by a delay device that includes
a three-stage inverter and connects the $V_{\mathrm{Spike}}$ to
$\mathrm{M}_{\mathrm{3}}$ and $\mathrm{M}_{\mathrm{5}}$, resetting
$V_{\mathrm{mem}}$ to zero. This process mimics the activation of K+ channels
in biological neurons, resulting in the outward flow of K+ ions and the
eventual return of the cell membrane to its resting state.

Synapses are essential modules in SNNs, as neurons are interconnected by them.
We have designed a neuron module for generating time-domain signals, and then
we need a transmission medium, i.e., a synapse, to transmit this time-domain
signal to other neurons. To compose a complete neural network, we design a
synapse module based on frequency signals, as shown in Fig. 2(c). The synapse
consists mainly of a voltage-controlled ring oscillator operating under a
leakage current, which is composed of a three-stage inverter
($\mathrm{M}_{\mathrm{6}},\mathrm{M}_{\mathrm{7}},\mathrm{M}_{\mathrm{8}},\mathrm{M}_{\mathrm{9}},\mathrm{M}_{\mathrm{10}}$
and $\mathrm{M}_{\mathrm{11}}$). The previous neuron circuit fires and
generates a spike $V_{\mathrm{Spike}}$, which is inverted by an inverter,
making $\mathrm{M}_{\mathrm{5}}$ open for a short time, and the current
flowing through $\mathrm{M}_{\mathrm{5}}$ charges $C_{\mathrm{SYN}}$, which
will increase $V_{\mathrm{SYN}}$. Once $V_{\mathrm{SYN}}$ reaches the voltage
that triggers the oscillation, the ring oscillator begins to oscillate (Fig.
2(b)(iv) and Fig. 2(b)(v)). If the preceding neuron does not fire for a long
time, $V_{\mathrm{SYN}}$ will leak until the initial state, at which point the
synapse becomes inactive again. Since $V_{\mathrm{SYN}}$ is equivalent to the
supply voltage of the ring oscillator, the current flowing out of
$\mathrm{M}_{\mathrm{5}}$ controls $V_{\mathrm{SYN}}$ and thus the frequency
of the ring oscillator.

SNNs achieve learning function by adjusting the weights; therefore, we propose
a weight module that is compatible with the proposed time-domain neuron and
synapse modules described above, as shown in Fig. 2(d). The proposed weight
module tunes the time-domain information, which is the width of the output
pulses. This module consists of a delay line, a multiplexer, and an AND gate.
$V_{\mathrm{Ring}}$ is the square wave signal from the synapse that will pass
through the delay line. $V_{\mathrm{Weight}}$ is the digital code that
represents weight, which is determined after learning and is used to control
the multiplexer. The width of the output pulse that corresponds to the time-
domain weight is adjusted according to which tap in the inverter chain is
selected by the multiplexer. As mentioned earlier, if the excitatory or
inhibitory pulse width is wide, the voltage $V_{\mathrm{mem}}$ in the
subsequent neuron is charged or discharged faster, respectively. This
corresponds to a large weight. In this study, we chose a multiplexer with 16
inputs, i.e. four bit weights (0000 to 1111). The output of the weight module
is connected to the input device of the subsequent neuron circuits. The
frequency of the pulse (pulse spacing) and the width of the pulse act
simultaneously on the neuron to change its activity. The frequency of the
pulse is determined by the output frequency of the previous synapse, while the
coupling strength depends on the width of the pulse output determined by the
weight module.

Figure 3(a) shows the experimental setup used to test the fabricated neural
structure chip (Fig. 1(b)), where the chip was placed on a probe station
Summit11000 and tested with probes in direct contact with it. In the
experiments, we assume that the inputs of the two weight modules is the pre-
stage synapses, which is emulated by the arbitrary function generators. The
output of the neuron is connected to the synapse module, and the output of
which will be varied in response to the change in the output of the neuron. We
used a Tektronix AFG31252 arbitrary function generator as a pre-stage synapse
to provide square wave signals for our fabricated neural circuits. At the same
time, we observed the output waveforms using oscilloscopes (Keysight MSOX6004A
and DSOX93304Q).

The experimental results are shown in Figs. 3(b), 3(c) and 3(d). To verify the
effect of weights on the firing rate of neurons, we fixed the frequency of the
pre-stage synapse output (function generator) at 100 Hz and observed the
change in the firing rate of neurons for 4 chips by adjusting the weights
module. We averaged the spike frequencies 1024 times over a time range of 100
ms to derive the corresponding neuronal firing frequency under each weight
setting, as shown in Fig. 3(b). The proposed neuron is basically firing with
the rate determined by the leakage currents into and out from
$C_{\mathrm{mem}}$ in balance, and input from the previous stage modulates it.
We can see that when the weights become larger, the neuron module firing
frequency becomes larger. Mainly due to the process variation of the FETs, the
firing frequency fluctuates in about $\pm$10–17 % over 4 chips. Especially for
the use in a reservoir, however, due to the random weights in its recurrent
connections, these random variations should be compensated for during the
learning process in the output weights.

Figure 3(c) compares the variation of neuron fire times depending on the
signal from the pre-stage synapse. The insets (i), (ii) and (iii) of Fig. 3(c)
show the cases with 100 Hz inhibitory input (weight is set to 1100), no input,
and with 100 Hz excitatory input (weight is set to 1100), respectively, from
which we can see that the inhibitory input decreases the fire frequency of the
neuron and increases the fire interval, while the excitatory input works as
the opposite of the inhibitory input. The experimental results show that the
firing interval of the proposed neuron is on the order of milliseconds, which
is in accordance with the feature of biological neurons having millisecond
time constants. When no signal is fed from the pre-stage synapse, the power
consumption is about 800 pW, generating about 20 spikes in a 100 ms cycle.
From this, it can be roughly estimated that each spike consumes about 4 pJ of
energy. Subsequently, the insets (i), (ii) and (iii) of Fig. 3(c) were used as
input signals to the synapse to influence $V_{\mathrm{Ring}}$. The measured
$V_{\mathrm{Ring}}$ waveforms in these three cases are shown in Fig. 3(d). The
average of the frequencies for each case measured in 5 s time period are 41
Hz, 90 Hz and 98 Hz, respectively. The feasibility of this synapse output
frequency range will be validated with system-level simulations in the
following discussion.

To facilitate the observation of the synchronous response of the synapse, we
also fabricated the structure of Fig. 4(a). Figure 4(b) is the experimental
results on Fig. 4(a). We used a Tektronix AFG31252 arbitrary function
generator to generate a 10 Hz square wave signal $V_{\mathrm{IN}}$ as shown in
Fig. 4(a) (i). After $V_{\mathrm{IN}}$ passes through a weight module, it
produces a spike signal $V_{\mathrm{OUT(Weight)}}$. The voltage
$V_{\mathrm{SYN}}$ is observed through an on-chip source follower as an analog
buffer. Though $V_{\mathrm{OUT(Weight)}}$ is not designed to be observed from
outside as it is a narrow pulse, with the arrival of the
$V_{\mathrm{OUT(Weight)}}$ after the falling edge of $V_{\mathrm{IN}}$, the
$V_{\mathrm{SYN}}$ voltage at the synapse rises instantaneously as shown in
Fig. 4(b)(ii), which in turn increases the frequency of $V_{\mathrm{Ring}}$.
If the $V_{\mathrm{OUT(Weight)}}$ does not arrive for a long time,
$V_{\mathrm{SYN}}$ decreases, which in turn affects the $V_{\mathrm{Ring}}$
frequency to become smaller.

Table 1: Performance Comparison of Stand-Alone Neuron Circuits. Ref. | Technology |  Energy per spike (pJ) |  Single neuron area ($\mu m^{2}$) |  Frequency (Hz) |  Neuron model |  Sim. or Meas.
---|---|---|---|---|---|---
Indiveri, Chicca, and Douglas,2006 |  350 nm CMOS | 900 | 2573 | 100 | IF | Meas.
Wu  _et al._ ,2015 |  180 nm CMOS | 9.3 | $10^{4}$ | N/A | LIF | Meas.
Joubert  _et al._ ,2012 |  65 nm CMOS | 41 | 538 | 300 | LIF | Sim.
Aamir  _et al._ ,2018a |  65 nm CMOS | 200 | 3363 | $1.9\times 10^{6}$ | AdEx-IF111Adaptive-exponential integrate-and-fire. | Meas.
Dutta  _et al._ ,2017 |  32 nm SOI MOSFET | 35 | 1.8 | $10^{6}$ | LIF | Meas.
Rubino, Payvand, and Indiveri,2019 |  22 nm FD-SOI | 14 | 900 | 30 | AdEx-IF111Adaptive-exponential integrate-and-fire. | Meas.
Saber Moradi,2017 |  22 nm CMOS | 1.17/0.36 | 70 | 30/100 | LIF | Sim.
this work |  65 nm CMOS | 4 | 127 | 230 | LIF | Meas.

Table 1 shows the performance comparison among stand-alone neuron circuits.
The proposed neuron circuit has advantages in terms of energy consumption and
area. The designs in Refs. Indiveri, Chicca, and Douglas, 2006; Wu  _et al._ ,
2015; Joubert  _et al._ , 2012; Aamir  _et al._ , 2018a used a continuous-time
or a clocked comparator, and these designs take up a large amount of chip area
as well as power consumption. The neuron fabricated in a non-CMOS process
proposed in Ref. Dutta  _et al._ , 2017 does not require a comparator, which
leads to an advantage in area. However, its energy consumption is relatively
high and these particular technologies are less mature and thus more costly
compared to standard CMOS processes. Both Ref. Rubino, Payvand, and Indiveri,
2019 and Ref. Saber Moradi, 2017 are being fabricated in an advanced process.
However, compared to this work Ref. Rubino, Payvand, and Indiveri, 2019 does
not have an advantage in terms of energy consumption and area. Although Ref.
Saber Moradi, 2017 shows better energy efficiency with simulation results,
when normalized by the technology node, the proposed neuron achieves better
area efficiency.

To demonstrate the feasibility of the proposed spiking neuron and the ring
oscillator-based synapse circuits, a behavioral simulation is carried out in
MATLAB environment as shown in Fig. 5(a). In this simulation, 100 neurons are
used with random recurrent connections with the proposed synapse and weighting
modules. The proposed weight modules are applied only in the reservoir layer
and their weights are assigned randomly in advance and fixed during the
learning process. Thus the random fluctuations in the reservoir are
compensated for during the learning process in the output weights. To
establish a realistic simulation, the output frequency range of each synapse
is set from $15$ Hz to $200$ Hz based on the actual measurement results. The
recursive least square (RLS) algorithm is used to train the output weights as
introduced in Ref. Sussillo and Abbott, 2009. A $10$ Hz sinewave, which
corresponds to the time-scale of human-activity-related information, is used
as an example of supervisory input signal. The supervisory and the trained
output signal are shown in Fig. 5(b)(i). The feedback signal from the output
is converted to excitatory and inhibitory pulse trains whose frequencies are
in proportion to the absolute value of the output amplitude as shown in Figs.
5(b)(ii) and (iii), respectively. After 5 periods of supervisory signal, the
output weights are fixed and the SNN generates the learned signal by itself,
which demonstrates the feasibility of the proposed neural structures for
learning function. We have also found from these simulations that to further
improve the learning capability the output frequency tuning range of the
synapse should be increased, which can be done by optimizing the synapse
circuit. For example, with the extended frequency tuning ranges from $15$
Hz$-2$ kHz and $15$ Hz$-20$ kHz, the learned signals become smoother to better
reproduce the supervisory signal as shown in Figs. 5(b)(iv) and 5(b)(v),
respectively.

In summary, we have proposed a neural structure for generating and
transmitting time-domain signals. The proposed neuron and synapse occupy an
area of 127 ${\mu m^{2}}$ and 231 ${\mu m^{2}}$, respectively. This structure
does not use op-amps and continuous-time or clocked comparators, while the
firing function is realized with an inverter-based comparator to provide
advantages in area and power consumption. The proposed time-domain neural
structure benefits from scaled process technologies compared to conventional
voltage/current-domain designs. Actual chip fabrication and measurement
results demonstrate the temporal signal communication function with
millisecond time constants. The proposed time-domain neural structure is well
suited for building spiking neural networks for processing real-time time-
series information for human-computer interaction.

This work is supported by Japan Science and Technology Agency (JST) CREST
Grant Number JPMJCR19K2.

## References

  * Zhang  _et al._ (2020) Y. Zhang, P. Qu, Y. Ji, W. Zhang, G. Gao, G. Wang, S. Song, G. Li, W. Chen, W. Zheng, F. Chen, J. Pei, R. Zhao, M. Zhao, and L. Shi, Nature 586, 378–384 (2020).
  * Shin and Yoo (2020) D. Shin and H.-J. Yoo, Proceedings of the IEEE 108, 1245–1260 (2020).
  * LeCun, Bengio, and Hinton (2015) Y. LeCun, Y. Bengio, and G. Hinton, Nature 521, 436–444 (2015).
  * Kohno and Aihara (2008) T. Kohno and K. Aihara, AIP Conference Proceedings 1028, 113–128 (2008).
  * Chicca and Indiveri (2020) E. Chicca and G. Indiveri, Appl. Phys. Lett. 116, 120501 (2020).
  * Bo  _et al._ (2020) Y. Bo, P. Zhang, Y. Zhang, J. Song, S. Li, and X. Liu, J. Appl. Phys. 127, 245101 (2020).
  * Yang and Sengupta (2020) K. Yang and A. Sengupta, Appl. Phys. Lett. 116, 043701 (2020).
  * Chen  _et al._ (2022) X. Chen, T. Yajima, I. H. Inoue, and T. Iizuka, Jpn. J. Appl. Phys. 61, SC1051 (2022).
  * Maass (1997) W. Maass, Neural. Netw. 10, 1659–1671 (1997).
  * Radhakrishnan  _et al._ (2021) S. S. Radhakrishnan, A. Sebastian, A. Oberoi, S. Das, and S. Das, Nat. Commun. 12, 2143 (2021).
  * Chen  _et al._ (2021) X. Chen, T. Yajima, I. H. Inoue, and T. Iizuka, Extended Abstract of International Conference on Solid State Devices and Materials (SSDM) , 682–683 (2021).
  * Jeong (2018) D. S. Jeong, J. Appl. Phys. 124, 152002 (2018).
  * Indiveri, Chicca, and Douglas (2006) G. Indiveri, E. Chicca, and R. Douglas, IEEE Trans. Neural. Netw. 17, 211–221 (2006).
  * Wu  _et al._ (2015) X. Wu, V. Saxena, K. Zhu, and S. Balagopal, IEEE Trans. Circuits Syst. II, Express Briefs 62, 1088–1092 (2015).
  * Joubert  _et al._ (2012) A. Joubert, B. Belhadj, O. Temam, and R. Héliot, The 2012 International Joint Conference on Neural Networks (IJCNN) , 1–5 (2012).
  * Aamir  _et al._ (2018a) S. A. Aamir, P. Müller, G. Kiene, L. Kriener, Y. Stradmann, A. Grübl, J. Schemmel, and K. Meier, IEEE Trans. Biomed. Circuits Syst. 12, 1027–1037 (2018a).
  * Basu and Hasler (2010) A. Basu and P. E. Hasler, IEEE Trans. Circuits Syst. I, Reg. Papers 57, 2938–2947 (2010).
  * Dutta  _et al._ (2017) S. Dutta, V. Kumar, A. Shukla, N. R. Mohapatra, and U. Ganguly, Scientific Reports 7, 8257 (2017).
  * Rubino, Payvand, and Indiveri (2019) A. Rubino, M. Payvand, and G. Indiveri, 2019 26th IEEE International Conference on Electronics, Circuits and Systems (ICECS) , 458–461 (2019).
  * Aamir  _et al._ (2018b) S. A. Aamir, Y. Stradmann, P. Müller, C. Pehle, A. Hartel, A. Grübl, J. Schemmel, and K. Meier, IEEE Trans. Circuits Syst. I, Reg. Papers 65, 4299–4312 (2018b).
  * Saber Moradi (2017) R. M. Saber Moradi, Sunil A. Bhave, 2017 IEEE Symposium Series on Computational Intelligence (SSCI) , 1–5 (2017).
  * Asada  _et al._ (2018) K. Asada, T. Nakura, T. Iizuka, and M. Ikeda, IEICE Electronics Express 15, 20182001–20182001 (2018).
  * Staszewski  _et al._ (2004) R. Staszewski, K. Muhammad, D. Leipold, C.-M. Hung, Y.-C. Ho, J. Wallberg, C. Fernando, K. Maggio, R. Staszewski, T. Jung, J. Koh, S. John, I. Y. Deng, V. Sarda, O. Moreira-Tamayo, V. Mayega, R. Katz, O. Friedman, O. Eliezer, E. de Obaldia, and P. Balsara, IEEE J. Solid-State Circuits 39, 2278–2291 (2004).
  * Gallicchio and Micheli (2011) C. Gallicchio and A. Micheli, Neural. Netw. 24, 440–456 (2011).
  * Abbott and Dayan (2005) L. F. Abbott and P. Dayan,  _Theoretical Neuroscience_ (The MIT Press, 2005).
  * Gerstner and Kistler (2012) W. Gerstner and W. M. Kistler,  _Spiking Neuron Models: Single Neurons, Populations_ (Cambridge University Press, 2012).
  * Yajima (2022) T. Yajima, Sci. Rep. 12, 1150 (2022).
  * Razavi (1995) B. Razavi,  _Principle of Data Conversion System Design_ (Wiley-IEEE Press, 1995).
  * Wunderlich  _et al._ (2019) T. Wunderlich, A. F. Kungl, E. Müller, A. Hartel, Y. Stradmann, S. A. Aamir, A. Grübl, A. Heimbrecht, K. Schreiber, D. Stöckel, C. Pehle, S. Billaudelle, G. Kiene, C. Mauch, J. Schemmel, K. Meier, and M. A. Petrovici, Front. Neurosci. 13, 1–15 (2019).
  * Sussillo and Abbott (2009) D. Sussillo and L. Abbott, “Generating coherent patterns of activity from chaotic neural networks,” Neuron 63, 544–557 (2009).