I am looking to accelerate the computation of Input-Output Polynomial Models by using parfor.

I have a function that uses a for loop, but it takes a lot of time to complete the model identification.

function [Idt_OE] = ident_oe(tt,nk)
Idt_OE=[];
count = 4^6;
for k = 1:count
    [nb1, nb2, nb3, nf1, nf2, nf3] = ind2sub([4, 4, 4, 4, 4, 4], k);
    Model_OE = oe(tt, [[nb1 nb2 nb3] [nf1 nf2 nf3] [nk nk nk]]);
    [a,b] = compare(tt,Model_OE);
    Idt_OE = [Idt_OE; b, nb1, nb2, nb3, nf1, nf2, nf3];
end

Idt_OE = sortrows(Idt_OE, 1);
end

And I tried using the code below, but something is wrong because the fit is not calculated very well.

 function [Idt_OE] = ident_oe(tt, nk)
     Idt_OE = [];
     count = 4^6; 
     Idt_OE_par = cell(count, 1); 
     parfor k = 1:count
         [nb1, nb2, nb3, nf1, nf2, nf3] = ind2sub([2, 2, 2, 2, 2, 2], k);
         try
             Model_OE = oe(tt, [[nb1 nb2 nb3] [nf1 nf2 nf3] [nk nk nk]]);
             [~, fit] = compare(tt, Model_OE);
             Idt_OE_par = [fit, nb1, nb2, nb3, nf1, nf2, nf3];
         catch
             Idt_OE_par{k} = [-Inf, nb1, nb2, nb3, nf1, nf2, nf3];
         end
     end
     Idt_OE = cell2mat(Idt_OE_par);
     Idt_OE = sortrows(Idt_OE, 1, 'descend');
 end
 

Many thanks,

[Walter Roberson]

Part of the slowness is that you are growing Idt_OE each iteration. Each time you grow the array, the previous data needs to be copied. The first iteration, 0 rows are copied. The second iteration, 1 row is copied -- total 1 row so far. The second iteration, 2 rows are copied -- total 3 rows so far. The third iteration, 3 rows are copied -- total 6 rows so far. The N'th interation, N*(N-1)/2 total rows will have been copied.

You can eliminate the copying by using

%Idt_OE=[];    %don't bother to initialize
count = 4^6;
for k = count:-1:1
    [nb1, nb2, nb3, nf1, nf2, nf3] = ind2sub([4, 4, 4, 4, 4, 4], k);
    Model_OE = oe(tt, [[nb1 nb2 nb3] [nf1 nf2 nf3] [nk nk nk]]);
    [a,b] = compare(tt,Model_OE);
    Idt_OE(k,:) = [b, nb1, nb2, nb3, nf1, nf2, nf3];
end

[Gabriel]

Thank you very much for the response. My problem is with the second code, when I use parfor. For some reason, the fit values with for and parfor are significantly different.

My problem is when i use the code below.

Many Thanks,

function [Idt_OE] = ident_oe(tt, nk)
     Idt_OE = [];
     count = 4^6; 
     Idt_OE_par = cell(count, 1); 
     parfor k = 1:count
         [nb1, nb2, nb3, nf1, nf2, nf3] = ind2sub([2, 2, 2, 2, 2, 2], k);
         try
             Model_OE = oe(tt, [[nb1 nb2 nb3] [nf1 nf2 nf3] [nk nk nk]]);
             [~, fit] = compare(tt, Model_OE);
             Idt_OE_par = [fit, nb1, nb2, nb3, nf1, nf2, nf3];
         catch
             Idt_OE_par{k} = [-Inf, nb1, nb2, nb3, nf1, nf2, nf3];
         end
     end
     Idt_OE = cell2mat(Idt_OE_par);
     Idt_OE = sortrows(Idt_OE, 1, 'descend');
 end

[Walter Roberson]

https://www.mathworks.com/help/ident/ref/oe.html Report property includes a RandState field, "State of the random number stream at the start of estimation. Empty, [], if randomization was not used during estimation. For more information, see rng."

This implies that randomness is typically used during estimation. The random state for parfor is different for parfor; see https://www.mathworks.com/help/parallel-computing/control-random-number-streams-on-workers.html

It is possible that by chance the seeds used for parfor happen to lead to worse results than the seeds used for serial processing.

[Gabriel]

The problem persists; even after fixing the seed, the difference between the calculations remains.

[Walter Roberson]

I am not sure why you have ind2sub([2, 2, 2, 2, 2, 2], k) in the parfor version, but ind2sub([4, 4, 4, 4, 4, 4], k) in the serial version? You are going to run into errors when k > 2^6 in the ind2sub

[Gabriel]

I did this solely to minimize interactions within the pafor, allowing for quicker verification during test runs of the solution.

Many Thanks,

[Walter Roberson]

parfor generally starts from the end of the range. So k=4^6 is going to be processed almost immediately. But that is going to generate

Error using sub2ind (line 71)

Out of range subscript.

soon after parfor starts.

[Gabriel]

Sory, I undestand now, yes we have k=4^6 to k=2^6 in this case.

[Walter Roberson]

As an experiment, try

Idt_OE_nonpar = cell(1,1);     
Idt_OE_par = cell(1, 1);

k = 1;
     [nb1, nb2, nb3, nf1, nf2, nf3] = ind2sub([2, 2, 2, 2, 2, 2], k);
     try
             rng(655321);
             Model_OE = oe(tt, [[nb1 nb2 nb3] [nf1 nf2 nf3] [nk nk nk]]);
             [~, fit] = compare(tt, Model_OE);
             Idt_OE_nonpar{k} = [fit, nb1, nb2, nb3, nf1, nf2, nf3];
     catch
             Idt_OE_nonpar{k} = [-Inf, nb1, nb2, nb3, nf1, nf2, nf3];
     end
     
 parfor k = 1
     nb1, nb2, nb3, nf1, nf2, nf3] = ind2sub([2, 2, 2, 2, 2, 2], k);
     try
             rng(655321);
             Model_OE = oe(tt, [[nb1 nb2 nb3] [nf1 nf2 nf3] [nk nk nk]]);
             [~, fit] = compare(tt, Model_OE);
             Idt_OE_par{k} = [fit, nb1, nb2, nb3, nf1, nf2, nf3];
     catch
             Idt_OE_par{k} = [-Inf, nb1, nb2, nb3, nf1, nf2, nf3];
     end

end

and see whether Idt_OE_nonpar and Idt_OE_par come out the same.

[Walter Roberson]

             Idt_OE_par = [fit, nb1, nb2, nb3, nf1, nf2, nf3];

Note that you are overwriting all of Idt_OE_par in that case. You should be assigning to Idt_OE_par{k}

[Gabriel]

I made this code, and the difference is still very big.

Many Thanks.

[Walter Roberson]

It is a mystery to me how parfor could be returning different values in that situation.

Hmmm... as an experiment try

Idt_OE_nonpar = cell(2,1);     
k = 1;
     [nb1, nb2, nb3, nf1, nf2, nf3] = ind2sub([2, 2, 2, 2, 2, 2], k);
     try
             rng(655321);
             Model_OE = oe(tt, [[nb1 nb2 nb3] [nf1 nf2 nf3] [nk nk nk]]);
             [~, fit] = compare(tt, Model_OE);
             Idt_OE_nonpar{k} = [fit, nb1, nb2, nb3, nf1, nf2, nf3];
     catch
             Idt_OE_nonpar{k} = [-Inf, nb1, nb2, nb3, nf1, nf2, nf3];
     end
     try
             rng(655321);
             Model_OE = oe(tt, [[nb1 nb2 nb3] [nf1 nf2 nf3] [nk nk nk]]);
             [~, fit] = compare(tt, Model_OE);
             Idt_OE_nonpar{k+1} = [fit, nb1, nb2, nb3, nf1, nf2, nf3];
     catch
             Idt_OE_nonpar{k+1} = [-Inf, nb1, nb2, nb3, nf1, nf2, nf3];
     end

This should return identical results.