How to make a multiple for-loop run faster, the function of which is to find the steady state?

Hi. I am programing a double for-loop to find the steady state of a iterative physical process (include fft/ifft of matrix). The program works well with default single thread calculation but takes to much time. Whereas the 'parfor' I've tried seems to be useless here, and even spent more time. I attach a part of the codes below. Much appreciate if anyone could give me advice to help improve the speed of my codes, or a better way to find the steady state.

%% input parameters and create variables
nt = 2^12;                            % sampling points
N = 2000;                             % max number of iteration
A = zeros(N,nt);                      % create initial A
A(1,:) = randn(1,nt)+1i*randn(1,nt);  
% ...and other parameters

%% Main loop
rmse = ones(N,1);                     % Create Root-mean-square error variable
B = fftshift(fft(ifftshift(A)));      
C = abs(B).^2;
flag = 0;                             % initial steady-state mark == 0

for m = 1:10 % to find the minimum m value which leads to steady state (linear, when m>=threshold, become steady)
    for n = 1:N-1                     % iterative process
        A(n+1,:) = Function(A(n,:),m);% a custom function, input A(n) and m, output A(n+1)
        B(n+1,:) = fftshift(fft(ifftshift(A(n+1,:))));
        C(n+1,:) = abs(B(n+1,:)).^2;
        rmse(n+1) = sqrt(mean((C(n+1,:)-C(n,:)).^2)); % calculate Root-mean-square error between C(n) and C(n+1)
        fin = n+1;                                    % store n value
        if (rmse(fin) < 0.02) && (fin > 2)            % check out steady state
            flag = 1;                                 % set mark == 1 when at steady state
            break
        end
    end
    if flag == 1                                      % junmp out the loop
        break
    end
end

%% Visualization
figure;
subplot(2,1,1);
mesh(A);
subplot(2,1,2);
mesh(C);

Answer 1

[Raymond Norris]

Yuu ought to be able to change the last several lines as such

    end
    temp(m,1) = rmsep;
end
rmse = temp;

With that said, I'm not going to go line by line, other than to say, I don't think what you originally posted maps to your parfor suggestion. Specifically, the original example might exit early, whereas your parfor example will not (well, only partially).

Answer 2

[Raymond Norris]

Can you elaborate a bit on the following:

• how are you measuring the time (walltime, tic/toc, profiler)
• how long it took to run serially and what speedup you're looking for (e.g. takes 2 days, lookng for 3 hours, etc.)
• what compute do you have access to that you're running your code on (e.g. 4 cores, 16 GB RAM, etc.)
• I'm assuming you're trying to parallelize the outter for-loop since rmse has a dependency on C(n+1) and C(n). Regardless of which one you chose, you'd need to re-write the for-loop as they're currently written. Can you post what your parfor loops look like.

steady-state problems (i.e. "needle in the haystack") won't work with parfor, where you can't preemptively end the loop. Instead, you're better off with parfeval, where you can see the results and cancel "futures" when you've met a threashold.

For instance

% Set the size of the pool how you choose (can be less then the number of
% iterations
p = parpool(10);
for m = 10:-1:1
    f(m) = p.parfeval(@calc_steady_state,3,m);
end

for m = 1:10
    [indx, flag, A, C] = f(m).fetchNext;
    if flag==true
        % Found our steady-state early
        f.cancel
    end
end

%% Visualization
figure;
subplot(2,1,1);
mesh(A);
subplot(2,1,2);
mesh(C);


function [flag, A, C] = calc_steady_state
%% input parameters and create variables
nt = 2^12;                            % sampling points
N = 2000;                             % max number of iteration
A = zeros(N,nt);                      % create initial A
A(1,:) = randn(1,nt)+1i*randn(1,nt);
% ...and other parameters
%% Main loop
rmse = ones(N,1);                     % Create Root-mean-square error variable
B = fftshift(fft(ifftshift(A)));
C = abs(B).^2;
flag = 0;

for n = 1:N-1                     % iterative process
    A(n+1,:) = Function(A(n,:),m);% a custom function, input A(n) and m, output A(n+1)
    B(n+1,:) = fftshift(fft(ifftshift(A(n+1,:))));
    C(n+1,:) = abs(B(n+1,:)).^2;
    rmse(n+1) = sqrt(mean((C(n+1,:)-C(n,:)).^2)); % calculate Root-mean-square error between C(n) and C(n+1)
    fin = n+1;                                    % store n value
    if (rmse(fin) < 0.02) && (fin > 2)            % check out steady state
        flag = 1;
        break
    end
end
end

[ChiQAQ]

Hi Norris. Thanks for your quick and professional answer. I will try parfeval as the example you demonstrated. Frankly speaking I'm not very familiar with parallel pool function. Let me learn more about it via Help Center first.

FYI,

• Use profiler to measure time.
• It takes 4 min to complete one internal for-loop consists of 2000 iterations. So the whole program would cost 40 min, counting the external loop. However, this is a simplified test code. The one I am handling includes at least m=1:380, and even need triple for-loop to address more parameters. Just considering the case of m=1:10, I prefer no more than 10 min, if possible.
• i7-6700@3.4GHz, 4 cores, 8 threads, 16GB RAM. (Also have GPU but didn't use it)
• You're right. I've tried to parallelize the outter loop. Below is the parfor loop I imitate from another case.

% create A, B, C, rmse, fin
A = zeros(N,nt,10);
A(1,:,:) = repmat(randn(1,nt)+1i*randn(1,nt),1,1,10);
B = fftshift(fft(ifftshift(A)));
C = abs(B).^2;
rmse = ones(10,1);
fin = zeros(10,1);

% parfor loop
parfor m = 1:10
    Ap = A;  % Cannot directly use A, cuz A is not recognized as sliced variable.
             % Or error happens "Valid indices for 'A' are restricted in PARFOR loops"
    Bp = B;  % B, C, rmse, fin have same reason
    Cp = C;
    rmsep = rmse;
    finp = fin;
    for n = 1:N-1
        Ap(n+1,:,m) = Function(Ap(n,:,m),m);
        Bp(n+1,:,m) = fftshift(fft(ifftshift(Ap(n+1,:,m))));
        Cp(n+1,:,m) = abs(Bp(n+1,:,m)).^2;
        rmsep(m) = sqrt(mean((Cp(n+1,:,m)-Cp(n,:,m)).^2));
        finp(m) = n+1;
        if (rmsep(m) < 0.02) && (finp > 2)
            break
        end
    end
    temp{m} = rmsep;
end

for m=1:10
    rmse(m) = temp{m};
end