Large matrix minimization using fminunc

Dear all,

I am relatively new to optimization problems in general.

I am trying to implement differential image integration solution proposed in this paper: Thüring, Thomas, et al. "Non-linear regularized phase retrieval for unidirectional X-ray differential phase contrast radiography." Optics express 19.25 (2011): 25545-25558.

To give you a background, I measured a differential image, where each row is the derivative on the integrated profile f along the rows. I want to find f. Line integration using cumsum doesn't provide a good solution, because noise introduces stripes artifacts (see the image).

The idea in the paper is that, given a differential image, they integrate it along the rows minimizing the difference between the differential experimental profile and the integrated profile. However, a regularization term taking into account the derivative along the columns is also added to impose some degree of continuity in the vertical direction, preventig stripes.

To better explain, they look for want to find :

where and are the derivatives along rows and columns, respectively; f is the unknown integratedimage I aim to find and ϕ is the experimental differential image. Note that in the second term you can also use the norm.

My first implementation of this idea is based on fminunc. Find it below.

p=phantom('Modified Shepp-Logan',512);
pdiff=imresize(diff(padarray(p,[0,1],'replicate','post'),1,2),[128,128]);
noise = 0.02 * randn(size(pdiff)); %sig is your standard dev.
pdiff_noise = pdiff + noise ;

% normal integration
integratedProfile=cumsum(pdiff_noise,2);

options=optimoptions(@fminunc,'Display','iter-detailed','UseParallel',true,'MaxIterations',35);

% regularized integration
startingPoint=zeros(size(pdiff_noise));
fun=@(x)costFunction(pdiff_noise,x);
integratedProfile_optmized=fminunc(fun,startingPoint,options);


function difference=costFunction(ep,op)
L=0.2;
dep_o=diff(padarray(op,[0,1],'replicate','post'),1,2);
dep_v=diff(padarray(op,[1,0],'replicate','post'),1,1);
difference=sum(sum((ep-dep_o).^2))+L*sum(sum(dep_v.^2));
end

It works, but it's slow even running in parallel using 24 workers. Eventually, it gets a nice solution on a 128x128 test image. I think the problem behind speed is the large number of involved variables. Basically, each pixel is a variable to optimize and because the objective function returns a scalar fminunc must call it 16384 times.

The main problem is that as soon as I try to scale up to 256x256 I run out of memory. I guess because now I am requiring a 65536x65536 hessian matrix that doesn't fit in the RAM. How can I solve this problem? In the paper they report they used non-linear Conjugate Gradient algorithm. What is the solver using such a method? I need to work with images that are approximately 500x300 pixels.

Do you have any other solution to solve this problem in a reasonable time? I have approximately 3000 images to process.