Unbiased Risk Estimates for Singular Value Thresholding and Spectral Estimators

Introduction

Suppose we have noisy observations ${\bm{Y}}$ about an $m\times n$ data matrix ${\bm{X}}^{0}$ of interest,

We wish to estimate ${\bm{X}}^{0}$ as accurately as possible. In this paper, we are concerned with situations where the estimand has some structure, namely, ${\bm{X}}^{0}$ has low rank or is well approximated by a low-rank matrix. This assumption is often met in practice since the columns of ${\bm{X}}^{0}$ can be quite correlated. For instance, these columns may be individual frames in a video sequence, which are typically highly correlated. Another example concerns the acquisition of hyperspectral images in which each column of ${\bm{X}}^{0}$ is a 2D image at a given wavelength. In such settings, images at nearby wavelengths typically exhibit strong correlations. Hence, the special low-rank regression problem (1.1) occurs in very many applications and is the object of numerous recent studies.

Recently, promoting low-rank has been identified as a promising tool for denoising series of MR images, such as those that arise in functional MRI (fMRI) , relaxometry , cardiac MRI , NMR spectroscopy , and diffusion-weighted imaging , among others. In dynamic applications like cine cardiac imaging, where “movies” of the beating heart are created, neighboring tissues tend to exhibit similar motion profiles through time due to their physical connectedness. The diversity of temporal behaviors in this settings will, by nature, be limited, and the Casorati matrix (a matrix whose columns comprise vectorized frames of the image series) formed from this data will be low-rank . For example, in breath-hold cine cardiac imaging, background tissue is essentially static and the large submatrix of the Casoratian corresponding to this region is very well approximated by a rank-1 matrix.

Whenever the object of interest has (approximately) low rank, it is possible to improve upon the naive estimate $\hat{{\bm{X}}}_{0}={\bm{Y}}$ by regularizing the maximum likelihood. A natural approach consists in truncating the singular value decomposition of the observed matrix ${\bm{Y}}$, and solve

where $\lambda$ a positive scalar. As is well known, if

is a singular value decomposition for ${\bm{Y}}$, the solution is given by retaining only the part of the expansion with singular values exceeding $\lambda$,

In other words, one applies a hard-thresholding rule to the singular values of the observed matrix ${\bm{Y}}$. Such an estimator is discontinuous in ${\bm{Y}}$ and a popular alternative approach applies, instead, a soft-thresholding rule to the singular values:

2 A SURE formula

The classical question is, of course, how much shrinkage should be applied. Too much shrinkage results in a large bias while too little results in a high variance. To find the correct trade-off, it would be desirable to have a method that would allow us to compare the quality of estimation for different values of the parameter $\lambda$. Ideally, we would like to select $\lambda$ as to minimize the mean-squared error or risk

Unfortunately, this cannot be achieved since the expectation in (1.6) depends on the true ${\bm{X}}^{0}$, and is thus unknown. Luckily, when the observations follow the model (1.1), it is possible to construct an unbiased estimate of the risk, namely, Stein’s Unbiased Risk Estimate (SURE) given by

when ${\bm{Y}}$ is simple—i.e., has no repeated singular values—and otherwise, say, is a valid expression for the weak divergence. Hence, this simple formula can be used in (1.7), and ultimately leads to the determination of a suitable threshold level by minimizing the estimate of the risk, which only depends upon the observed data.

In MR applications, observations can take on complex values. Indeed, MRI scanners employ a process known as quadrature detection (two orthogonal phase-sensitive detectors) to observe both the magnitude and phase of rotating magnetization that is induced by radio-frequency (RF) excitation . Quadrature detection both increases SNR and allows for the encoding of motion information like flow. MRI noise, which is Gaussian, is also complex valued. In general, MRI noise can also be assumed iid, noting that inter-channel correlation in parallel MRI data can be removed via Cholesky pre-whitening. Thus, model (1.1) has to be modified as

where the real and imaginary parts are also independent. In this case, SURE becomes

We also provide an expression for the weak divergence in this context, namely,

when ${\bm{Y}}$ is simple, and 0 otherwise. This formula can be readily applied to MR data as we shall see in Section 2. The decomposition into real and imaginary parts would suggest that the divergence would be proportional to twice the divergence in the real case. However, this is not the case. The most significant difference is that there is a contribution of the inverse of the singular values even when the matrix is square.

3 Extensions

When the observations ${\bm{Y}}$ follow model (1.1), it is possible to write the degrees of freedom as

Therefore, the expression we provide is also useful to estimate or calculate the degrees of freedom of singular value thresholding.

and that under mild assumptions there exists a closed-form for their divergence:

is spectral. Hence, (1.13) can be used broadly.

4 Connections to other works

During the preparation of this manuscript, the conference paper came to our attention and we would like to point out some differences between this work and ours. In , the authors propose to recursively estimate the divergence of an estimator given by (1.13). This is done by using proximal splitting algorithms. To this end, they provide an expression for the directional subdifferential of a matrix-valued spectral function. The divergence is then estimated by averaging subdifferentials taken along random directions at each iteration of the proximal algorithm. Our approach is obviously different since we provide closed-form formulas for the divergence that are simple to manipulate and easy to evaluate. This has clear numerical and conceptual advantages. In particular, since we have a closed-form expression for SURE, it becomes easier to understand the risk of a family of estimators. Finally, we also address the case of complex-valued data that seems out of the scope of .

5 Content

The structure of the paper is as follows. In Section 2 we present applications in MRI illustrating the advantages of choosing the threshold in a disciplined fashion. Section 3 provides precise statements, and a rigorous justification for (1.8) and (1.10). Section 4 deals with the differentiability of spectral functions (1.11), and thus supports Section 3. We conclude with a short discussion of our results, and potential research directions in Section 5.

Applications in Magnetic Resonance Imaging (MRI)

SVT is a computationally-straightforward, yet, powerful denoising strategy for MRI applications where spatio-temporal/parameteric behavior is either a priori unknown or else defined accordingly to a complicated nonlinear model that may be numerically challenging to work with. A practical challenge in most rank-based estimation problems in MRI (and inverse problems in general), however, lies in the selection of the regularization or threshold parameter. Defining an automated, disciplined, and consistent methodology for selecting this parameter inherently improves clinical workflow both by accelerating the tuning process through the use of optimization-based, rather than heuristic, search strategies and by freeing the MRI scanner technician so that they can focus on other patient-specific tasks. Moreover, eliminating the human element from the denoising process mitigates inter- and intra-operator variability, and thus raises diagnostic confidence in the denoising results since images are wholly reproducible.

Here, we demonstrate a potential use of the SVT unbiased risk estimate for automated and optimized denoising of dynamic cardiac MRI series. Dynamic cardiac imaging is performed either in cine or real-time mode. Cine MRI, the clinical gold-standard for measuring cardiac function/volumetrics , produces a movie of roughly 20 cardiac phases over a single cardiac cycle (heart beat). However, by exploiting the semi-periodic nature of cardiac motion, it is actually formed over many heart beats. Cine sampling is gated to a patient’s heart beat, and as each data measurement is captured it is associated with a particular cardiac phase. This process continues until enough data has been collected such that all image frames are complete. Typically, an entire 2D cine cardiac MRI series is acquired within a single breath hold (less than 30 secs).

In real-time cardiac MRI, an image series covering many heart beats is generated. Although the perceived temporal resolution of this series is coarser than that of a cine series, the “temporal footprint” of each frame is actually shorter since it is only comprised of data from one cardiac cycle. Real-time cardiac MRI is of increasing clinical interest for studying transient functional processes such as first-pass myocardial perfusion, where the hemodynamics of an intravenously-injected contrast bolus are visualized. First-pass myocardial perfusion proffers visualization of both damage to heart muscle as well as coronary artery blockage. Typically, a real-time cardiac MRI study will exceed feasible breath-hold time and respiratory motion may be visible.

Simultaneously achieving high spatial resolution and SNR is challenging in both cine and real-time cardiac MRI due to the dynamic nature of the target signal. Signal averaging (NEX $\geq 1$), a standard MRI techniques for noise reduction, is infeasible in real-time imaging and undesirable in cine imaging due to potential misregistration artifacts, since it cannot be executed within a single breath-hold. Real-time acquisitions, which are generally based on gradient recalled echo (GRE) protocols, have inherently low SNR due to their use of very short repetition (TR) and echo times (TE). Cine acquisitions use either standard GRE or balanced steady-state free precession (bSSFP) protocols. When imaging with 1.5 T (Tesla) magnets, bSSFP cine sequences can often yield sufficient SNR. However, unlike many other MRI protocols, SSFP sequences do not trivially gain SNR when moved to higher-field systems ($\geq$ 3.0 T) which is an emerging clinical trend. As magnetic field strength is increased, the RF excitation flip angle used in a bSSFP sequence must be lowered to adhere to RF power deposition (SAR) safety limits. This results in weaker signal excitation which can mitigate gains in bulk magnetization. Poor receiver coil sensitivity at the heart’s medial location and signal loss due to iron overload (hemochromatosis), among other factors, can further reduce SNR in both imaging strategies. Beyond confounding visual radiological evaluation, even moderate noise levels in a cardiac image can also degrade the performance of automated segmentation methods that are used for quantitative cardiac function evaluation. Therefore, effective denoising techniques that preserve the morphological and dynamic profiles of cardiac image series are clinical valuable.

Consider a series of $t$ separate $n\times{n}$ 2D MR images. To utilize SVT to denoise this data, it must first be transformed into a $n^{2}\times{t}$ Casorati matrix, ${\bm{Y}}={\bm{X}}^{0}+{\bm{W}}$. In many MRI scenarios, spatial dimensionality will greatly exceed temporal/parametric dimensionality ($n^{2}\gg t$) and ${\bm{Y}}$ typically is a very thin matrix. Due to limited degrees-of-freedom, even optimally-parameterized SVT may result in temporal/parametric blurring. One solution to this problem is to analyze the image series in a spatial block-wise manner (see Figure 2) rather than globally .

Let ${\bm{R}}_{b}$ be a binary operator that extracts $k^{2}$ rows from a matrix corresponding to a $k\times{k}$ spatial block, specified by index $b$, within each image. Block-wise SVT can be defined as

where $\Omega$ denotes a set of (potentially overlapping) blocks that uniformly tiles the image domain, i.e., $\sum_{b\in\Omega}{\bm{R}}_{b}^{*}{\bm{R}}_{b}=c{\bm{I}}_{n^{2}\times{n^{2}}}$, $c>0$. In words, (2.1) performs SVT on a family of submatrices of ${\bm{Y}}$ and accumulates a weighted sum of the results. Of course, for $k=n$ and $|\Omega|=1$, (2.1) resorts to standard SVT. The unbiased risk estimator developed in the previous section readily extends for this generalized SVT model. By linearity,

Extending the identity (10) from for matrices then asserts

and the singular value decomposition ${\bm{R}}_{b}{\bm{Y}}={\bm{U}}_{b}{\bm{\Sigma}}_{b}{\bm{V}}_{b}^{*}$. An unbiased estimator of (2.4) is

We now present three MRI examples, one on simulated data and two on real clinical data, and demonstrate the utility of the developed unbiased risk estimators for automatic parameterization of SVT-based denoising.

Initial evaluation of SURE for SVT was performed on the physiologically-improved NCAT (PINCAT) numerical phantom , which simulates a first-pass myocardial perfusion real-time MRI series. In particularly, the free-breathing model ($n=128$, $t=50$) available in the kt-SLR software package was adapted to include a spatially-smooth and temporally-varying phase field such that the target signal was complex valued. Complex iid Gaussian noise ($\tau$ = 30) was then added to the image data. Both standard, or global, and block-wise SVT were each executed at 101 values equispaced over $[10^{-1},10^{7}]$. Block-wise SVT was performed with $k=7$ and $\Omega$ comprising one block for each image pixel, under periodic boundary conditions, such that $|\Omega|=n^{2}$.

Figure 3 shows early, middle, and late time frames from the PINCAT series for the noise-free (truth), noisy, and SVT-denoising results. The threshold values used to generate the SVT results were selected as the MSE/SURE-minimizers in Figure 4a. Also observe in Figure 4a that SURE provides a precise estimate of MSE for both the global and block-wise SVT models. The high accuracy exhibited in this case can be attributed to the high dimensionality of the MRI series data. Note that both global and block-wise SVT yield strong noise reduction generally preserve both morphology and contrast. However, block-wise SVT simultaneously demonstrates a greater degree of noise removal and fidelity to ground truth. In particular, note the relative contrast of the various components of the heart in late frame results. The first observation is corroborated by Figure 4a, which shows that block-wise SVT is able to achieve a lower MSE than global SVT. The second observation is corroborated by Figures 4b-c, which show the worst-case absolute error (compared to ground truth) through time for the two SVT setups. Clearly, global SVT exhibits higher residual error than block-wise SVT, particularly in areas of high motion near the myocardium. The difference between these results can be attributed to the matrix anisotropy problem discussed earlier in this section. Thus, SURE can also be used to automatically-determine the block-size setting as well as the threshold value.

Example 2: Cine Cardiac Imaging

In the second experiment, a bSSFP long-axis cine sequence ($n=192$, $t=19$) acquired at 1.5 T using an phased-array cardiac receiver coil ($\l=8$ channels) was denoised via SVT. The Casorati matrices for individual data channels, which have undergone pre-whitening, are stacked to form a single $ln^{2}\times{t}$ matrix. Assuming that the spatial sensitivity of the receiver channels does not vary substantially through time, the rank of this composite matrix will be equal to that for any individual channel. For visualization purposes, multi-channel denoising results were compressed by calculating the root-sum-of-squares image across the channel dimension. Background noise was determined using a manually-drawn region-of-interest (ROI) to have $\tau$ = 0.67. In this example, block-wise SVT was performed with $k=5$ and used the same block set as in Example 1. The parameter sweep was executed for 101 values equispaced over $[10^{-3},10^{5}]$. In this example, the ground truth is unknown, and thus the MSE cannot be computed to be compared against SURE. However, inspection of Figure 6a reveals that the same qualitative behavior seen for SURE in the numerical phantom example is observed here as well.

Figure 5 demonstrates the effect of sub-optimally selecting the threshold value for both global and block-wise SVT. In particular, over- and under-estimates, along with the SURE-optimal values are applied within SVT. For both SVT strategies, under-estimation of $\lambda$ fails to remove noise as expected. Conversely, over-estimation of $\lambda$ leads to spatio-temporal blurring by both strategies albeit in different manners. In global SVT, temporal blurring occurs near areas of high motion like the tricuspid valve of the heart (indicated via the red arrow). However, less active areas like the pulmonary branching vessels, seen just above the heart, are undistorted. Some background noise also remains visible. In the block-wise SVT result, these characteristics are essentially reversed—minimal temporal blurring is induced but there is noticeable spatial blurring and loss of low-contrast static features. Noise is, however, strongly reduced. A compromise between these extremes is made by the SURE-optimized results. As suggested by Figure 6a, block-wise SVT offers stronger noise suppression (see Figures 6b-c) without inducing spatial or temporal blur, the latter which is seen even in the optimized global SVT result. This example highlights the sensitivity of SVT-denoising performance on parameter selection, and the importance of having a disciplined framework for choosing threshold values. Also of note is that, even after empirical pre-whitening, multi-channel MRI noise may not be exactly iid Gaussian. Nonetheless, SURE allows production of extremely reliable results.

Example 3: First-Pass Myocardial Perfusion

The third denoising experiment was performed on the single-channel, first-pass myocardial perfusion real-time MRI sequence ($n_{x}=190$, $n_{y}=90$, $t=70$) that is also provided with the kt-SLR software package . To simulate a low-field acquisition (1.0 T), such as with an open or interventional MRI system, complex Gaussian noise was added to the data originally acquired at 3.0 T using a GRE sequence. Following this addition, the noise level was estimated at $\tau=0.105$. Since the duration of this exam exceeded feasible breath-hold time, there is substantial respiratory motion at the end of the series. For this example, only block-wise SVT ($k=6$) denoising was performed, and executed at 101 threshold values equispaced over $[10^{-3},10^{5}]$. As in Example 2, $\Omega$ comprises one block for each image pixel, and only SURE can be computed due to the absence of ground truth.

Mirroring Figure 5, Figure 7 demonstrates the effect of threshold selection on different image frames from the perfusion series. In particular, images showing the transition of contrast from the right to left vertical, as well as later-stage onset of myocardial blush are shown. Of the 101 tested threshold values, 5 threshold settings corresponding to 2 over-estimation, 2 under-estimation, and the SURE-optimal value are depicted (see Figure 8). Under all display settings, temporal contrast dynamics are well preserved; however, threshold setting has a marked effect on noise level, spatial resolution, and contrast. As before, employing SVT with an under-estimated threshold fails to remove noise, as evident in the first two columns of Figure 7. Conversely, employing SVT with an over-estimated threshold will remove noise but also induce both spatial blurring and contrast loss. This is particularly evident in the 4th image row ($t=20$). At high threshold values, the contrast of pulmonary vasculature is diminished, and there is visible blurring of the papillary muscles, which appear as dark spots in the contrast-enhanced left ventricle. Finally, the SVT result obtained with the SURE-optimal threshold represents an ideal balance between these extremes, offering strong noise reduction without degradation of important anatomical features. As suggested by Figure 8, there may only exist a narrow parameter window in which effective denoising can be achieved, and the presented unbiased risk estimators can greatly aid in the identification of these optimal settings.

2 Extensions and generalizations

The three examples in the previous subsection demonstrate the utility of the unbiased risk estimation for SVT-based denoising of cardiac MRI series data. Of course, this methodology can also be applied to any other dynamic or parametric MRI series, including those discussed earlier such as functional MRI. Although exhaustive search over a wide range of threshold values was performed for the sake of exposition, the observed unimodality of the risk functional suggests that practical automated parameter selection for SVT denoising could be accelerated using a bisection strategy such as golden section search.

In addition to optimizing over a range of threshold values, unbiased risk estimation can also be used to identify the optimal block-size when using the local SVT model in (2.1). For example, one could imagine performing a sweep of $\lambda$ over a prescribed range of values for a collection of block-size values, and computing the lower envelope of this set of risk measures to investigate best achievable performance as a function of block-size. In addition to simply optimizing a denoising setup, this type of information can be used to guide the development of rank-based reconstruction frameworks for undersampled dynamic and/or parallel MRI.

SURE Formulas for Singular Value Thresholding

Roughly speaking, the derivatives can fail to exist over regions of Lebesgue measure zero.

for $(i,j)\in\mathcal{I}:=\{1,\ldots,m\}\times\{1,\ldots n\}$. Then

Proof Since all norms are equivalent in finite dimensions, $g$ is also Lipschitz with respect to the Frobenius norm. Hence,

Regarding the integrability, since $g({\bm{0}})={\bm{0}}$, (3.2) yields $\|g({\bm{Y}})\|_{F}^{2}\leq L^{2}\|{\bm{Y}}\|_{F}^{2}$ and we deduce

Furthermore, the Cauchy-Schwarz inequality gives

Finally, (3.2) asserts that the derivatives of $g_{ij}$—whenever they exist—are bounded by $L$. Hence,

2 The complex case

In the complex case, we need to analyze the conditions for SURE with respect to model (1.9), which can be expressed in vector form as

We thus need to study the existence of the weak partial derivatives

and whether they satisfy our integrability conditions, as these are the only ones that play a role in the divergence.

As in the real case, we can discard sets of Lebesgue measure zero. Consider then the analogue of $\mathcal{F}$, whose complement has Lebesgue measure zero.

Differentiability of Spectral Functions: General SURE Formulas

Differentiation obeys the standard rules of calculus and we have

Standard analysis arguments immediately establish that the function $f$ is differentiable at ${\bm{\Sigma}}$ whenever ${\bm{\Sigma}}$ is simple, as in this case there is no ambiguity as to which function $f_{i}$ is being applied to the singular values. Furthermore, the differentiability of the SVD of a simple matrix with full-rank is a consequence of Theorem 1 and Theorem 2 in applied to ${\bm{X}}^{*}{\bm{X}}$ and ${\bm{X}}{\bm{X}}^{*}$ (see also and ). We summarize these arguments in the following lemma.

Let ${\bm{X}}={\bm{U}}{\bm{\Sigma}}{\bm{V}}^{*}$ be a simple and full-rank matrix and $f$ be a spectral function.

The factors of the SVD are differentiable in a neighborhood of ${\bm{X}}$.

If $f$ is differentiable at ${\bm{\Sigma}}$, then it is differentiable at ${\bm{X}}$.

We now focus on determining the divergence of $f$ in closed-form. Since there are differences between the real- and complex-valued cases, we treat them separately.

Lemma 4.1 asserts that the SVD is differentiable at any simple matrix ${\bm{X}}$ with full-rank. The following result shows that we can determine the differentials of each factor in closed-form.

Let ${\bm{X}}$ be simple and full-rank. Then the differentials of ${\bm{U}}$ and ${\bm{V}}$ are given by

which says that ${\bm{\Omega}}_{\bm{V}}[{\bm{\Delta}}]$ is anti-symmetric. Similarly,

In particular, (4.5) and (4.6) can be summarized in the system of linear equations

Since ${\bm{X}}$ is simple, the coefficient matrix is invertible, and thus

which is well defined since the orthogonal projector ${\bm{Q}}{\bm{Q}}^{*}$ is independent of the choice of ${\bm{Q}}$.

Since Lemma 4.2 provides a closed-form expression for the differentials, computing the value of the divergence is now a matter of calculus.

By the same arguments, from (4.8) we obtain

Now follow (4.11) and decompose the divergence as

while it follows from (4.9) that the second equals

2 The complex case

Let ${\bm{X}}$ be simple and full-rank. Then the differentials of ${\bm{U}}$ and ${\bm{V}}$ are given by

Furthermore, it is not hard to see that the system (4.8) becomes

and the expression equivalent to (4.9) is

Proof Since ${\bm{X}}$ is simple and has full-rank, Lemma 4.4 holds, and

for $(i,j)\in\mathcal{I}$ and $k=1,\ldots,n$. By the same arguments, (4.15) gives

for $(i,j)\in\mathcal{I}$, $k,l=1,\ldots,n$, and $k\neq l$. Next, (4.16) gives

for $(i,j)\in\mathcal{I}$, $k=n+1,\ldots,m$, and $l=1,\ldots,n$. Finally, (4.14) implies

3 An extension to all matrix arguments

Theorems 4.3 and 4.5 provide the divergence of a matrix-valued spectral function in closed-form for simple, full-rank arguments. Here, we indicate that these formulas have a continuous extension to all matrices. Before continuing, we introduce some notation. For a given ${\bm{X}}$, we denote as $\kappa\leq n$ the number of distinct singular values while $s_{i}$ and $d_{i}$, $i=1,\ldots,\kappa$, denote the $i$-th distinct singular value and its multiplicity.

As previously discussed, when ${\bm{\Sigma}}$ is simple there is no ambiguity as to which $f_{i}$ is being applied to a singular value. However, when ${\bm{\Sigma}}$ is not simple, the ambiguities make $f$ nondifferentiable unless $f_{i}\equiv f$ for $i=1,\ldots,n$.

The first terms are easy to analyze. In fact,

We now formalize the notion of non-tangential limit. Consider a pair $(i,j)$ with $i\neq j$ and $\sigma_{i}=\sigma_{j}:=\sigma$, and define

We say that ${\bm{X}}_{k}$ approaches ${\bm{X}}$ non-tangentially if there exists a constant $c_{ij}>0$ such that $|\rho^{k}_{ij}-1|>c_{ij}$. Roughly speaking, the matrices ${\bm{X}}_{k}$ approach ${\bm{X}}$ fast enough, so that they remain simple, and have full-rank. We have

We distinguish two cases. If $\sigma>0$, we see that

The second sum $I_{1}$ in (4.21) is easier to analyze, as it can be seen that

Discussion

Beyond MRI, SVT denoising also has potential application for medical imaging modalities outside of MRI. For example, in perfusion X-ray computed tomography (CT) studies, where some anatomical region (e.g., the kidneys) is repeatedly scanned to visualize hemodynamics, X-ray tube energy is routinely lowered to limit the ionizing radiation dose delivered to the patient. However, this can result in a substantial increase of image noise. Recent works (e.g., ) have shown that retrospective denoising of low-dose CT data can often yield diagnostic information comparable to full-dose images. It may be possible to adapt SVT, and the proposed risk estimators, for the CT dose reduction problem as well.

E. C. is partially supported by AFOSR under grant FA9550-09-1-0643, by ONR under grant N00014-09-1-0258 and by a gift from the Broadcom Foundation. C. S-L. is partially supported by a Fulbright-CONICYT Scholarship and by the Simons Foundation. J. D. T. and the Mayo Clinic Center for Advanced Imaging Research are partially supported by National Institute of Health grant RR018898. E. C. would like to thank Rahul Mazumder for fruitful discussions about this project.

References