Question about Parfor-loop

I do nonlinear fitting within a nested for loop as shown below:

   for i=1:200
       for j=1:170
            for k=1:50
                %do nonlinear fitting
                [p,res] = lsqnonlin()
                params(i,j,k,1:5)=[p,fix_par(i,j,k),res];
            end
        end
    end

The above loop is extremely slow and I am trying to find ways to speed it up. I've read about the parlour-loop (which I am not familiar at all - I'm new to programming). I tried to use it but I got the following error "The variable params in a parfor cannot be classified." I tried something as below:

parfor i=1:200
       parfor j=1:170
            parfor k=1:50
                %do nonlinear fitting
                [p,res] = lsqnonlin()
                params(i,j,k,1:5)=[p,fix_par(i,j,k),res];
            end
        end
    end

[Adam]

This page should help.

You cannot use nested parfor loops though, you can only use it on one loop (usually the outer loop is best if it is possible to parallelise at this level). This is intuitive from the fact that if the outer loop is using all workers then the loop nested inside it clearly cannot spread the work to all workers within each worker as they are all already being used.

[Steven Lord]

Why do you need to perform 1.7 million individual fits? Perhaps there's a way to achieve your ultimate goal without calling lsqnonlin almost 2 million times. If you describe what you're trying to do we may be able to offer suggestions about a faster alternative approach.

[Mar]

My data are 4D images with dimensions [200 170 50 31]. I am doing pixel-wise fitting that's why I am calling lsqnonlin so many times. Is there any alternative way to suggest?

Answer 1

[Abhishek]

Hi @Mar,

You can definitely speed up your nonlinear fitting using 'parfor', and it’s great that you’re exploring parallel computing.

In MATLAB, ‘parfor’ can only be used on a single loop level (typically the outermost one), but you can still parallelize all elements by flattening your 3D loop into a single index using ‘ind2sub’.

The below code will help to accomplish the workaround:

N = [200, 170, 50]; N1 = 7;
M = prod(N);
params = zeros(M, N1);  
parfor idx = 1:M
         [i, j, k] = ind2sub(N, idx); 
         [p, res] = lsqnonlin(@(x)  ModelFunctions(x, i, j, k), x0, lb, ub, options); 
         params(idx, :) = [p, fix_par(i, j, k), res]; 
end
params = reshape(params, [N,N1]);

This way, the workload is nicely distributed across workers, and you'll likely see a significant speed-up depending on your system resources. Just ensure that ‘fix_par’ is accessible in the correct format, ideally sliced by indices inside the loop.

Also make sure ‘ModelFunction’ is written to accept ‘(i, j, k)’ as context if needed. Moreover, avoid using shared/global variables; pass all required data explicitly.

This workaround efficiently distributes the fitting tasks across workers, often resulting in substantial runtime reductions for large data sets.

You can find additional information about ‘ind2’sub’ in the MATLAB’s official documentation: https://www.mathworks.com/help/releases/R2024b/matlab/ref/ind2sub.html

[Walter Roberson]

params(i, j, k, :) = [p, fix_par(i, j, k), res]; 

That code will not work. Output variables must have one position indexed by a simple expression involving the parfor index, and all other indices to the output variable must be constants or the : operator. Furthermore, the output variable can be assigned to only once.

You can do something like

N = [200, 170, 50]; N1 = 7;
M = prod(N);
params = zeros(M, N1);  
parfor idx = 1:M
         [i, j, k] = ind2sub(N, idx); 
         [p, res] = lsqnonlin(@(x)  ModelFunctions(x, i, j, k), x0, lb, ub, options); 
         params(idx, :) = [p, fix_par(i, j, k), res]; 
end
params = reshape(params, [N,N1]);

[Abhishek]

Updated, thanks for the suggestion.