Conjugate gradient acceleration of iteratively re-weighted least squares methods

Iteratively Re-weighted Least Squares (IRLS) is a method for solving minimization problems involving non-quadratic cost functions, perhaps non-convex and non-smooth, which however can be described as the infimum over a family of quadratic functions. This transformation suggests an algorithmic scheme that solves a sequence of quadratic problems to be tackled efficiently by tools of numerical linear algebra. Its general scope and its usually simple implementation, transforming the initial non-convex and non-smooth minimization problem into a more familiar and easily solvable quadratic optimization problem, make it a versatile algorithm. However, despite its simplicity, versatility, and elegant analysis, the complexity of IRLS strongly depends on the way the solution of the successive quadratic optimizations is addressed. For the important special case of $\textit{compressed sensing}$ and sparse recovery problems in signal processing, we investigate theoretically and numerically how accurately one needs to solve the quadratic problems by means of the $\textit{conjugate gradient}$ (CG) method in each iteration in order to guarantee convergence. The use of the CG method may significantly speed-up the numerical solution of the quadratic subproblems, in particular, when fast matrix-vector multiplication (exploiting for instance the FFT) is available for the matrix involved. In addition, we study convergence rates. Our modified IRLS method outperforms state of the art first order methods such as Iterative Hard Thresholding (IHT) or Fast Iterative Soft-Thresholding Algorithm (FISTA) in many situations, especially in large dimensions. Moreover, IRLS is often able to recover sparse vectors from fewer measurements than required for IHT and FISTA.Comment: 40 page

Joint Total Variation ESTATICS for Robust Multi-Parameter Mapping

Quantitative magnetic resonance imaging (qMRI) derives tissue-specific parameters -- such as the apparent transverse relaxation rate R2*, the longitudinal relaxation rate R1 and the magnetisation transfer saturation -- that can be compared across sites and scanners and carry important information about the underlying microstructure. The multi-parameter mapping (MPM) protocol takes advantage of multi-echo acquisitions with variable flip angles to extract these parameters in a clinically acceptable scan time. In this context, ESTATICS performs a joint loglinear fit of multiple echo series to extract R2* and multiple extrapolated intercepts, thereby improving robustness to motion and decreasing the variance of the estimators. In this paper, we extend this model in two ways: (1) by introducing a joint total variation (JTV) prior on the intercepts and decay, and (2) by deriving a nonlinear maximum \emph{a posteriori} estimate. We evaluated the proposed algorithm by predicting left-out echoes in a rich single-subject dataset. In this validation, we outperformed other state-of-the-art methods and additionally showed that the proposed approach greatly reduces the variance of the estimated maps, without introducing bias.Comment: 11 pages, 2 figures, 1 table, conference paper, accepted at MICCAI 202

Symmetrization Techniques in Image Deblurring

This paper presents a couple of preconditioning techniques that can be used to enhance the performance of iterative regularization methods applied to image deblurring problems with a variety of point spread functions (PSFs) and boundary conditions. More precisely, we first consider the anti-identity preconditioner, which symmetrizes the coefficient matrix associated to problems with zero boundary conditions, allowing the use of MINRES as a regularization method. When considering more sophisticated boundary conditions and strongly nonsymmetric PSFs, the anti-identity preconditioner improves the performance of GMRES. We then consider both stationary and iteration-dependent regularizing circulant preconditioners that, applied in connection with the anti-identity matrix and both standard and flexible Krylov subspaces, speed up the iterations. A theoretical result about the clustering of the eigenvalues of the preconditioned matrices is proved in a special case. The results of many numerical experiments are reported to show the effectiveness of the new preconditioning techniques, including when considering the deblurring of sparse images

Semi-local scaling exponent estimation with box-penalty constraints and total-variation regularisation

We here establish and exploit the result that 2-D isotropic self-similar fields beget quasi-decorrelated wavelet coefficients and that the resulting localised log sample second moment statistic is asymptotically normal. This leads to the development of a semi-local scaling exponent estimation framework with optimally modified weights. Furthermore, recent interest in penalty methods for least squares problems and generalised Lasso for scaling exponent estimation inspires the simultaneous incorporation of both bounding box constraints and total variation smoothing into an iteratively reweighted least-squares estimator framework. Numerical results on fractional Brownian fields with global and piecewise constant, semi-local Hurst parameters illustrate the benefits of the new estimators

Cygnus A super-resolved via convex optimisation from VLA data

We leverage the Sparsity Averaging Reweighted Analysis (SARA) approach for interferometric imaging, that is based on convex optimisation, for the super-resolution of Cyg A from observations at the frequencies 8.422GHz and 6.678GHz with the Karl G. Jansky Very Large Array (VLA). The associated average sparsity and positivity priors enable image reconstruction beyond instrumental resolution. An adaptive Preconditioned Primal-Dual algorithmic structure is developed for imaging in the presence of unknown noise levels and calibration errors. We demonstrate the superior performance of the algorithm with respect to the conventional CLEAN-based methods, reflected in super-resolved images with high fidelity. The high resolution features of the recovered images are validated by referring to maps of Cyg A at higher frequencies, more precisely 17.324GHz and 14.252GHz. We also confirm the recent discovery of a radio transient in Cyg A, revealed in the recovered images of the investigated data sets. Our matlab code is available online on GitHub.Comment: 14 pages, 7 figures (3/7 animated figures), accepted for publication in MNRA