PaperlayerView exact revision on arXiv

Source-linked AI summary

An Augmented Lagrangian Approach to the Constrained Optimization Formulation of Imaging Inverse Problems

Manya V. Afonso, José M. Bioucas-Dias, Mário A. T. Figueiredo

arXiv:0912.3481v1math.OCmath.NA

TL;DR

Ill-posed imaging inverse problems combine huge dimensionality with non-smooth regularization, limiting direct use of general optimization tools. The paper introduces C-SALSA, which transforms the constrained problem through variable splitting and applies ADMM. Across image recovery benchmarks, C-SALSA is usually faster than prior state-of-the-art methods while supporting several convex regularizers.

  • Problem

    Huge dimensions and non-smooth convex regularizers make constrained image-recovery problems difficult for off-the-shelf optimization methods.

  • Method

    C-SALSA transforms the constrained problem using an indicator function and variable splitting, then applies ADMM and efficiently computed Moreau proximity operators.

  • Results

    C-SALSA is usually faster than previous state-of-the-art methods on standard deconvolution, MRI reconstruction, and inpainting problems.

  • Takeaways & Limitations

    The method supports convex wavelet-based, total-variation, and more general regularizers across several imaging inverse problems.

Abstract

from arXiv · show

We propose a new fast algorithm for solving one of the standard approaches to ill-posed linear inverse problems (IPLIP), where a (possibly non-smooth) regularizer is minimized under the constraint that the solution explains the observations sufficiently well. Although the regularizer and constraint are usually convex, several particular features of these problems (huge dimensionality, non-smoothness) preclude the use of off-the-shelf optimization tools and have stimulated a considerable amount of research. In this paper, we propose a new efficient algorithm to handle one class of constrained problems (often known as basis pursuit denoising) tailored to image recovery applications. The proposed algorithm, which belongs to the family of augmented Lagrangian methods, can be used to deal with a variety of imaging IPLIP, including deconvolution and reconstruction from compressive observations (such as MRI), using either total-variation or wavelet-based (or, more generally, frame-based) regularization. The proposed algorithm is an instance of the so-called "alternating direction method of multipliers", for which convergence sufficient conditions are known; we show that these conditions are satisfied by the proposed algorithm. Experiments on a set of image restoration and reconstruction benchmark problems show that the proposed algorithm is a strong contender for the state-of-the-art.

I. INTRODUCTION

Imaging inverse problems estimate an image from noisy indirect observations despite ill-conditioning, using constrained regularization with convex but often non-smooth penalties. The paper considers synthesis and analysis formulations based on wavelet-like frames and total variation.

  • Ill-posed inverse problems require regularization because the forward operator may be singular or very ill-conditioned.
  • The constrained formulation minimizes a regularizer subject to an observation-error bound determined by the noise level.
  • With an ℓ1 regularizer, this formulation is known as basis pursuit denoising, while ε = 0 yields basis pursuit.
  • Typical non-smooth regularizers include total variation and the ℓ1 norm.
  • Wavelet-like representations may use orthogonal bases or overcomplete dictionaries, and analysis operators are assumed to arise from 1-tight frames.
  • Synthesis approaches estimate frame coefficients β and reconstruct x = Wβ, whereas analysis approaches regularize the image directly, commonly with total variation.

C. Previous Algorithms

Earlier methods address constrained or unconstrained inverse formulations with specialized first-order, shrinkage, Bregman, and splitting techniques. Their applicability varies with the regularizer, constraint parameter, operator representation, and problem scale, motivating C-SALSA.

  • C. Previous Algorithms: Very high dimensions and operator-only access to B, W, or P prevent direct use of many off-the-shelf optimization algorithms.Matrix-vector products remain efficiently computable through transforms such as FFTs and wavelet algorithms.
  • C. Previous Algorithms: Shrinkage and thresholding methods rely on the Moreau proximal mapping, which becomes soft thresholding for the ℓ1 norm.
  • C. Previous Algorithms: SALSA uses variable splitting and augmented Lagrangians for unconstrained problems and was reported faster than FISTA, TwIST, and SpaRSA on several benchmarks.
  • C. Previous Algorithms: The constrained formulation has a noise-linked parameter ε, whereas the equivalent unconstrained formulation uses τ, whose correct value may require solving multiple instances.
  • C. Previous Algorithms: NESTA applies smoothed accelerated first-order methods and supports ℓ1 or total-variation regularization in synthesis or analysis formulations.
  • C. Previous Algorithms: C-SALSA transforms the constrained problem through an indicator function and variable splitting before applying ADMM, supporting convex regularizers with efficiently computable proximal mappings.
  • C. Previous Algorithms: C-SALSA was experimentally reported faster than SPGL1 and NESTA on MRI reconstruction and wavelet- or TV-based image deblurring.

II. BASIC INGREDIENTS

Variable splitting separates terms by introducing an auxiliary variable, while augmented Lagrangian methods enforce the resulting equality constraint without requiring an infinite penalty. ADMM alternates minimization and multiplier updates.

  • A. Variable Splitting: Variable splitting introduces v as the argument of f2 and imposes v = Gu, potentially making the constrained problem easier to solve.
  • A. Variable Splitting: Quadratic-penalty splitting alternates over u and v while increasing α, but the subproblems become increasingly ill-conditioned as α grows.
  • B. Augmented Lagrangian: Augmented Lagrangian methods alternate minimization over the primal variable with multiplier updates for linear equality constraints.
  • B. Augmented Lagrangian: Unlike quadratic-penalty methods, augmented Lagrangian methods do not require µ to approach infinity to guarantee convergence to the constrained solution.
  • B. Augmented Lagrangian: Completing the square yields an equivalent ALM/MM form with a quadratic residual term and a scaled dual-variable update.

C. ALM/MM for Variable Splitting and ADMM

ADMM solves the split constrained problem by alternating approximate minimizations over u and v, followed by a dual update. Under full-column-rank and convexity assumptions, the resulting sequence converges when a solution exists.

  • The split problem combines f1(u) and f2(v) while enforcing the linear relation represented by G.
  • The joint minimization can contain a non-separable quadratic term and non-smooth terms, motivating nonlinear block-Gauss-Seidel alternation.
  • ADMM alternates minimization over u and v, then updates d using the residual G u_k+1 − v_k+1.
  • The convergence theorem assumes full column rank for G and closed, proper, convex objective functions.
  • If problem (6) has a solution, the sequence {u_k} converges to a solution; otherwise, at least one of {u_k} or {d_k} diverges.
  • Exact minimization is unnecessary when the errors in the subproblems are absolutely summable, and exact ADMM satisfies the theorem with zero errors.

D. A Variant of ADMM

The paper generalizes ADMM to problems with multiple convex functions by splitting variables so the minimizations decouple, yielding an implementable ADMM-2 framework with convergence under a full-column-rank condition.

  • The generalized problem uses J closed, proper, convex functions connected through arbitrary matrices H(j).
  • The formulation maps the multi-function problem into the standard two-function ADMM structure by stacking auxiliary variables.
  • The f1 = 0 choice makes one ADMM step a quadratic minimization with a unique solution when G has full column rank.
  • The minimizations over u(1), ..., u(J) decouple and can therefore be solved independently.
  • Convergence is guaranteed when G has full column rank, while the minimizations may be solved inexactly if their errors are absolutely summable.
  • The other updates use Moreau proximal mappings, and the resulting ADMM-2 framework is presented as an algorithm.

III. PROPOSED METHOD

The paper instantiates ADMM-2 for constrained inverse problems using a regularizer and an ellipsoidal data-fidelity set, producing C-SALSA-1 with a convergence guarantee.

  • The algorithm is applied to the basic constrained problem and its analysis formulation.
  • The feasible set is E(ε, B, y) = {x ∈ R^n : ∥Bx − y∥2 ≤ ε}, an ellipsoid that may be unbounded when B is singular.
  • The constraint is represented through the indicator function ιE(ε,I,y), with E(ε, I, y) being a closed ε-radius Euclidean ball centered at y.
  • The regularizer step uses an efficiently computable Moreau proximal map, including soft thresholding for the ℓ1 norm and fast TV denoising methods.
  • The constraint step is projection onto the closed ε-radius ball, and convergence follows when the associated matrix has full column rank.
  • C-SALSA-1 results from instantiating ADMM-2 for the constrained problem.
  • For problem (2), C-SALSA-1 is obtained by replacing B with BW.

B. Problem (3)

The analysis formulation is handled by a closely related C-SALSA-2 algorithm, whose convergence depends on a full-column-rank condition supported by tight-frame assumptions.

  • Problem (3) is rewritten as an unconstrained problem and solved using an ADMM algorithm called C-SALSA-2.
  • C-SALSA-2 assumes that P is the analysis operator of a 1-tight Parseval frame, satisfying PHP = I.
  • The algorithm’s linear-system step is similar to the corresponding step in C-SALSA-1 under this frame assumption.
  • Convergence holds when φ and ιE(ε,I,y) are closed, proper, convex functions and the associated matrix has full column rank.

C. Solving (35)

The paper exploits matrix structure to solve C-SALSA linear systems efficiently, using FFT-based diagonalization for convolution and frame identities or SMW for analysis operators.

  • Explicitly forming B may be infeasible in imaging problems, but the required inverse can sometimes be computed with O(n log n) cost.
  • For deconvolution, B is diagonalized by the DFT as B = UᴴDU, reducing inversion to diagonal operations and FFT products.
  • The diagonal inversion costs O(n), while products by U and Uᴴ cost O(n log n) using the FFT.
  • For synthesis formulations with a Parseval frame W, the Sherman-Morrison-Woodbury lemma avoids directly diagonalizing BW.
  • The resulting computation costs O(n log n) or the cost of multiplying by W and Wᴴ.
  • Many tight frames used in image processing provide fast transforms, including translation-invariant wavelets and curvelets with O(n log n) cost.
  • For pixel-loss analysis formulations, B is a row-subset of the identity and satisfies BBᴴ = I, enabling another SMW-based simplification.

3) Missing Pixels: Image Inpainting:

C-SALSA exploits structured observation and frame operators to keep inpainting and related image-recovery iterations computationally efficient. Across broad regularizer and frame choices, its per-iteration complexity is O(n log n).

  • Missing-observation structure makes the relevant linear system diagonal, so a key C-SALSA update costs O(n).The diagonal entries are 1 or 1/2.
  • For frame-based synthesis, the same update uses the frame and observation operators through the SMW formula.The derivation uses BB^H = I and WW^H = I.
  • Inpainting updates reduce to synthesizing, masking, and recomputing representation coefficients, with cost typically O(n log n).The cost is governed by products with W and W^H.
  • Partial Fourier observations support the same complexity because the sampling mask is diagonal and FFT products cost O(n log n).This applies to both analysis and synthesis formulations.
  • Other algorithm lines use fast operator products, linear-time additions and projections, and proximal maps that are linear-time in common cases.The ℓ1 proximal map is a soft threshold with O(n) cost; TV proximal computation uses iterative algorithms that are usually O(n) per iteration.
  • For a wide choice of regularizers and frame representations, C-SALSA has O(n log n) computational complexity.

IV. EXPERIMENTS

The experiments compare C-SALSA with available state-of-the-art methods across deconvolution, inpainting, and partial-Fourier reconstruction. In deconvolution benchmarks, C-SALSA is reported as almost always clearly faster than the fastest competitor for both wavelet- and TV-based regularization.

  • IV. EXPERIMENTS: The study compares speed on deconvolution, inpainting, and partial-Fourier reconstruction problems, including MRI-related acquisition.
  • IV. EXPERIMENTS: Experiments average operator calls, iterations, CPU times, and MSE over 10 runs, using a protocol that matches C-SALSA against competitors’ stopping outcomes.
  • A. Image Deconvolution with wavelets: Five Cameraman deblurring problems compare C-SALSA with SPGL1 and NESTA for synthesis, and with NESTA for analysis.The experiments focus on computational speed because restored-image quality and SNR improvements were visually and numerically similar to prior results.
  • A. Image Deconvolution with wavelets: The wavelet experiments use redundant four-level Haar frames and also evaluate orthogonal Haar bases in synthesis and analysis settings.
  • B. Image Deblurring with Total Variation: The same five deblurring problems are evaluated with isotropic discrete TV regularization, whose proximal map uses five iterations of Chambolle’s algorithm.
  • B. Image Deblurring with Total Variation: C-SALSA is almost always clearly faster than the fastest competing algorithm for wavelet- and TV-regularized image deconvolution.

C. MRI Image Reconstruction

The MRI experiment reconstructs a noisy 128 × 128 Shepp–Logan phantom from 22 radial Fourier lines using TV regularization. C-SALSA is reported as faster than NESTA while achieving comparable reconstructed-image mean squared error.

  • MRI reconstruction: The experiment reconstructs a 128 × 128 Shepp–Logan phantom from 22 radial lines of its 2D discrete Fourier transform.The observations contain circular complex Gaussian noise with variance σ2 = 0.5 × 10−6.
  • MRI reconstruction: TV regularization is implemented through a Moreau proximal mapping computed with 10 iterations of Chambolle’s algorithm.
  • MRI reconstruction: C-SALSA produces an estimate that is visually indistinguishable from NESTA’s estimate.
  • MRI reconstruction: C-SALSA is faster than NESTA while achieving comparable values of mean squared error for the reconstructed image.The comparison includes operator calls, iterations, CPU time, and mean squared error.
  • Additional imaging experiments: The paper evaluates the same C-SALSA framework on deconvolution, high-dynamic-range TV reconstruction, and image inpainting benchmarks.The inpainting experiment removes 40% of the Cameraman image’s pixels and adds Gaussian noise at 40 dB SNR.
Loading 0912.3481v1…