title: Implement a fast optimization algorithm using fixed point method in matlab
body: <p>I am implementing a fast optimization algorithm using fixed point method in matlab. The goal of that method is that find optimal value of <code>u</code>. Denote <code>u={u_i,i=1..2}</code>. The optimal value of u can be obtained as following steps: </p>
  
  <p><img src="https://i.stack.imgur.com/PtUzH.png" alt="enter image description here"></p>
  
  <p>Sorry about my image because I cannot type mathematics equation in here.
  To do that task, I tried to find u follows above steps. However, I don't know how to implement the term \sum_{j!=i} (u_j-1) in equation 25. This is my code. Please see it and could you give me some comment or suggestion about my implementation to correct them. Currently, I tried to run that code but it give an incorrect answer. </p>
  
  <pre><code>function u = compute_u_TV(Im0, N_class)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Initialization
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  theta=0.001;
  gamma=0.01;
  tau=0.1;
  sigma=0.1;
  N_class=2; % only have u1 and u2
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Iterative segmentation process
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  for i=1:N_class
      v(:,:,i) = Im0/max(Im0(:)); % u between 0 and 1.   
      qxv(:,:,i) = zeros(size(Im0));
      qyv(:,:,i) = zeros(size(Im0));
      u(:,:,i) = v(:,:,i);
  for iteration=1:10000
      u_temp=u;
      % Update v
      Divqi = ( BackwardX(qxv(:,:,i)) + BackwardY(qyv(:,:,i)) );
      Term = Divqi - u(:,:,i)/ (theta*gamma);
      TermX = ForwardX(Term);
      TermY = ForwardY(Term); 
      Norm = sqrt(TermX.^2 + TermY.^2);
      Denom = 1 + tau*Norm;
      %Equation 24
      qxv(:,:,i) = (qxv(:,:,i) + tau*TermX)./Denom;
      qyv(:,:,i) = (qyv(:,:,i) + tau*TermY)./Denom;  
      v(:,:,i) = u(:,:,i) - theta*gamma* Divqi;  %Equation 23   
      % Update u  
      u(:,:,i) = (v(:,:,i) - theta* gamma* Divqi -theta*gamma*sigma*(sum(u(:))-u(:,:,i)-1))./(1+theta* gamma*sigma);
      u(:,:,i) = max(u(:,:,i),0); 
      u(:,:,i) = min(u(:,:,i),1);
      check=u_temp(:,:,i)-u(:,:,i);
      if(abs(sum(check(:)))&lt;=0.1)
          break;  
      end
  end
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Sub-functions- X.Berson
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function [dx]=BackwardX(u);
  [Ny,Nx] = size(u);
  dx = u;
  dx(2:Ny-1,2:Nx-1)=( u(2:Ny-1,2:Nx-1) - u(2:Ny-1,1:Nx-2) );
  dx(:,Nx) = -u(:,Nx-1);
  
  function [dy]=BackwardY(u);
  [Ny,Nx] = size(u);
  dy = u;
  dy(2:Ny-1,2:Nx-1)=( u(2:Ny-1,2:Nx-1) - u(1:Ny-2,2:Nx-1) );
  dy(Ny,:) = -u(Ny-1,:);
  
  function [dx]=ForwardX(u);
  [Ny,Nx] = size(u);
  dx = zeros(Ny,Nx);
  dx(1:Ny-1,1:Nx-1)=( u(1:Ny-1,2:Nx) - u(1:Ny-1,1:Nx-1) );
  
  function [dy]=ForwardY(u);
  [Ny,Nx] = size(u);
  dy = zeros(Ny,Nx);
  dy(1:Ny-1,1:Nx-1)=( u(2:Ny,1:Nx-1) - u(1:Ny-1,1:Nx-1) );
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % End of sub-function
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  </code></pre>