I'm working with non-linear optimization for imaging, such as MRI and CT.

Our problem is of the form $\|Ax-b \|_2^2+\lambda \|Wx\|_1$. $A$ is never formed explicitly, so we're limited to approaches that only use $Ax$ and $A^Hf$. Even so, these evaluations are quite expensive, whereas $Wx$ and $W^Tx$ are quite inexpensive to compute.

I've been using a non-linear conjugate gradient (NCG) for this problem, but my backtracking line search tends to fail as $A^H$ is not the exact conjugate transpose of $A$, but an approximation.

Ignoring the non-linear term and using the conjugate gradient on the normal equations, everything works beautifully, and so does the Split-Bregman algorithm. The latter is however something like an order of magnitude slower than the NCG.

I was wondering if anyone knew if an iterative solver that allowed for the gradient ($A^H(Ax-f)$ to be approximate?

• The complicating factor is having the weights in the $\ell_1$ norm here. $W$ wouldn't be cheap to invert by any chance, would it? That is to say, can you compute $W^{-1}y$ efficiently? Mar 5 '14 at 2:58
• We do have some special cases where $W^T=W^{-1}$, which opens up for stuff like the fast iterative shrinkage (FISTA) method, but typically $W$ calculates the elementwise partial derivative of $x$ and is not invertible. Mar 5 '14 at 9:46
• Well, the smoothed conic dual approach that I co-developed with Becker and Candes in [this paper](Now, with that, you can consider the smoothed conic dual approach that is documented in a paper I co-wrote with Becker and Candes can handle problems like these for exact $A^H$, but I'm not sure how well it would hold up for approximate adjoints. Mar 7 '14 at 3:04

If your gradients are inexact, they are no longer guaranteed to be descent directions, and so any algorithm that uses both functional values and gradients (via line searches or trust regions) and thus relies on that property will run into problems sooner or later.

Splitting methods, on the other hand, allow for inexact evaluation and can be combined with extrapolation and preconditioning to improve the convergence. The classical Chambolle-Pock algorithm (which is a variant of Douglas-Rachford splitting and can be interpreted as a proximal point algorithm) for your problem can be written as

1. $x^k = x^{k-1} - \tau (W^T z^{k-1} + A^T y^{k-1})$
2. $\bar x^k = 2x^k - x^{k-1}$
3. $y^k = \frac{1}{1+\sigma}(y^{k-1}+\sigma (A \bar x^k-b))$
4. $w^k = z^{k-1} + \sigma W\bar x^k$
5. $z^k = \lambda \frac{w^k}{\max\{\lambda,|w^k|\}}$

where the step sizes should satisfy $\sigma\tau < \|A\|^2 +\|W\|^2$. Step 3 and 5 correspond to the proximal mappings of the $l^2$ norm and the $l^1$ norm, respectively.

Preconditioners are discussed in this preprint by Bredies and Sun, where $K = (A,W)^T$. The case of inexact solution of the normal equations (which is not quite what you are asking about, but allows for errors to remain bounded from below) is discussed after Remark 2.2.

In general, proximal point methods converge as long as the evaluation error of the proximal mapping remains summable (i.e., evaluation has to become exact sufficiently quickly as the iteration progresses), see for example http://arxiv.org/pdf/1109.2415v2.pdf and the references therein.