Why does performance of functions saturate with number of cores using parfeval but not with parfor?

I am developing an application that MUST take advantage of parallelization, and ideally offer real-time updates after each iteration, which makes use of parfeval prefarable. I believe the algorithm that I have developed is highly parallelizable (see attached for performance of 'WT_Ex_2_b' as a function of number of cores used in parfeval function). From 1 to 8 cores, the speedup factor agrees with theoretical expectation (Amdahl's Law with p=0.95), however, performance of my application saturates at 8 cores. This led me to create a dummy function (see attached script) to compare the performance of using parfor and parfeval as a function of number of cores. I discovered that the parfor version behaves quite similarly to theoretical expectation (Ahmdal's Law, also with p=0.95), however the parfeval version continues to show strange saturation behavior, even for the dummy function. Notice how the Speedup factor improves with core number upto 12 cores, then suddenly no further improvement is observed. I have attached the script in case you want to reproduce this behavior on your end.

Is there a fundamental limitation to the number of cores the parfeval function can leverage? Or is there an obvious mistake I am making in the way I am using the parfeval function? Why does the performance behavior of the dummy algorithm suddenly saturate at 12 cores? Any recommendation how to use the parfeval function to perform as well as parfor?

I would like to emphasize that I have already developed my application to use parfeval, so converting to parfor would be time-consuming and prevent me from utilizing the update-after-iteration feature of parfeval.

Thank you for your help on this critical matter.

[Rik]

I don't think I can personally help, but it seems relevant to know your number of physical and logical cores.

[Joseph Smalley]

Hi Rik,

There are 32 physcial cores and 64 threads. Here is the result from entering ver -support

-----------------------------------------------------------------------------------------------------

MATLAB Version: 9.8.0.1323502 (R2020a)

MATLAB License Number: ___

Operating System: Microsoft Windows Server 2019 Datacenter Version 10.0 (Build 17763)

Java Version: Java 1.8.0_202-b08 with Oracle Corporation Java HotSpot(TM) 64-Bit Server VM mixed mode

-----------------------------------------------------------------------------------------------------

MATLAB Version 9.8 (R2020a) License __

Image Processing Toolbox Version 11.1 (R2020a) License __

MATLAB Compiler Version 8.0 (R2020a) License __

Parallel Computing Toolbox Version 7.2 (R2020a) License __

[Rik]

I'm not sure what people could do with it, but I think I would redact that license number. I'm on mobile now, so it's a pain to edit it away for you.

[Joseph Smalley]

Thanks, done.

Answer 1

[Edric Ellis]

The main difference between parfor and parfeval is that in the parfeval case, you are responsible for scheduling the work on the workers. parfor has an advantage over parfeval in that it knows how many loop iterations there are, and so what it does is schedule a fixed number of chunks of work per worker (see the documentation for parforOptions - the chunks are referred to as "sub-ranges"). So, in your case, parfeval will incur more overhead since each parfeval request is sent on its own to a worker, where as parfor groups things together, and this will generally be more efficient in the case where the request durations are of a similar duration to the overheads of making a single remote request.

So, parfeval doesn't have a fundamental limitation in this regard, but you might need to amalgamate your requests if they are too short to match parfor performance. Another option might be to use parfor together with DataQueue which would let you perform updates at the client after each parfor iteration completes.

[Joseph Smalley]

Thank you for that explanation Edric. Based on your response, it seems that parfor would be the better choice. However, I have another issue in that the function that I call in each loop iteration takes two input arguments, each of which is a unique user-defined handle object. While parfeval updates the properties of these objects, use of parfor does not. Further I cannot use the reduction variable construct, because I need to update both user-defined handle objects. From my understanding, the reduction variable and additional arguments passed into a function must be associative. This would not be the case if the two variables are objects of different handle classes. Is there a workaround here? Returning to parfeval then, is there any advantage to use fetchOutputs rather than fetchNext?

[Edric Ellis]

I'm not sure I understand your problem here. Neither parfeval nor parfor can update handle objects at the client - handle objects always get copied to workers. Do you have a simple minimal reproduction for this that you can show?

[Joseph Smalley]

Below is the code snippet: W_temp and inputRays are array objects of two different user-defined handle classes. Effectively this algorithm updates their properties on each iteration. Using parfeval I am able to output both objects into the objListFuture variable on each iteration, then use their updated properties. Perhaps I may have used the wrong terminology in my previous comment, but this algorithm enables effective updating of both objects. It is unclear to me how I could achieve the same result via parfor, but based on your initial answer, it seems preferable to use parfor.

    objListFuture(1:N) = parallel.FevalFuture;
    W_temp(N,1) = W;
    for i=1:N
        W_temp(i) = W;
    end
    for i=1:N
        [objListFuture(i)] = parfeval(@IterTrace_oneParent_par,2,W_temp(i),inputRays(i));
        disp(['Iteration #', num2str(i), ' of ', num2str(N) ' is initialized.'])
    end
    
    objListNow           = cell(N,W.numObj);  
    maxSegPerRay         = W.maxSegments*W.maxBranches;
    rayListAll_origin    = zeros(3,N,maxSegPerRay);
    rayList_length       = zeros(N,1);
    
    for i=1:N
        % return index of fetch and object list
        [cindx,valOut_1,valOut_2]   = fetchNext(objListFuture);
        objListNow(cindx,:)         = valOut_1;
        rayList                     = valOut_2;    
        rayList_length(cindx)       = length(rayList);
        % convert updated property of inputRays object to numerical array
        rayListAll_origin(:,cindx,1:rayList_length(cindx))     = [rayList.origin];
        disp(['Iteration #', num2str(i), ' of ', num2str(N)  ' complete. (Ray #', num2str(cindx)]);

        % update property of W object using objListNow 
        for j=1:W.numObj    
            if class(W.objects{j})=="Detector"   
               if ~isempty(objListNow{cindx,j}.hitList)
                   W.objects{j}.hitList = objListNow{cindx,j}.hitList;
               end
            end
        end           
    end

[Edric Ellis]

From a brief inspection of that code, the main thing that looks un-parfor-able is the "update W" loop - that has order-dependent behaviour (or it appears to). It might be possible that you can recast that into an order-independent form. But the main computation ought to work something like this:

parfor i = 1:N
    [objListNow(i, :), rayList] = IterTrace_oneParent_par(W_temp(i), inputRays(i));
    rayList_length(i) = length(rayList);
    % to get the right "sliced" indexing form, need to be a little careful here
    rayList_tmp = zeros(1, maxSegPerRay);
    rayList_tmp(1:rayList_length(i)) = rayList;
    rayListAll_origin(:, i, :) = rayList_tmp;
end

I think the simplest thing would be to perform the "W" updates after the parfor loop, back on the client, hoping that that isn't the time-consuming piece of the operation.

[Joseph Smalley]

Edric, thank you for this snippet as it was helpful and I now have the parfor loop updating my objects. HOWEVER, unfortunately some unexpected behavior appears that does not appear when using parfeval or serialized for-loop. Namely, when the number of iterations of the parfor loop exceeds the number of workers in my paralell pool, the results deviate from the expected result, and the problem becomes worse as the number of iterations increases. Without going into too much application detail, the result becomes increasingly "unphysical" as the number of iterations exceeds the number of workers.

After some debugging I have found that because objList contains handle objects, the properties of these objects do appear to be updated on each "worker cycle", such that the objects contained in objList(7,:) are updates, rather than copies, of the objects contained in objList(1,:), assuming use of 6 workers. From your earlier comment, I thought this would not be the case.

Does this seem expected to you?

Is there anything in the new code (below) that appears as a red flag to you?

    %% parfor algorithm

    objList              = cell(N,W.numObj);
    maxSegPerRay         = W.maxSegments*W.maxBranches;
    rayListAll_origin    = zeros(3,N,maxSegPerRay);
    rayList_length       = zeros(N,1);

    parfor i=1:N
        [objList(i,:),rayList]      = IterTrace_oneParent_par(W,inputRays(i)); 
        rayList_length(i)           = length(rayList);
        zeroList_length             = maxSegPerRay - rayList_length(i);
        % convert updated ray objects into double arrays with consistent size; zeros are removed later
        rayListAll_origin(:,i,:)    = [rayList.origin, zeros(3,zeroList_length)];
    end
    
    % update property of detector objects contained in W object
    for i=1:N
        for j=1:W.numObj    
            if class(W.objects{j})=="Detector"   
               if ~isempty(objList{i,j}.hitList)
                   W.objects{j}.rays = objList{i,j}.rays;
               end
            end
        end
    end

[Edric Ellis]

As N exceeds the number of workers in the pool, the parfor machinery breaks the loop up into "subranges" to send to the workers. So, each worker will get a subrange of a number of loop iterations. objList here is a cell array, so I don't see how elements of that can interact with each other unless they're created that way. The other suspect here is the "broadcast" variable W. This will get copied once to each worker, and then the same instance will be used for multiple "subranges". So, if this has handle behaviour, that might explain it. Here's the sort of thing I'm thinking of.

h = containers.Map();
out = cell(1, 10);
parfor i = 1:10
    % Following line is to fool parfor into letting me modify 'h'
    hh = h;
    % Check the current contents of 'h'
    out{i} = hh.keys();
    % Modify 'h'
    hh(string(i)) = magic(i);
end

[Joseph Smalley]

Edric, after taking a break and returning to this problem, your last above recommended code still behaves in the non-physical way from my previous comment. The disp() line within my parfor loop (code below) checks a property of both the W_temp and W object. All values of "Src Intensity" should be equal to ~1. However after exceeding a multiple of the number of workers in my pool (6), the properties experience a "step-like" behavior in that one of the workers sees an updated object rather than the original object (before the parfor loop begins). Below shows Src Intensity jumping to 2 on the 7th worker and 3 on the 13th worker.

N=7

Iteration #3 of 7 complete. Src Intensity(temp)=1, Src Intensity(main)=1

Iteration #5 of 7 complete. Src Intensity(temp)=1, Src Intensity(main)=1

Iteration #4 of 7 complete. Src Intensity(temp)=1, Src Intensity(main)=1

Iteration #2 of 7 complete. Src Intensity(temp)=1, Src Intensity(main)=1

Iteration #1 of 7 complete. Src Intensity(temp)=1, Src Intensity(main)=1

Iteration #7 of 7 complete. Src Intensity(temp)=1, Src Intensity(main)=1

Iteration #6 of 7 complete. Src Intensity(temp)=2.0002, Src Intensity(main)=2.0002

----

N=13

Iteration #2 of 13 complete. Src Intensity(temp)=1, Src Intensity(main)=1

Iteration #3 of 13 complete. Src Intensity(temp)=1, Src Intensity(main)=1

Iteration #5 of 13 complete. Src Intensity(temp)=1, Src Intensity(main)=1

Iteration #4 of 13 complete. Src Intensity(temp)=1, Src Intensity(main)=1

Iteration #1 of 13 complete. Src Intensity(temp)=1, Src Intensity(main)=1

Iteration #6 of 13 complete. Src Intensity(temp)=1, Src Intensity(main)=1

Iteration #8 of 13 complete. Src Intensity(temp)=2, Src Intensity(main)=2

Iteration #11 of 13 complete. Src Intensity(temp)=2, Src Intensity(main)=2

Iteration #10 of 13 complete. Src Intensity(temp)=2, Src Intensity(main)=2

Iteration #9 of 13 complete. Src Intensity(temp)=2, Src Intensity(main)=2

Iteration #7 of 13 complete. Src Intensity(temp)=2, Src Intensity(main)=2

Iteration #13 of 13 complete. Src Intensity(temp)=2, Src Intensity(main)=2

Iteration #12 of 13 complete. Src Intensity(temp)=3.0001, Src Intensity(main)=3.0001

It then appears that W is updated, WITHIN the parfor loop, after a completed cycle of 6 workers. However upon completion of the parfor loop, only the W_temp object is updated. Hence I need to "manually" update the properties of W with a serial for loop, which is OK. The problem is that I do not want W_temp or W to be updated within the parfor loop after completion of a multiple of the size of the parallel pool. I want all workers to see the original W object for all iterations. Is this possible? Thank you for your continued assistance.

    % pre-allocation
    maxSegPerRay         = W.maxSegments*W.maxBranches;
    rayListAll_origin    = zeros(3,N,maxSegPerRay);
    rayList_length       = zeros(N,1);
    % W is the handle object whose properties include other handle classes that we want to update
    W_temp(N,1) = W;
    for i=1:N
        W_temp(i) = W;
    end

    parfor i=1:N
        % Convert broadcast variable into temporary variable
        W_temp2 = W;
        % Call main function
        [~,rayList]                 = IterTrace_oneParent_par(W_temp2,inputRays(i));
        % Convert updated rayList properties to numeric array (not a problem)
        rayList_length(i)           = length(rayList);
        zeroList_length             = maxSegPerRay - rayList_length(i);
        rayListAll_origin(:,i,:)    = [rayList.origin,      zeros(3,zeroList_length)];
        % Update W_temp object and display Src Intensity of W_temp and W (should always be ~1)
        W_temp(i) = W_temp2;
        disp(['Iteration #', num2str(i), ' of ', num2str(N)  ' complete.  Src Intensity(temp)=', num2str(sum([W_temp(i).objects{5}.rays.intensity])),...
            ',  Src Intensity(main)=', num2str(sum([W.objects{5}.rays.intensity]))]);
    end
    % Note: W_temp and W are both updated WITHIN the parfor loop after subRange is complete, but only W_temp is updated on completion of the parfor loop

    % "Manually" update properties of detector objects contained in W
     for i=1:N
        for j=1:W.numObj    
            if class(W.objects{j})=="Detector"                  
               if ~isempty(W_temp(i).objects{j}.rays)
                   W.objects{j}.rays    = W_temp(i).objects{j}.rays;
               end
            end
        end
     end
    

[Joseph Smalley]

Just wanted to say that I accepted this answer because, overall, the problem is addressed more easily by switching to a parfor loop, as Edric first proposed. Additionally I switched all my classes to value classes over handle classes. The latter is a compromise for my application, and was first motivated by requirements of codegen for MEX files. Nonetheless the combination of parfor with values classes has been working for several weeks now, with pretty good scalability of 12x at 24 cores.