A comparison of some deconvolution priors

Continuing with experiments based on ADMM I implemented and experimented with a few different deconvolution priors.  I chose several $\ell_2$ and $\ell_1$ priors based on sparse gradients, sparse curvature and simply the norm of the solution vector as well as an $\ell_{0.5}$ HyperLaplacian prior.

For all results the ground-truth signal is shown in black, the blurred and corrupted signal in red and the deblurred signal in blue.  As a ground-truth signal I'm using a simple piecewise linear signal with a slanted region, some constant regions and discontinuous jumps.  These features help to show the issues with different regularizers.  The input signal is 128 samples long, blurred by a Gaussian filter with $\sigma=4$ and then corrupted by Gaussian noise with $\sigma=2.5%$ and $\mu=0$. 

For the $\ell_2$ results I'm solving for the optimal solution directly in the Fourier domain, while for the $\ell_1$ and $\ell_0.5$ I'm using ADMM, solving the data-subproblem in the Fourier domain and the splitting variable subproblem with shrinkage/proximal operators.  For all problems I've played with the prior weights a bit to try to get good results, but don't claim they are optimal.

$\ell_2$ norm of solution vector

First up is the $\ell_2$ norm of the solution.  This does sharpen up the signal somewhat, but has bad ringing artifacts.  These are expected for the $\ell_2$ solution, particularly given the high noise levels in the corrupted signal.  The relatively large blur also leaves lots of freedom for low-frequency ringing.

$\ell_2$ norm of solution first derivative

Results for the $\ell_2$-norm of the signal derivatives do a bit better, particularly for the slanted segment, but still show low-frequency ringing.  The same is true for the $\ell_2$-norm of the second derivative, although the ringing is a bit more pronounced in the flat regions:

$\ell_2$ norm of solution second derivative

Moving on to the $\ell_1$ priors, the total-variation prior is simply the $1$-norm of the signal first derivatives. It favours piecewise constant solutions which do well on the large flat regions but introduce staircase artifacts for the slanted region:

$\ell_1$ norm of solution first derivative

The total variation prior produces a much sharper signal than the $\ell_2$ priors, however for the slanted region it also introduces spurious features in the slanted region that are hard to distinguish from (correct) features elsewhere in the signal.  This occurs because the assumption of the total variation prior is that the input is piecewise constant, which is not true for this signal.  A better assumption is that the solution is piecewise linear, leading to the sparse-curvature prior.  The sparse-curvature assumes that the second derivative of the signal is non-zero in only a few locations:

$\ell_1$ norm of solution second derivative

The sparse-curvature prior does much better in the slanted region than the total variation prior and still improves the signal, but the edges are not as sharp at discontinuities.  Often people combine these priors to try to achieve a compromise between the two behaviours.

The final prior is the Hyper-Laplacian prior, which is simply the $\ell_{0.5}$-norm of the signal first derivatives:

$\ell_{0.5}$ norm of solution first derivative

The Hyper-Laplacian prior is non-convex, making it difficult to optimize for.  Here it appears that the optimization got stuck: some edges are extremely sharp, however the trailing edge of the high-amplitude constant region is lost entirely, possible because of the slanted region adjacent to it.  I played with the parameters quite a bit but did not get a result better than this. 

Matlab example code generating these results is available at my GitHub repository:

https://github.com/jamesgregson/Matlab-deblurring-example

Matlab/MEX utility for loading VTK Rectilinear Grid files (.vtr files)

The title says it all, I've written some basic code for loading VTK rectilinear grid files into Matlab.  The code supports uniformly spaced meshes in up to four dimensions for both point and cell data.

You can download the code from my Github repository:
https://github.com/jamesgregson/matlab_vtr

C++/Mex Image Deblurring using ADMM

I've posted some sample code on Github for performing image deblurring in Matlab using Mex.  I wrote it as a way to play around with the ADMM algorithm for sparse signal reconstruction, as described in Stephen Boyd's ADMM paper, as well as to get some experience using C++ code from Matlab.

The code uses matrix-free linear operators to evaluate the forward and reverse measurement (blurring) process, solving the data-term subproblem with gradient descent.  While you can do this much faster using FFTs for deconvolution, the code should generalize quite well to tomography with non-orthographic projections.  I've tried to comment it well so that I can figure out what I did if I ever need to pick it up again.

In the image above, the left plot is the ground-truth image and the right is blurred by a point-spread function and corrupted by Gaussian noise.  After running the demo, the following results are obtained:

The reconstructed image is shown on the right, while a slice through the center scanline shows the reconstruction (blue dots), ground-truth (solid black line) and input (red line).  From the results it's clear that the method does well at preserving discontinuities even for fairly large blurs and noise levels.

The demo is pretty quick to run and is helpful for choosing parameters since you can see how the behavior of the method changes with the regularizer weight and the penalty term weight.

You can download the code from https://github.com/jamesgregson/MexADMMDeblurDemo

Matlab Signal Deblurring & Denoising Example

To date my research has been largely focused on inverse problem such as tomography or image deblurring.  These problems are often highly under-determined and so must include strong priors to obtain good solutions and finding efficient solvers for these priors is challenging.  A labmate recently pointed me towards the ADMM method, which splits the full problem into coupled sub-problems.  It has a number of advantages:

• Flexibility - ADMM handles a number of non-trivial problems in a common framework
• Efficiency - Often the subproblems have efficient, often embarassingly-parallel, solvers

An excellent introduction to the method can be found in the paper: Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers by Boyd et al.. It's a pretty gentle and practical introduction, with numerous example problems.

I've quickly implemented ADMM for combined deblurring and denoising of 1D input signals using the total-variation regularization in a Generalized-Lasso problem definition.  Unlike the Boyd paper, I've chosen to use Landweber iterations to solve the data subproblem as these are commonly used in large-scale deblurring and tomography.  My well-commented sample implementation (see the end of this post) allows all the method parameters to be tweaked to see the effect of using inexact solves, different penalty parameters and different regularization weights.

As an example, the image below shows a 100 sample square-wave signal, blurred by a Gaussian with standard deviation of 3 and corrupted by Gaussian noise with a sigma of 5%.  Only two inner Landweber iterations were used for each outer iteration.  The red line shows the original blurred and noisy input, while the blue line shows the reconstructed result after 100 iterations.

Overall I find that the method works as advertised. It's fairly easy to implement, either in Matlab or in C/C++, can handle 1-norms of general linear regularizers and converges quickly.  It also allows matrix-free solvers to be used for the data-subproblem, while the solve for the prior is embarrassingly-parallel, so it should scale quite well too.

You can get the implementation that produced this plot below. Note that you may get slightly different results since the noise is generated randomly.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                                                                     %%%
%%% Sample signal denoising & deblurring using ADMM                 %%%
%%% code written by James Gregson (james.gregson@gmail.com), 2013       %%%
%%% Use the code for whatever you'd like                                %%%
%%%                                                                     %%%
%%% Demonstrates performing Total-Variation (TV) denoising and          %%%
%%% deblurring of a 1D input signal using the ADMM iterative scheme     %%%
%%% applied to a Generalized-Lasso problem.  More information on the    %%%
%%% approach can be found in [1].  The approach from [1] is modified    %%%
%%% slightly by using gradient-descent (Landweber iterations) to solve  %%%
%%% the first subproblem in place of a direct solver or iterative       %%%
%%% method such as conjugate gradient.  Matrix-free Landweber           %%%
%%% iterations are commonly used for large-scale linear inverse         %%%
%%% problems and this sample code allows the accuracy of the            %%%
%%% subproblem solves to be adjusted to see the effect on the final     %%%
%%% reconstructions.                                                    %%%
%%%                                                                     %%%
%%% [1] Boyd et al., Distributed Optimization and Statistical Learning  %%%
%%%     via the Alternating Direction Method of Multipliers, 2010       %%%
%%%                                                                     %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all;
clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Problem Parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N           = 100;     % Signal sample count
lambda      = 0.1;     % Total-variation weight
sigma       = 3.0;     % PSF sigma for generating blurred input
noise       = 0.05;    % Gaussian-noise sigma

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ADMM Parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rho         = 1.0;     % ADMM constraint weight 0 < rho <= 2
outer_iters = 100;     % Number of ADMM iterations
inner_iters = 2;       % Number of Landweber steps per ADMM iteration
relax       = 1.9;     % Under-relaxation factor for Landweber steps [0,2]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Construct System PSF & Image Formation Model %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

off = -ceil(sigma*3):ceil(sigma*3);  % psf pixel offsets
psf = exp( -off.^2/sigma^2 );        % psf values for offsets
psf = psf/sum(psf);                  % normalize to 1.0

% generate psf matrix by setting diagonals based on 1D psf above
M = zeros( N, N );                   
for i=1:numel(psf),
    M = M + diag( psf(i)*ones(N-abs(off(i)),1), off(i) );
end
M = sparse( M );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Generate noisy and blurred synthetic input %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

input = zeros( N, 1 );                  % start with zeros
input(floor(N/4):ceil(3*N/4),:) = 1.0;  % make a center region at 1.0
blur = M * input + noise*(randn(N,1));  % blur the input by the psf


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Construct the difference matrix that computes image gradients %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

A = sparse( -diag( ones(N,1), 0 ) + diag( ones(N-1,1), 1 ) );
A(N,N-1)=1; % use a backward difference for the final point

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Setup ADMM and perform iterations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute the eigenvalues of the first subproblem
% system matrix and use the largest to define a
% Landweber step size, for large system a a power-
% iteration should be used instead.
tmp = eigs(M'*M + rho*A'*A);
step = relax/tmp(1);

% initialize the ADMM variables
x = blur;             % intrinsic (sharp) signal
z = zeros( N, 1 );    % splitting variable
u = zeros( N, 1 );    % scaled Lagrange multipliers

% define an anonymous shrinkage operator to implement
% the second sub-problem solve
shrink = @(kappa,x) max( abs(x)-kappa, 0 ).*sign(x);

% perform outer_iters outer iterations, plot the
% solver progress as the iterations proceed
for k=1:outer_iters,
    fprintf( 1, 'iteration %d of %d\n', k, outer_iters );
    plot( x, 'b-+' );
    drawnow;
 
    % define an anonymous function returning the gradient
    % of the first subproblem w.r.t. x, holding z and u
    % fixed, then perform gradient-descent (Landweber 
    % iterations).
    gradF = @(x) M'*(M*x) - M'*blur + rho*A'*( A*x - z + u );
    x = gradient_descent( gradF, x, step, inner_iters );
    
    % solve the second sub-problem using the anonymous
    % shrinkage operator
    z = shrink( lambda/rho, A*x + u );
    
    % update the scaled Lagrange multipliers
    u = u + A*x - z;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Setup ADMM and perform iterations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

hold on;
plot( blur, 'r-o' );
plot( x, 'b-o' );

Symmetric Discretization of Linear Boundary Conditions

When solving large-scale PDEs, it's generally preferable (where possible) to have a symmetric discretization so that the conjugate gradient method can be used.  Conjugate gradient uses relatively little memory compared to more general purpose methods (like GMRES) and can also take advantage of symmetric preconditioners (like incomplete LDL^T) that are faster and use lower memory than preconditioners such as incomplete LU.

However applying boundary conditions to the problem can be an issue, since it is easy to end up with a non-symmetric matrix unless the constrained variables are eliminated from the system matrix.  This means that solvers like GMRES have to be used, which is costly in terms of computation and memory.

The following Matlab script gives an example for how to apply linear boundary conditions via an auxiliary boundary condition matrix suggested by a labmate of mine.  It shows both pin (Dirichlet) constraints and gradient (Neumann) boundary conditions at the domain boundary and in the interior.  The only downside is that Neumann constraints are a bit peculiar to specify since one of the variables must be eliminated from the boundary condition matrix.  However this can usually be specified fairly easily, and certainly more easily than eliminating the variable from the system matrix.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Example for showing symmetric discretization  %%
%% of general linear boundary conditions.        %%
%% (c) James Gregson 2013                        %%
%% Please send questions or comments to:         %%
%% james.gregson@gmai.com                        %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% solve a 1D Laplace equaiton with 8 mesh points
% subject to the constraint that the leftmost 
% grid point has value 1.0 and the rightmost 
% grid point has value 2.0 and the fifth is 
% constrained to the value of the fourth minus 0.2
%
% then solves the same problem, but with a biharmonic
% operator, to show that the result is different 
% from just pinning the values; i.e. the biharmonic
% energy is applied to the constrained grid points. 
% If this weren't the case, both problems would have 
% the same solution.

% the approach is to use a standard symmetric 
% discretization of the Laplace operator A for all
% grid points (including those constrained by
% boundary conditions) and a second boundary
% condition matrix that enforces the constraints
% specifically, the boundary conditions are
% represented as:
%
% u = B v + c
%
% where u has the same dimensions as the Laplace
% operator and B has one column for every
% unconstrained grid-point.  c is a vector of
% constants. this allows the problem: 
%
% A u = f
%
% to be written as:
%
% A*B*v = f - A*c
%
% however this is system is rectangular and 
% non-symmetric.  This can be remedied by 
% multiplying through by B', giving:
%
% B'*A*B*v = B'*(f-A*c)
% 
% which is square and symmetric.

% define the objective function, in this case a 
% 6-variable 1D Laplace equation. The Laplacian
% is discretized with a standard 3-point stencil
A = [  1, -1,  0,  0,  0,  0,  0,  0;
      -1,  2, -1,  0,  0,  0,  0,  0;
       0, -1,  2, -1,  0,  0,  0,  0;
       0,  0, -1,  2, -1,  0,  0,  0;
       0,  0,  0, -1,  2, -1,  0,  0;
       0,  0,  0,  0, -1,  2, -1,  0;
       0,  0,  0,  0,  0, -1,  2, -1;
       0,  0,  0,  0,  0,  0, -1,  1 ];

% we are solving the Laplace equation, so the
% right hand side is zero
f = [ 0, 0, 0, 0, 0, 0, 0, 0 ]';

% the boundary condition matrix.  
B = [  0,  0,  0,  0,  0; % first grid point, set to c[1]
       1,  0,  0,  0,  0;
       0,  1,  0,  0,  0;
       0,  0,  1,  0,  0;
       0,  0,  1,  0,  0; % fifth grid point, set to c[4]-0.2
       0,  0,  0,  1,  0;
       0,  0,  0,  0,  1;
       0,  0,  0,  0,  0 ]; % last grid point, set to c[8]

% the boundary condition constant matrix
c = [ 1, 0, 0, 0, -0.2, 0, 0, 2 ]';

% form the system with boundary conditions
M = B'*A*B
rhs = B'*( f - A*c );

% solve for the primary variables
x = M\rhs;

% use the primary variables to solve for the
% grid values
u1 = B * x + c;

% now solve a biharmonic problem with the same constraints
% in the same manner, note the A' factor in the  system
% matrix and right-hand-side.
M = B'*A'*A*B;
rhs = B'*( f - A'*A*c );
x = M\rhs;
u2 = B * x + c;

% plot the results
figure
hold on
plot( u1, '-ro' );
plot( u2, '-b+' );
legend( 'Laplace', 'Biharmonic' );  

The plot produced by this script is shown below:

Both solves reproduce the boundary conditions correctly, but the biharmonic solve finds a solution that is smoother than the Laplace solution around the interior constraints.