# Preconditioned ADMM with Nonlinear Operator Constraint

Florian Knoll, Carola-Bibiane Schönlieb, Tuomo Valkonen, Martin Benning

► To cite this version:

Florian Knoll, Carola-Bibiane Schönlieb, Tuomo Valkonen, Martin Benning. Preconditioned ADMM with Nonlinear Operator Constraint. 27th IFIP Conference on System Modeling and Optimization (CSMO), Jun 2015, Sophia Antipolis, France. pp.117-126, 10.1007/978-3-319-55795-3_10. hal-01626884

HAL Id: hal-01626884

https://inria.hal.science/hal-01626884

Submitted on 31 Oct 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Preconditioned ADMM with nonlinear
operator constraint

Martin Benning<fnref n="1" />, Florian Knoll<fnref n="2" />, Carola-Bibiane Schönlieb<fnref n="1" />, and
Tuomo Valkonen<fnref n="3" />

**Abstract.** We are presenting a modification of the well-known Alternat-
ing Direction Method of Multipliers (ADMM) algorithm with additional
preconditioning that aims at solving convex optimisation problems with
nonlinear operator constraints. Connections to the recently developed
Nonlinear Primal-Dual Hybrid Gradient Method (NL-PDHGM) are pre-
sented, and the algorithm is demonstrated to handle the nonlinear inverse
problem of parallel Magnetic Resonance Imaging (MRI).

**Keywords:** ADMM, primal-dual, nonlinear inverse problems, parallel MRI, proximal point method, operator splitting, iterative Bregman method

# 1 Introduction

Non-smooth regularisation methods are popular tools in the imaging sciences. They allow to promote sparsity of inverse problem solutions with respect to specific representations; they can implicitly restrict the null-space of the forward operator while guaranteeing noise suppression at the same time. The most prominent representatives of this class are total variation regularisation [19] and ℓ<fnref n="1" />-norm regularisation as in the broader context of compressed sensing [10, 8].

In order to solve convex, non-smooth regularisation methods with linear op-
erator constraints computationally, first-order operator splitting methods have
gained increasing interest over the last decade, see [11, 12, 3, 9] to name just a
few. Despite some recent extensions to certain types of non-convex problems [16,
14, 7, 15] there has to our knowledge only been made little progress for nonlinear
operators constraints [2, 22].

In this paper we are particularly interested in minimising non-smooth, con-
vex functionals with nonlinear operator constraints. This model covers many
interesting applications; one particular application that we are going to address
is the joint reconstruction of the spin-proton-density and coil sensitivity maps
in parallel MRI [21, 13].

<footnote n="1">¹ University of Cambridge, Department of Applied Mathematics and Theoretical Physics, Wilberforce Road, Cambridge CB3 0WA, United Kingdom
({mb941, cbs31}@cam.ac.uk)</footnote>

<footnote n="2">² New York University, Center for Advanced Imaging, Innovation and Research, New York 4th Floor 660 First Avenue New York, NY 10016, USA
(florian.knoll@nyumc.org)</footnote>

<footnote n="3">³ University of Liverpool, Department of Mathematical Sciences, Mathematical Sciences Building, L69 7ZL Liverpool, United Kingdom
(tuomo.valkonen@liverpool.ac.uk)</footnote>
The paper is structured as follows: we will introduce the generic problem formulation, then address its numerical minimisation via a generalised ADMM method with linearised operator constraints. Subsequently we will show connections to the recently proposed NL-PDHGM method (indicating a local convergence result of the proposed algorithm) and conclude with the joint spin-proton-density and coil sensitivity map estimation as a numerical example.

## 2 Problem formulation

We consider the following generic constrained minimisation problem:

$$ (\hat{u}, \hat{v}) = \arg \min_{u,v} \{H(u) + J(v) \text{ subject to } F(u, v) = c\}. \quad (1) $$

Here $H$ and $J$ denote proper, convex and lower semi-continuous functionals, $F$ is a nonlinear operator and $c$ a given function. Note that for nonlinear operators of the form $F(u, v) = G(u) - v$ and $c = 0$ problem (1) can be written as

$$ \hat{u} = \arg \min_{u} \{H(u) + J(G(u))\}. \quad (2) $$

In the following we want to propose a strategy for solving (1) that is based on simultaneous linearisation of the nonlinear operator constraint and the solution of an inexact ADMM problem.

## 3 Alternating direction method of multipliers

We solve (1) by alternating optimisation of the augmented Lagrange function

$$ \mathcal{L}_{\delta}(u, v; \mu) = H(u) + J(v) + \langle \mu, F(u, v) - c \rangle + \frac{\delta}{2} \|F(u, v) - c\|_{2}^{2}. \quad (3) $$

Alternating minimisation of (3) in $u, v$ and subsequent maximisation of $\mu$ via a step of gradient ascent yields this nonlinear version of ADMM [11]:

$$ u^{k+1} \in \arg \min_{u} \left\{ \frac{\delta}{2} \|F(u, v^k) - c\|_2^2 + \langle \mu^k, F(u, v^k) \rangle + H(u) \right\}, \quad (4) $$

$$ v^{k+1} \in \arg \min_{v} \left\{ \frac{\delta}{2} \|F(u^{k+1}, v) - c\|_{2}^{2} + \langle \mu^{k}, F(u^{k+1}, v) \rangle + J(v) \right\}, \quad (5) $$

$$ \mu^{k+1} = \mu^k + \delta(F(u^{k+1}, v^{k+1}) - c). \quad (6) $$

Not having to deal with nonlinear subproblems, we replace $F(u, v^k)$ and $F(u^{k+1}, v)$ by their Taylor linearisations around $u^k$ and $v^k$, which yields $F(u, v^k) \approx F(u^k, v^k) + \partial_u F(u^k, v^k)(u - u^k)$ and $F(u^{k+1}, v) \approx F(u^{k+1}, v^k) + \partial_v F(u^{k+1}, v^k)(v - v^k)$, respectively. The updates (4) and (5) modify to

$$ u^{k+1} \in \arg \min_{u} \left\{ \frac{\delta}{2} \|A^k u - c_1^k\|_2^2 + \langle \mu^k, A^k u \rangle + H(u) \right\}, \quad (7) $$

$$ v^{k+1} \in \arg \min_{v} \left\{ \frac{\delta}{2} \|B^k v - c_2^k\|_2^2 + \langle \mu^k, B^k v \rangle + J(v) \right\}, \quad (8) $$
with $A^k := \partial_u F(u^k, v^k)$, $B^k := \partial_v F(u^{k+1}, v^k)$, $c_1^k := c + A^k u^k - F(u^k, v^k)$ and $c_2^k := c + B^k v^k - F(u^{k+1}, v^k)$. Note that the updates (7) and (8) are still implicit, regardless of $H$ and $J$. In the following, we want to modify the updates such that they become simple proximity operations.

## 4 Preconditioned ADMM

Based on [23], we modify (7) and (8) by adding the surrogate terms $\|u^{k+1} - u^k\|_{Q_1^k}^2/2$ and $\|v^{k+1} - v^k\|_{Q_2^k}^2/2$, with $\|w\|_Q := \sqrt{\langle Qw, w \rangle}$ (note that if $Q$ is chosen to be positive definite, $\| \cdot \|_Q$ becomes a norm). We then obtain

$$ 
\begin{aligned}
u^{k+1} &\in \arg\min_u \left\{ \frac{\delta}{2} \|A^k u - c_1^k\|_2^2 + \langle \mu^k, A^k u \rangle + H(u) + \frac{1}{2} \|u - u^k\|_{Q_1^k}^2 \right\}, \\
v^{k+1} &\in \arg\min_v \left\{ \frac{\delta}{2} \|B^k v - c_2^k\|_2^2 + \langle \mu^k, B^k v \rangle + J(v) + \frac{1}{2} \|v - v^k\|_{Q_2^k}^2 \right\}.
\end{aligned}
$$

If we choose $Q_1^k := \tau_1^k I - \delta A^{k*} A^k$ with $\tau_1^k \delta < 1/\|A^k\|^2$ and $Q_2^k := \tau_2^k I - \delta B^{k*} B^k$ with $\tau_2^k \delta < 1/\|B^k\|^2$ and if we define $\bar{\mu}^k := 2\mu^k - \mu^{k-1}$ we obtain

$$ u^{k+1} = (I + \tau_1^k \partial H)^{-1} (u^k - \tau_1^k A^{k*} \bar{\mu}^k), \quad (9) $$

$$ v^{k+1} = (I + \tau_2^k \partial J)^{-1} (v^k - \tau_2^k B^{k*} (\mu^k + \delta(F(u^{k+1}, v^k) - c))), \quad (10) $$

with $(I + \alpha \partial E)^{-1}(w)$ denoting the proximity or resolvent operator

$$ (I + \alpha \partial E)^{-1}(w) := \arg\min_u \left\{ \frac{1}{2} \|u - w\|_2^2 + \alpha E(u) \right\}. $$

The entire proposed algorithm with updates (9), (10) and (6) reads as

**Algorithm 1 Preconditioned ADMM with nonlinear operator constraint**

**Parameters:** $H, J, F, c$

**Initialization:** $u^0, v^0, \mu^0, \delta$

$$ \bar{\mu}^0 = \mu^0 $$

**while** convergence criterion is not met **do**

$$ A^k = \partial_u F(u^k, v^k) $$

Set $\tau_1^k$ such that $\tau_1^k \delta < 1/\|A^k\|^2$

$$ u^{k+1} = (I + \tau_1^k \partial H)^{-1} (u^k - \tau_1^k A^{k*} \bar{\mu}^k) $$

$$ B^k = \partial_v F(u^{k+1}, v^k) $$

Set $\tau_2^k$ such that $\tau_2^k \delta < 1/\|B^k\|^2$

$$ v^{k+1} = (I + \tau_2^k \partial J)^{-1} (v^k - \tau_2^k B^{k*} (\mu^k + \delta(F(u^{k+1}, v^k) - c))) $$

$$ \mu^{k+1} = \mu^k + \delta(F(u^{k+1}, v^{k+1}) - c) $$

$$ \bar{\mu}^{k+1} = 2\mu^{k+1} - \mu^k $$

**end while**

**return** $u^k, v^k, \mu^k, \bar{\mu}^k$
## 5 Connection to NL-PDHGM

In the following we want to show how the algorithm simplifies in case the nonlinear operator constraint is only nonlinear in one variable, which is sufficient for problems of the form (2). Without loss of generality we consider constraints of the form $F(u, v) = G(u) - v$, where $G$ represents a nonlinear operator in $u$. Then we have $A^k = \mathcal{J}G(u^k)$ (with $\mathcal{J}G(u^k)$ denoting the Jacobi matrix of $G$ at $u^k$), $B^k = -I$ and if we further choose $\tau_2^k = 1/\delta$ for all $k$, update (10) reads

$$v^{k+1} = \left(I + \frac{1}{\delta}\partial J\right)^{-1} \left(G(u^{k+1}) + \frac{1}{\delta}\mu^k\right).$$

Applying Moreau's identity [18] $b = (I + \frac{1}{\delta}\partial J)^{-1}(b) + \frac{1}{\delta}(I + \delta\partial J^*)^{-1}(\delta b)$ yields

$$\mu^{k+1} = (I + \delta\partial J^*)^{-1} (\mu^k + \delta G(u^{k+1})).$$

If we further change the order of the updates, starting with the update for $\mu$, the whole algorithm reads

$$\mu^{k+1} = (I + \delta\partial J^*)^{-1} (\mu^k + \delta G(u^k)),$$

$$\bar{\mu}^{k+1} = 2\mu^{k+1} - \mu^k,$$

$$u^{k+1} = (I + \tau_1^k \partial H)^{-1} (u^k - \tau_1^k \mathcal{J} G(u^k)^* \bar{\mu}^{k+1}).$$

Note that this algorithm is almost the same as NL-PDHGM proposed in [22] for $\theta = 1$, except that the extrapolation step is carried out on the dual variable $\mu$ instead of the primal variable $u$. In the following we want to briefly sketch how to prove convergence for this algorithm in analogy to [22]. We define

$$N(\mu^{k+1}, u^{k+1}) := \left( \frac{\partial J^*(\mu^{k+1}) - \nabla G(u^k)u^{k+1} - c^k}{\partial H(u^{k+1}) + \mathcal{J}G(u^k)^*\mu^{k+1}} \right),$$

$$L^k := \begin{pmatrix} \frac{1}{\delta}I & \mathcal{J}G(u^k) \\ \mathcal{J}G(u^k)^* & \frac{1}{\tau_1^k}I \end{pmatrix},$$

with $c^k := G(u^k) - \mathcal{J}G(u^k)u^k$. Now the algorithm is: find $(\mu^{k+1}, u^{k+1})$ such that

$$N(\mu^{k+1}, u^{k+1}) + L^k(\mu^{k+1} - \mu^k, u^{k+1} - u^k) \ni 0.$$

If we exchange the order of $\mu$ and $u$ here, i.e., reorder the rows of $N$, and the rows and columns of $L^k$, we obtain almost the “linearised” NL-PDHGM of [22]. The difference is that the sign of $\mathcal{J}G$ in $L^k$ is inverted. The only points in [22] where the exact structure of $L^k$ ($M_{xk}$ therein) is used, are Lemma 3.1, Lemma 3.6 and Lemma 3.10. The first two go through exactly as before with the negated structure. Reproducing Lemma 3.10 demands bounding actual step lengths $\|u^k - u^{k+1}\|$ and $\|\mu^k - \mu^{k+1}\|$ from below, near a solution for arbitrary $\epsilon > 0$. A proof would go beyond the page limit of this proceeding. Let us just point out that this can be done, implying that the convergence results of [22] apply for this algorithm as well. This means that under somewhat technical conditions, which for TV type problems amount to Huber regularisation, local convergence in a neighbourhood of the true solution can be guaranteed.
6 Joint estimation of the spin-proton density and coil sensitivities in parallel MRI

We want to demonstrate the numerical capabilities of Algorithm 1 by applying it to the nonlinear problem of joint estimation of the spin-proton density and the coil sensitivities in parallel MRI. The discrete problem of joint reconstruction from sub-sampled k-space data on a rectangular grid reads

$$ \begin{pmatrix} \hat{u} \\ \hat{c}_1 \\ \vdots \\ \hat{c}_2 \end{pmatrix} \in \arg \min_{\mathbf{v}=(u,c_1,\dots,c_n)} \left\{ \frac{1}{2} \sum_{j=1}^{n} \| S \mathcal{F}(G(\mathbf{v}))_j - f_j \|_2^2 + \alpha_0 R_0(u) + \sum_{j=1}^{n} \alpha_j R_j(c_j) \right\}, $$

where $\mathcal{F}$ is the 2D discrete Fourier transform, $f_j$ are the k-space measurements for each of the $n$ coils, $S$ is the sub-sampling operator and $R_j$ denote appropriate regularisation functionals. The nonlinear operator $G$ maps the unknown spin-proton density $u$ and the different coil sensitivities $c_j$ as follows [21]:

$$ G(u, c_1, \dots, c_n) = (uc_1, uc_2, \dots, uc_n)^T. \qquad (11) $$

In order to compensate for sub-sampling artefacts in sub-sampled MRI it is common practice to use total variation as a regulariser [6, 17]. Coil sensitivities are assumed to be smooth, cf. Figure 1, motivating a reconstruction model similar to the one proposed in [13]. We therefore choose the discrete isotropic total variation, $R_0(u) = \|\nabla u\|_{2,1}$, and the smooth 2-norm of the discretised gradient, i.e. $R_j(c_j) := \|\nabla c_j\|_{2,2}$, for all $j > 0$, following the notation in [4]. We further introduce regularisation parameters $\lambda_j$ in front of the data fidelities and rescale all regularisation parameters such that $\alpha_0 + \frac{1}{n} (\sum_{j=1}^n \lambda_j + \sum_{j=1}^n \alpha_j) = 1$. In order to realise this model via Algorithm 1 we consider the following operator splitting strategy. We define $F(u_0, \dots, u_n, v_0, \dots, v_{2n})$ as

$$ F(u_0, \dots, u_n, v_1, \dots, v_n) := \begin{pmatrix} G(u_0, \dots, u_n) \\ \nabla u_0 & 0 & \dots & 0 \\ 0 & \nabla u_1 & \dots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \dots & 0 & \nabla u_n \end{pmatrix} - \begin{pmatrix} v_0 \\ \vdots \\ v_n \\ \vdots \\ v_{2n} \end{pmatrix}, $$

set $H(u_0, \dots, u_n) \equiv 0$, and $J(v_0, \dots, v_{2n}) = \sum_{j=0}^{2n} J_j(v_j)$ with $J_j(v_j) := \frac{\lambda_j}{2} \| S \mathcal{F} v_j - f_j \|_2^2$ for $j \in \{0, \dots, n-1\}$, $J_n(v_n) = \alpha_0 \| v_n \|_{2,1}$ and $J_j(v_j) = \alpha_{j-n} \| v_j \|_{2,2}$ for $j \in \{n+1, \dots, 2n\}$. Note that with these choices of functions, all the resolvent
operations can be carried out easily. In particular, we obtain

$$(I + \tau_1^k \partial H)^{-1}(w) = w,$$

$$(I + \tau_2^k \partial J_j)^{-1}(w) = \mathcal{F}^{-1} \left( \frac{\mathcal{F} w_j + \tau_2^k \lambda_j S^T f_j}{1 + \tau_2^k \lambda_j \mathrm{diag}(S^T S)} \right) \quad \text{for } j \in \{0, \dots, n-1\},$$

$$(I + \tau_2^k \partial J_n)^{-1}(w) = \frac{w_n}{\|w_n\|_2} \max (\|w_n\|_2 - \alpha_0 \tau_2^k, 0),$$

$$(I + \tau_2^k \partial J_j)^{-1}(w) = \frac{w_j}{\|w_j\|_{2,2}} \max (\|w_j\|_{2,2} - \alpha_{j-n} \tau_2^k, 0) \quad \text{for } j \in \{n+1, \dots, 2n\}.$$

Moreover, as $B_k = -I$ (and thus, $\|B^k\| = 1$) for all $k$, we can simply eliminate $\tau_2^k$ by replacing it with $1/\delta$, similar to Section 5.

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 1: Figure 1a shows the brain phantom as described in Section 6.1. Figure 1c - 1j show visualisations of the measured coil sensitivities of a water bottle. Figure 1b shows the simulated, spiral-shaped sub-sampling scheme used to sub-sample the k-space data.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 2: Reconstructions for noise with low noise level $\sigma = 0.05$. Despite the sub-sampling, features of the brain phantom are very well preserved. In addition, the coil sensitivities seem to correspond well to the original ones, despite a slight loss of contrast. Note that coil sensitivities remain the initial value where the signal is zero.</figcaption></figure>

## 6.1 Experimental setup

We now want to discuss the experimental setup. We want to reconstruct the synthetic brain phantom in Figure 1a from sub-sampled k-space measurements. The numerical phantom is based on the design in [1] with a matrix size of 190 × 190. It consists of several different tissue types like cerebrospinal fluid (CSF), gray matter (GM), white matter (WM) and cortical bone. Each pixel is assigned a set of MR tissue properties: Relaxation times $T_1(x, y)$ and $T_2(x, y)$ and spin density $\rho(x, y)$. These parameters were also selected according to [1]. The MR signal $s(x, y)$ in each pixel was then calculated by using the signal equation of a fluid attenuation inversion recovery (FLAIR) sequence [5]:

$$s(x, y) = \rho(x, y)(1 - 2 e^{-TI/T_1(x,y)})(1 - e^{-TR/T_1(x,y)}) e^{-TE/T_2(x,y)}$$

The sequence parameters were selected: TR = 10000 ms, TE = 90 ms. TI was set to 1781 ms to achieve signal nulling of CSF ($T_1^{csf} \log(2)$ with $T_1^{csf} = 2569$ms).

In order to generate artificial k-space measurements for each coil, we proceed as follows. First, we produce 8 images of the brain phantom multiplied by the measured coil sensitivity maps shown in Figure 1c - 1j. The coil sensitivity maps were generated from the measurements of a water bottle with an 8-channel head coil array. Then we produce artificial k-space data by applying the 2D discrete Fourier-transform to each of those individual images. Subsequently, we
<figure><img src="image_5.png" /><figcaption>Fig. 3: Reconstructions for noise with high noise level $\sigma = 0.95$. Due to the large amount of noise, higher regularisation parameters are necessary. As a consequence, fine structures are smoothed out and in contrast to the case of little noise, compensation of sub-sampling artefacts is less successful.</figcaption></figure>

sub-sample only approx. 25% of each of the k-space datasets via the spiral shown in Figure 1b. Finally, we add Gaussian noise with standard deviation $\sigma$ to the sub-sampled data.

## 6.2 Computations

For the actual computations we use two noisy versions $f_j$ of the simulated k-space data; one with small noise ($\sigma = 0.05$) and one with a high amount of noise ($\sigma = 0.95$). As stopping criterion we simply choose a fixed number of iterations; for both the low noise level as well as the high noise level dataset we have fixed the number of iterations to 1500. The initial values used for the algorithm are $u_j^0 = \mathbf{1}$ with $\mathbf{1} \in \mathbb{R}^{l \times 1}$ being the constant one-vector, for all $j \in \{0, \dots, n\}$. All other initial variables ($v^0, \mu^0, \bar{\mu}^0$) are set to zero.

**Low noise level** We have computed reconstructions from the noisy data with noise level $\sigma = 0.05$ via Algorithm 1, with regularisation parameters set to $\lambda_j = 0.0621$, $\alpha_0 = 0.0062$ and $\alpha_j = 0.9317$ for $j \in \{1, \dots, n\}$. We have further created a naïve reconstruction by averaging the individual inverse Fourier-transformed images obtained from zero-filling the k-space data. The modulus images of the results are visualised in Figure 2. The PSNR values for the averaged zero-filled reconstruction is 10.2185, whereas the PSNR of the reconstruction with the proposed method is 24.5572.
**High noise level** We proceeded as in the previous section, but for noisy data with noise level $\sigma = 0.95$. The regularisation parameters were set to $\lambda_j = 0.0149$, $\alpha_0 = 0.0135$ and $\alpha_j = 0.9716$ for $j \in \{1, ..., n\}$. The modulus images of the results are visualised in Figure 3. The PSNR values for the averaged zero-filled reconstruction is 9.9621, whereas the PSNR of the reconstruction with the proposed method is 16.672.

# 7 Conclusions & outlook

We have presented a novel algorithm that allows to compute minimisers of a sum of convex functionals with nonlinear operator constraint. We have shown the connection to the recently proposed NL-PDHGM algorithm which implies local convergence results in analogy to those derived in [22]. Subsequently we have demonstrated the computational capabilities of the algorithm by applying it to a nonlinear joint reconstruction problem in parallel MRI.

For future work, the convergence of the algorithm in the general setting has to be verified, and possible extensions to guarantee global convergence have to be studied. Generalisation of stopping criteria such as a linearised primal-dual gap will also be of interest as well. With respect to the presented parallel MRI application, exact conditions for the convergence (like the exact norm of the bounds) have to be verified. The impact of the algorithm- as well as the regularisation-parameters on the reconstruction has to be analysed, and a rigorous study with artificial and real data would also be desirable. Moreover, future research will focus on alternative regularisation functions, e.g. based on spherical harmonics motivated by [20]. Last but not least, other applications that can be modelled via (1) should be considered in future research.

**Acknowledgments** MB, CS and TV acknowledge EPSRC grant EP/M00483X/1. FK acknowledges National Institutes of Health grant NIH P41 EB017183.

**EPSRC data statement:** the corresponding code and data are available for download at https://www repository.cam.ac.uk/handle/1810/256221

# References

1. Berengere Aubert-Broche, Alan C Evans, and Louis Collins. A new improved version of the realistic digital brain phantom. *NeuroImage*, 32(1):138-145, 2006.

2. Markus Bachmayr and Martin Burger. Iterative total variation schemes for nonlinear inverse problems. *Inverse Problems*, 25(10):105004, 2009.

3. Amir Beck and Marc Teboulle. Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems. *Image Processing, IEEE Transactions on*, 18(11):2419-2434, 2009.

4. Martin Benning, Lynn Gladden, Daniel Holland, Carola-Bibiane Schönlieb, and Tuomo Valkonen. Phase reconstruction from velocity-encoded mri measurements—a survey of sparsity-promoting variational approaches. *Journal of Magnetic Resonance*, 238:26-43, 2014.
5. Matt A Bernstein, Kevin F King, and Xiaohong Joe Zhou. *Handbook of MRI pulse sequences*. Elsevier, 2004.

6. Kai Tobias Block, Martin Uecker, and Jens Frahm. Undersampled radial mri with multiple coils. iterative image reconstruction using a total variation constraint. *Magnetic resonance in medicine*, 57(6):1086–1098, 2007.

7. S. Bonettini, I. Loris, F. Porta, and M. Prato. Variable metric inexact line-search based methods for nonsmooth optimization. *Siam Journal on Optimization*, 2015. Submitted.

8. Emmanuel J Candes et al. Compressive sampling. In *Proceedings of the international congress of mathematicians*, volume 3, pages 1433–1452. Madrid, Spain, 2006.

9. Antonin Chambolle and Thomas Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. *Journal of Mathematical Imaging and Vision*, 40(1):120–145, 2011.

10. David L Donoho. Compressed sensing. *Information Theory, IEEE Transactions on*, 52(4):1289–1306, 2006.

11. Daniel Gabay. Applications of the method of multipliers to variational inequalities. *Studies in mathematics and its applications*, 15:299–331, 1983.

12. Tom Goldstein and Stanley Osher. The split bregman method for 11-regularized problems. *SIAM Journal on Imaging Sciences*, 2(2):323–343, 2009.

13. Florian Knoll, Christian Clason, Kristian Bredies, Martin Uecker, and Rudolf Stollberger. Parallel imaging with nonlinear reconstruction using variational penalties. *Magnetic Resonance in Medicine*, 67(1):34–41, 2012.

14. Thomas Möllenhoff, Evgeny Strekalovskiy, Michael Möller, and Daniel Cremers. The primal-dual hybrid gradient method for semiconvex splittings. *SIAM Journal on Imaging Sciences*, 8(2):827–857, 2015.

15. Michael Möller, Martin Benning, Carola Schönlieb, and Daniel Cremers. Variational depth from focus reconstruction. *Image Processing, IEEE Transactions on*, 24(12):5369–5378, Dec 2015.

16. Peter Ochs, Yunjin Chen, Thomas Brox, and Thomas Pock. ipiano: inertial proximal algorithm for nonconvex optimization. *SIAM Journal on Imaging Sciences*, 7(2):1388–1419, 2014.

17. Sathish Ramani, Jeffrey Fessler, et al. Parallel mr image reconstruction using augmented lagrangian methods. *Medical imaging, IEEE Transactions on*, 30(3):694–706, 2011.

18. R Tyrrell Rockafellar. Convex analysis (princeton mathematical series). Princeton University Press, 46:49, 1970.

19. Leonid I Rudin, Stanley Osher, and Emad Fatemi. Nonlinear total variation based noise removal algorithms. *Physica D: Nonlinear Phenomena*, 60(1):259–268, 1992.

20. Alessandro Sbrizzi, Hans Hoogduin, Jan J Lagendijk, Peter Luijten, and Cornelis AT den Berg. Robust reconstruction of b1+ maps by projection into a spherical functions space. *Magnetic Resonance in Medicine*, 71(1):394–401, 2014.

21. Martin Uecker, Thorsten Hohage, Kai Tobias Block, and Jens Frahm. Image reconstruction by regularized nonlinear inversion joint estimation of coil sensitivities and image content. *Magnetic Resonance in Medicine*, 60(3):674–682, 2008.

22. Tuomo Valkonen. A primal-dual hybrid gradient method for nonlinear operators with applications to mri. *Inverse Problems*, 30(5):055012, 2014.

23. Xiaoqun Zhang, Martin Burger, and Stanley Osher. A unified primal-dual algorithm framework based on bregman iteration. *Journal of Scientific Computing*, 46(1):20–46, 2011.