Use Parfoor Loop for Parameter Sweep Optimization

Hello, I have a general optimization problem that the solution looks something like this (sweeping across parameter p):

maxval = 0;
for a = 1:length(p);
    val = somefunction(p(a));
    if val > maxval
        pmax = p(a);
        maxval = val; 
    end
end

Where I would like to convert into a Parfor loop to find the optimal p that maximizes val. I am struggling to avoid converting into a parfor loop by using only reduction variables. The best I could do looks something like this:

parfor a = 1:length(p);
    listval = [val;somefunction(p(a))];
    listp = [listp,p(a)];
end

imax = find(val == max(val));
pmax = listp(imax);

So as you can see I am storing the value of every iteration when I am only interested in finding the max value. Also, I am not preallocating memory and my lists (listval & listp) are increasing in size in every iteration. I have the intuition there must be a better solution, as this seems very inefficient, but I have not been able to find anything on the documentation. I have also tried:

maxval = max(maxval,val) instead of the if statement

but of course this only yields the value for the best solution but not the parameter that achieves it. I would appreciate any pointers.

Answer 1

[Edric Ellis]

You can do this using a custom reduction function in parfor. Like this:

p = 1:10;
someFcn = @sin;
% Use a 2-element vector to store maximum
% value of someFcn(x), as well as the input x.
maxVal = [-Inf, NaN];
parfor a = 1:length(p)
    val = someFcn(p(a));
    % Next, use our "custom reduction" to 
    % update maxVal if the value piece is larger
    % than previous values.
    maxVal = iMax(maxVal, [val, p(a)]);
end

Starting parallel pool (parpool) using the 'Processes' profile ...
Connected to the parallel pool (number of workers: 2).
Analyzing and transferring files to the workers ...done.

disp(maxVal)

    0.9894    8.0000

function v = iMax(v, in)
if in(1) > v(1)
    v = in;
end
end

(Only after writing out this solution did I notice that the doc example for for custom reduction is precisely this case!)

[Andres Morales]

Thank you. This is exactly what I was looking for, since I have no interest in saving anything else besides the optimal solution.

Answer 2

[Walter Roberson]

List concatenation is one of the permitted reduction strategies, so you can assume that it will handle memory allocation appropriately.

But you could also use an array, something like

N = 1000000;
pvals = zeros(2,N);

parfor a = 1 : N
    thisp = p(a);
    val = somefunction(thisp);
    pvals(:,a) = [val;thisp];
end
[maxval, idx] = max(pvals(1,:));
pmax = pvals(2,idx);

Question: your current code would find all the places that have the same (identical) max. Is that deliberate? (There can be good reason to want to know them all, but it does mean there can be instabilities that are due to floating point round off. Though I guess the instabilities exist if you just want one out of all the locations that are maximal.)

[Andres Morales]

Yes, this will work, thank you! I tried using indexing on a preallocated array, but didn't think of using a for the indexing. The use of it as a sliced variable is clever!

Yeah, I didn't consider the possibility of having identical maximums. In my application, knowing one is sufficient, but that is a good point.

P.S. - I edited my question to show a = 1:length(p) in both code snippets as it was not consistent. Thus, what I need is N = length(p).

_____________________________

P.S.S. - This is my first question, and I am trying to find in the community what is the etiquette for asking side questions and/or follow up questions. If this is not the place, let me know and I will fix it.

My side question is why is [maxval, idx] = max(pvals(1,:)) preferred over imax = find(val == max(val)) ?

I have a follow-up question but that one could have enough conversation to be a standalone post. I will link this question in there.

[Walter Roberson]

[maxval, idx] = max(pvals(1,:)) only involves a single pass over the data, with it internally keeping track of the location of the maximum as it goes.

find(val == max(val)) involves a pass over the data to find the maximum, then a second pass over the data to compare each element to the maximum (producing a logical vector), and then a pass over the logical vector to locate all of the indices of the matches. It is therefore roughly 3 times as expensive compared to the single max() call.

A reason to use code you had would be if you needed the indices of all of the locations where the value matched the maximum.

pmax = listp(imax);

Note that in your original code the only thing you did with the indices returned by find(), was to index listp -- you did not display the indices or use them for other processing.

In that situation, where you want to know the locations so that you can index another array the same size, then it is less expensive to do logical indexing:

imax = val == max(val);
pmax = listp(imax);

This saves one of the three passes over the data, which is one factor -- but it also saves some internal calculations that compare each index to the subscripts (and rememember that in general the subscripts need not be sorted and need not be unique.) So even though with logical indexing there must surely be an internal pass to count the number of matches to know how large of an output will be needed, logical indexing is still faster than using explicit subscripts.

[Andres Morales]

Yeah, that's very clear. Thank you for providing great answers.