How to remove for loops in this code?

Hi, I am new to matlab. i am running the following code .It is taking a very long time to process. Can anyone suggest how to rewrite it to reduce the computation time. arr1 = 4d array with size (512 512 122 4)

parfor n=1:slice
  n    
 for i=1:width
    for j=1:height
      options = optimset('Algorithm', 'levenberg-marquardt', 'display','off');
      warning('off');            
        St = squeeze(arr1(i,j,n,:));
        if St(1)>50
            St = St./St(1);
            St = St(:)';
            bat = 1;
            c_bat = St(1);
            %pl_parm=[TTP,slop1,slop2];
            pl_parm=[2;0;.1];
            xdata = [1:length(St)];
            lb = [2;0;-1];
            ub = [4;10;1];

              [par_PL,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(@PLmodel,pl_parm,xdata,St,lb,ub,options);
              TTP(i,j,n)=par_PL(1);
              Slope1(i,j,n)=par_PL(2);
              Slope2(i,j,n)=par_PL(3);

          end
      end
   end
end

Answer 1

[Marc Jakobi]

If lsqcurvefit() can't be run on matrices, you may not be able to reduce the loop nesting. But if it is a lot of data, you can decrease time by a bit by taking variables that don't change out of the loop. Setting them every time is not efficient. Also, try to avoid displaying things (like n) in the command window. Command window outputs significantly increase the time it takes to perform computations.

options = optimset('Algorithm', 'levenberg-marquardt', 'display','off');
warning('off');  
bat = 1;
%pl_parm=[TTP,slop1,slop2];
pl_parm=[2;0;.1];
lb = [2;0;-1];
ub = [4;10;1];
  parfor n=1:slice 
   for i=1:width
      for j=1:height    
          St = squeeze(arr1(i,j,n,:));
          if St(1)>50
              St = St./St(1);
              St = St(:)';
              c_bat = St(1);
              xdata = [1:length(St)];
                [par_PL,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(@PLmodel,pl_parm,xdata,St,lb,ub,options);
                TTP(i,j,n)=par_PL(1);
                Slope1(i,j,n)=par_PL(2);
                Slope2(i,j,n)=par_PL(3);
            end
        end
     end
  end

Answer 2

[Walter Roberson]

lsqcurvefit is not itself vectorizable.

You have bounds constraints, so lsqcurvefit is going to turn this into a lsqnonlin call with trust-region-reflective selected.

You could pull out all the hyperpages where arr1(:,:,:,1) > 50, and arrange those as a series of rows or columns in a 2D array. The question then becomes whether there is a way that you could ask to perform the fitting over the entire 2D array at once. I would have to think more about whether there would be a practical way to transform a least-squared model to have one more dimension; it is not obvious to me that it is possible in general. However, if we had the code for your PLmodel and we could see what it is you are fitting, we just might be able to think of a different approach.

[Prativa sahoo]

please find code for PLmodel below. it is a piecewise linear model fits into two line

%par_PL=[besline,BAT,TTP,slop1,slop2]
  %y=f(x) = c                     if x<BAT
  %         c+m1(x-BAT)           if BAT<=x<TTP
  %         c+m1(TTP-BAT)+m2(x-TTP) if x>TTP

    function [fx] = PLmodel(parm,X)

    C=1;
    BAT=1;
    TTP = parm(1);
    slop1 = parm(2);
    slop2 = parm(3);
        for i=1:length(X)
            if X(i)<BAT 
                fx(i)=C;
            elseif X(i)<TTP
               fx(i)=C+slop1*(X(i)-BAT);
            else fx(i)=C+slop1*(TTP-BAT)+slop2*(X(i)-TTP);
            end
        end
    end

[Walter Roberson]

Vectorize that function to improve performance.

function [fx] = PLmodel(parm,X)

    C=1;
    BAT=1;
    TTP = parm(1);
    slop1 = parm(2);
    slop2 = parm(3);
    fx = zeros(size(X));
    mask = X < BAT;
    fx(mask) = C;
    mask = BAT <= X & X < TTP;
    fx(mask) = C + slop1 * (X(mask) - BAT);
    mask = X >= TTP;
    fx(mask) = C + slop1 * (TTP-BAT) + slop2 * (X(mask) - TTP);

[Prativa sahoo]

Thank you for your help.