1 Introduction
We consider a complex data vector
formed of noisy Fourier coefficients of a signal of interest :(1) 
where
is a multidimensional discrete Fourier transform and
is a diagonal undersampling mask with on the th diagonal entry if and otherwise, where is a sampling set with for . Here, where is the identity matrix anddenotes the distribution with independent real and imaginary parts that are normally distributed with mean
and covariance matrix . A wellstudied approach is to seek a solution of(2) 
where is a modelbased penalty function. Compressed sensing [Donoho2006CompressedSensing, Candes2006RobustInformation] concerns the reconstruction of signals of interest from underdetermined measurements, where sparsity in is promoted by solving Eq. 2 with for sparse weighting and sparsifying transform .
A prominent success of compressed sensing with Fourier measurements is accelerated Magnetic Resonance Imaging (MRI) [Lustig2007, Otazo2010CombinationMRI, Jaspan2015CompressedLiterature, Ye2019CompressedPerspective, Donoho2017HowIndustry], where the signal of interest is a natural image. Images of interest typically have a highly nonuniform spectral density that is concentrated at low frequencies. Accordingly, it is wellknown that better image restoration is possible if the sampling set
is generated with variable density, so that there is a higher probability of sampling low frequencies
[Puy2011OnSamplingb, wangVDS, Krahmer2014StableImaging, Chauffert2013VariableStrategies]. This work considers anwith elements drawn independently from a Bernoulli distribution with generic nonuniform probability, so that
.1.1 Approximate Message Passing
The Approximate Message Passing (AMP) algorithm [Donoho2009MessagepassingSensing.] is an iterative method that, for certain sensing matrices , efficiently estimates in problems of the form At iteration , AMP implements a denoiser on with meansquared error estimate . For instance, for penalty function , is the proximal operator associated with :
(3) 
which is equal to soft thresholding in the case of and orthogonal . For certain sensing matrices and given mild conditions on , AMP’s state evolution guarantees that in the large system limit , , vector is the original signal corrupted by zeromean complex Gaussian noise with a covariance matrix that is proportional to the identity:
(4) 
where
is a scalar iterationdependant standard deviation. When the covariance matrix is proportional to the identity, as in
Eq. 4, the noise and state evolution are referred to in this work as white. The white state evolution of AMP has been proven for real i.i.d. Gaussian measurements in [Bayati2011TheSensing] and i.i.d. subGaussian measurements in [Bayati2015UniversalityAlgorithms]. It has also been shown empirically that it holds for uniformly undersampled Fourier measurements of an artificial i.i.d. signal [Donoho2009MessagepassingSensing.]. When state evolution holds, AMP is known to exhibit very fast convergence. However, for generic , the behavior of AMP is not well understood and it has been noted by a number of authors [Rangan2014OnMatrices, Caltagirone2014OnPassingb, Guo2015ApproximateTransformation, Rangan2016FixedMatrices] that it can encounter convergence problems. The recent Orthogonal AMP (OAMP) [Ma2017OrthogonalAMP] and related Vector Approximate Message Passing (VAMP) [Rangan2019VectorPassing] algorithm obey a white state evolution for a broader class of measurement matrices , and were found to perform very well on certain reconstruction tasks. For VAMP, white state evolution was proven for sensing matrices that are ‘rightorthogonally invariant’: see [Rangan2019VectorPassing] for details.1.2 Colored aliasing
AMP, OAMP and VAMP assume that the entries of
are i.i.d. and that the sensing matrix is sufficiently random to ensure that the aliasing is white. However, when an image is sampled in the Fourier domain the aliasing is innately colored. To see this, consider the natural initialization for an approximate message passing algorithm: the unbiased estimator
, where is the diagonal matrix formed from sampling probabilities . Denoting , the power spectrum of the aliasing of is equationparentequation(5a)  
(5b)  
(5c) 
where is the entrywise absolute squared and is the dimensional vector of ones. Equation Eq. 5 depends on and , which are nonuniform and anisotropic in general. Note that although the specific case for constant does lead to white aliasing, it requires knowledge of the ground truth spectral density so is not feasible in practice.
A visual example of colored aliasing is shown in Fig. 1. Here, the unbiased estimate of a 512x512 synthetic SheppLogan uniformly sampled with for all is shown, with . The power spectrum of the aliasing is proportional to the power spectrum of the signal, so is highly colored. Sampling with variable density can partially but not entirely compensate for the power spectrum of the signal.
causes the horizontal, vertical and diagonal variances to differ, even at the same scale.
The intrinsically colored aliasing of variable density sampling from a nonuniform spectral density implies that the white state evolution of AMP Eq. 4 is not feasible. The primary development of this work is based on the use of the Discrete Wavelet Transform (DWT) to compute a multiresolution decomposition of the power spectrum of the aliasing [Goossens2011WaveletBasedNoise, Johnstone1997WaveletNoise]. In the wavelet domain colored aliasing has a structure that resembles a state evolution. To illustrate this, consider again the unbiased initialization . The corresponding estimator in the wavelet domain is
(6) 
Fig. 2 shows , and the entrywise absolute difference for the same image and sampling set as Fig. 1, where is a Haar wavelet with 4 decomposition scales. Fig. (c)c indicates that the power spectrum given by Eq. 5 has a simple structure in the wavelet domain; in particular, each wavelet subband behaves as the ground truth subband corrupted by zeromean white Gaussian noise. The variance of the Gaussian noise depends on the power spectrum of the aliasing and is different for each subband, even for subbands at the same scale.
1.3 Colored state evolution
Herein we present a new method for undersampled signal reconstruction that we term the Variable Density Approximate Message Passing (VDAMP) algorithm, see Algorithm 1. We present empirical evidence that VDAMP preserves the subbanddependent noise structure illustrated in Fig. (c)c for all iterations. Explicitly, the of VDAMP behaves as
(7) 
where is an orthogonal DWT and the covariance matrix is diagonal so that for a with decomposition scales
(8) 
where and refer to the variance and dimension of the th subband respectively. We refer to Eqs. 8 and 7 as the colored state evolution of VDAMP and as the effective noise of VDAMP. To our knowledge, VDAMP is the first algorithm for variable density Fourier sampling of natural images where a state evolution has been observed.
The joint timefrequency localization provided by the wavelet transform decomposes the color of the effective noise while retaining incoherence. In the numerical examples presented in Section 3 a sparse model on wavelet coefficients is considered, however, we emphasize that VDAMP uses the wavelet transform only as a tool for decomposing the aliasing, and is not constrained to models directly on wavelet coefficients [Metzler2016FromSensing, Xue2016DOAMP:Sensing].
2 Description of algorithm
The VDAMP algorithm is stated in Algorithm 1, where is the entrywise absolute value of a vector or matrix and is entrywise multiplication. The effective noise covariance matrix is represented by a vector that models the diagonal of Eq. 7. The function is a denoiser with a colored effective noise model . The notation in line 7 of Algorithm 1 is defined as the function with th entry
(9) 
where and are the real and imaginary parts respectively. In line 8, the notation is an operator that empirically averages over subbands, so that has the structure
(10) 
with
(11) 
where is the set of indices associated with subband and .
2.1 Weighted gradient descent
To ensure that is an unbiased estimate of , the sensing matrix must be correctly normalized. In VDAMP this is manifest in the gradient step of lines 34, which features a crucial weighting by that is absent in previous applications of AMP to variable density sampling [Sung2013LocationMRI, Eksioglu2018DenoisingBM3DAMPMRI], where a state evolution was not observed. This provides the correct normalisation in expectation over : , which implies that for all . VDAMP’s is the unbiased estimator from Eq. 6. Such a rescaling is referred to as ‘density compensation’ in the MRI literature [Pipe1999SamplingSolution, Pruessmann2001AdvancesTrajectories], and was used in the original compressed sensing MRI paper with zerofilling to generate a unregularized, noniterative baseline [Lustig2007].
Scaling by increases the variance of the measurement noise at frequencies sampled with low probability, and its inclusion in the gradient step will lead to a with higher meansquared error than an unweighted gradient step. However, as shown in Section 3, a careful choice of denoiser that leverages VDAMP’s state evolution can cause lines 69 to be very effective, leading to faster overall convergence than competing methods.
The final step of VDAMP, line 11, is a gradient step without a weighting, which generates a biased image estimate with high data fidelity.
2.2 Colored effective noise model
Line 5 of Algorithm 1 computes an estimate of the colored effective noise covariance matrix from Eq. 7. Through a similar calculation to Eq. 5, the power spectrum of the aliasing of is
(12) 
Using importance sampling with importance distribution , an unbiased estimate of Eq. 12 that does not require knowledge of the ground truth is
(13) 
The computation of in line 5 of Algorithm 1
is a linear transform of
to the wavelet domain. is the power spectrum of , so has unique columns; line 5 therefore requires inner products. We assume the estimator concentrates around its expectation, and leave the study of the constraints this imposes on for future works. Note that for fixed the complexity of VDAMP is governed by and , whose fast implementations have complexity .2.3 Orthogonality between errors
This section motivates the form of lines 69 of Algorithm 1. Consider the function formed by merging these lines:
(14) 
In the following we show that when behaves as Eq. 7, this construction establishes approximate orthogonality between error vectors and . Orthogonality between errors is used in OAMP as a nonrigorous alternative to independence between errors, a requirement for state evolution; see [Ma2017OrthogonalAMP] for details.
Consider a random variable
where and is defined in Eq. 8. Since ,(15) 
which, by Stein’s Lemma [Stein1972AVariables, Eldar_2009], is
(16) 
where the factor of arises due to the covariance matrix of the real and imaginary parts of . If VDAMP’s are instances of the random variable , is an estimate of , and the righthand side of Eq. 16 is approximately zero. Therefore, for large ,
(17) 
Note that this property holds for arbitrary . Two choices for are discussed in Section 3.2.
3 Numerical experiments
This section illustrates the performance of VDAMP compared with FISTA and SUREIT and presents empirical evidence for VDAMP’s state evolution.
3.1 Complex soft threshold tuning with SURE
Selecting appropriate regularisation parameters such as is a notable challenge in realworld compressed sensing applications. We present an approach to parameterfree compressed sensing reconstruction that leverages VDAMP’s state evolution by applying Stein’s Unbiased Risk Estimate (SURE) [Stein1981EstimationDistribution], building on work on AMP in [Mousavi2013ParameterlessPassing, Guo2015NearPassing].
A strength of automatic parameter tuning via SURE is that it is possible to have a richer regularizer than would usually be feasible for a handtuned . In the experiments presented in this section, the denoiser was the complex soft thresholding operator with a subbanddependant threshold that is tuned automatically with SURE. In other words, SURE with a effective noise model given by was used to approximately solve
(18) 
where denotes entrywise division, is the entrywise square root of and has the structure of Eq. 10.
Equation Eq. 18 was solved using a procedure related to SureShrink [Donoho1995AdaptingShrinkage] but for Gaussian noise that is complex and colored [Johnstone1997WaveletNoise]. Consider a vector corrupted by white complex Gaussian noise: . Let be an estimator of . SURE for complex variables has the form
(19) 
where is the empirical averaging operator. cSURE is an unbiased estimate of the risk:
(20) 
Consider the case where is the complex soft thresholding operator with threshold , which acts componentwise as
(21) 
In this case, Eq. 19 is
(22) 
In the following experiments the optimal threshold for each subband was estimated with
(23) 
by evaluating Eq. 22 for trial thresholds . For large dimension
one would expect by the law of large numbers that cSURE is close to the true risk, and for the threshold to be almost optimal. Since a larger number of decomposition scales
give subbands with lower dimensionality, there is a tradeoff between the size of and the quality of threshold selection with cSURE.SURE has previously been employed for parameterfree compressed sensing MRI in [Khare2012AcceleratedParameters], where the Fast Iterative ShrinkingThresholding Algorithm (FISTA) algorithm [Beck2009AProblems] was used with Eq. 18 in place of the usual shrinkage step. This algorithm is herein referred to as SUREIT. The iterationbyiteration error of SUREIT is highly nonGaussian, so deviates from a proper theoretical basis for using SURE and does not give nearoptimal threshold selection.
3.2 Updating
In the experiments in this section two possibilities for the update in line 8 of Algorithm 1 are considered. First,
(24) 
In this case lines 69 of Algorithm 1 are a colored version of the ‘denoising’ phase of VAMP when written in LMMSE form (see [Rangan2019VectorPassing], Algorithm 3). VDAMP with updated with Eq. 24 is herein referred to as VDAMP.
Secondly, as in [Xue2016DOAMP:Sensing], we suggest using SURE for a second time to estimate the that minimises the mean squared error of :
(25) 
Since , optimizing Eq. 19 reduces to a series of minimisation problems with a closedform solution, so that for the th subband equationparentequation
(26a)  
(26b) 
Vector is formed from the so that it has the structure of Eq. 10. VDAMP with updated with is herein referred to as VDAMPS.
3.3 Experimental method
We considered the reconstruction of the SheppLogan shown in Fig. (a)a. The undersampled data was artificially corrupted with complex Gaussian noise so that . We assumed that was known a priori. A Haar wavelet at decomposition scales was used. All experiments were conducted in MATLAB 9.4 and can be reproduced with code available at https://github.com/charlesmillard/VDAMP.
3.4 Comparison with FISTA and SUREIT
To establish comparative reconstruction quality and convergence speed, VDAMP and VDAMPS were compared with FISTA and SUREIT. The scalar sparse weighting of FISTA was tuned with an exhaustive search so that the meansquared error was minimised at . FISTA with a finetuned subbanddependent was not considered as it is not feasible in practise. For SUREIT the meansquared error estimate was updated using the ground truth: , and Eq. 18 with was employed as the denoiser. All algorithms were initialised with a vector of zeros.
Three variable density sampling schemes were considered:

Uniform sampling with .

Twolevel sampling with for a central region and uniform sampling with otherwise, so that .

Polynomial sampling generated using the variable density sampling function from the Sparse MRI package^{1}^{1}1available at https://people.eecs.berkeley.edu/~mlustig/Software.html, with .
A visualization of an instance of in each case is shown in Fig. 3.
Fig. 4 shows the normalized meansquared error (NMSE) versus iteration number for each algorithm for all three sampling sets. VDAMP and VDAMPS were run for iterations, and FISTA and SUREIT were run for 1000 iterations, except in the uniformly sampled case, where FISTA was run for 10000 iterations. All NMSE values were calculated after an additional data consistency step, as in line 11 of Algorithm 1.
In all cases VDAMP and VDAMPS converge within 50 iterations, illustrating the benefit of iterationdependant nearoptimal regularization that is made feasible by SURE. The comparatively poor performance of SUREIT highlights the need for zeromean Gaussian effective noise for effective automatic parameter tuning.
VDAMP and VDAMPS’s reduction in required number of iterations is particularly significant in the uniformly sampled case, Fig. (a)a, where FISTA did not converge even after 10000 iterations. This reflects the highly colored aliasing that arises when the sampling distribution does not compensate for the power spectrum of the signal, and highlights the benefit of modeling the aliasing as colored. Fig. 5 shows , the ground truth Fourier coefficients , and FISTA, SUREIT, VDAMP and VDAMPS reconstructions for uniform sampling at , with entrywise errors .
For twolevel and polynomial sampling the variable density distribution partially compensates for the power spectrum of the signal, so the aliasing is less strongly colored. Figs. (c)c and (b)b show that this leads to a less substantial but still significant improvement in convergence speed for VDAMP and VDAMPS. There is an NMSE improvement in VDAMPS over VDAMP, which indicates that VDAMPS’s update successfully leads to a with lower NMSE than VDAMP. The performance of VDAMP in these figures also highlights that VDAMP can converge to an estimate with a higher NMSE than FISTA, despite the richer regularization. This is due to the scaling in the gradient step, line 4 of Algorithm 1, which effectively increases the measurement noise for coefficients sampled with low probability.
In our implementation the mean periteration compute time over all experiments was 0.15s for FISTA, 0.16s for SUREIT, 0.18s for VDAMP and 0.20s for VDAMPS. The small additional computational cost of SUREIT over FISTA is due to the calculations required for threshold selection with SURE. The further periteration cost of VDAMP is due to the inner products required for the update in line 7 of Algorithm 1, and in the case of VDAMPS, the inner products required for updating by Eq. 26. Despite the additional periteration computational cost, the effectiveness of the denoiser significantly reduces the required computation time. For instance, the time required to reach an NMSE of 35dB in Fig. (a)a was 406s for FISTA, 105s for SUREIT, 3.6s for VDAMP and 3.4s for VDAMPS, and in Fig. (b)b was 12s for FISTA, 11s for SUREIT, 2.5s for VDAMP and 2.0s for VDAMPS.
3.5 Empirical evidence of state evolution
This section presents empirical evidence that VDAMP obeys the colored state evolution given by Eq. 7, using the illustrative case of VDAMPS with polynomial sampling at , as in Fig. (c)c. In Fig. 6, is shown for , indicating that the unbiased subbanddependant effective noise structure shown at in Fig. 2 is preserved. For , Fig. 7 shows quantilequantile plots against a Gaussian of the three illustrative subbands of : the diagonal detail at scale 1, the horizontal detail at scale 2 and the vertical detail at scale 4, where scale 1 is the finest and scale 4 is the coarsest. The linearity of the blue points provides evidence that the effective noise is Gaussian, so that VDAMP has the effective noise given by Eqs. 8 and 7.
4 Conclusions
Based on the observation that Fourier sampling from a nonuniform spectral density leads to colored aliasing, we propose VDAMP, an algorithm based on OAMP that obeys a colored state evolution. State evolution provides an effective way to tune model parameters via SURE, implying that a single algorithm can be used for arbitrary variable density scheme without the need for manual adjustment. More degrees of freedom are feasibly allowed in the model, enabling higher order prior information such as anisotropy, variability across scales and structured sparsity, without the need to estimate the structure a priori such as in modelbased compressed sensing
[Baraniuk2010ModelBasedSensing].VDAMP was motivated by the application of compressed sensing to accelerated MRI. Developments are required for VDAMP to be applicable to MRI data, which is acquired across a number of coils that possess a sensitivity profile[Liang2009AcceleratingSensing, Pruessmann1999SENSE:MRI.]. Further developments are required for VDAMP to be applicable to Fourier sampling with 1D readout curves, where elements of are generated dependently.
It is known that the state evolution of OAMP holds for a wide range of denoisers [Metzler2016FromSensing, Xue2016DOAMP:Sensing]. In [Ahmad2019PlugImaging], a number of standard compressed sensing algorithms that leverage image denoisers designed for Gaussian noise were shown to perform well on MRI reconstruction tasks, despite the mismatch between the aliasing and its model. A sophisticated denoiser equipped to deal with wavelet coefficients corrupted with known colored Gaussian noise would be expected to perform well in conjunction with VDAMP.
Comments
There are no comments yet.