Preconditionnate conjugate gradient

Preconditionner Conjugate gradient optimisation

PCONJGRAD - Preconditionner Conjugate gradient optimisation

[m] = conjGrad('calcAx', init, b, options, ...)

provide a minimizer 'm' of the criterion J(x) = 1/2 x^tAx + b^tx as
the solution of the linear system Ax = b, computed by a conjugate
gradient descente algorithm. This implementation is adapted, but
absolutely not restricted, to inverse problems where the criterion
take the form

J(x) = ||y - Hx||^2 + l||Dx||^2, and a algorithm is available to
compute Hx or H^te (x is the unkown, H the direct model, y is data
and D the regularisation.

In these cases, A = 2(H^tH + lD^tD) and b = H^ty.

Anyway, the code is quite general since it accept any function that
compute the product Ax, whatever A.

PARAMETERS

'calcAx' - the name of the function to compute the matrix vector
product Ax. Consequently, it's not necessary to compute and store
matrix A, only the product is necessary (think about the
convolution). The name can be anything, but must correspond to a
callable function. The first argument of 'calcAx' MUST BE the
vector x.

init - the starting point of the optimisation x^(0)

b - the term b = H^t y. Must be provided.

options - the options of the algorithm. See section below.

... - all the remaning argument are passed to the function 'calcAx'
which is call like this calcAx(x,...)

OPTIONS

The variable options is matlab structure. The field are listed
below. Options with '*' are necessary. Others are optional.

thresold (*) - the stoping criterion

maxIter (*) - the maximum number of iteration (the algorithm is
automatically stopped when the iteration equal the dimension of x)

numfig - this option must be a integer. In these case, the current
minimizer and residual are displayed as image of real number in the
figure 'numfig'.

An example is

cgoptions.thresold = 1e-6; cgoptions.maxIter = 50;

Copyright (C) 2013 by François Orieux <[e-mail]>