title: How to minimize a cost function with Matlab when input variable is a large image: increase speed and prevent memory crash
body: <p>I am trying to implement a differential phase integration method described in this paper:
  <em>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</em>. </p>
  
  <p>Basically, it's a way to integrate a differential image across the columns only, while imposing some constraints on continuity across the rows to prevent stripe noise. 
  From a mathematical point of view, I want to minimize the following equation:</p>
  
  <p><a href="https://i.stack.imgur.com/fzc5n.gif" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/fzc5n.gif" alt="$${||D_xA-\phi||_2}^2+\lambda{||D_yA||_2}^2$$"></a> </p>
  
  <p>where ||.|| is the L2 norm, Dx is the derivative along the columns, Dy is the derivative across the rows, A is the unknown integrated matrix, lambda is a user-defined parameter and phi is the differential profile I measured. Note that for the Dy operator the L1  norm can also be used.</p>
  
  <p>I wrote down a code using <em>fminunc</em> as Matlab solver</p>
  
  <pre><code>pdiff=imresize(diff(padarray(p,[0,1],'replicate','post'),1,2),[128,128]);
  noise = 0.02 * randn(size(pdiff)); 
  
  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
  </code></pre>
  
  <p>It works using a 128x128 differential image. 
  The problem arises as soon as I try to work with a larger image. In particular, when I use a 256x256 matrix takes forever to make each iteration even using the parallel option and takes almost the entire RAM. 
  When I move to a matrix that is 512x512 I get this error </p>
  
  <pre><code>Requested 262144x262144 (512.0GB) array exceeds maximum
  array size preference. 
  
  Error in fminusub (line 165)
      H = eye(sizes.nVar);
  
  Error in fminunc (line 446)
     [x,FVAL,GRAD,HESSIAN,EXITFLAG,OUTPUT] =
     fminusub(funfcn,x, ...
  
  Error in Untitled (line 13)
  integratedProfile_optmized=fminunc(fun,startingPoint,options); 
  </code></pre>
  
  <p>Unfortunately, my final goal is to process approximately 3000 images of 500x500 size.</p>
  
  <p>I think I have understood that the crash problem is related to the size of the matrix and to the fact that each pixel is a variable. Therefore, Matlab needs to calculate a huge hessian that doesn't fit into the memory. 
  However, I don't really know how to solve it while also speeding up the processing. 
  Do you have any suggestions on how to work with large images? Is there another solver that may work in a faster way? Any mathematical approach to making the problem easier?</p>
  
  <p>Thanks!</p>