improving the speed of parallel optimization

Hi, I am trying to optimize in parallel but the speed is increased just slightly by using parfor. Do you have any further recommendations? Thanks!

    parfor i = 1:M

    options = optimset('MaxFunEvals',Inf,'MaxIter',10,...
    'Algorithm','interior-point','Display','iter'); 
    startTime = tic; 
    [x(:,i),fval(:,i)] = fmincon(@(x)revenue(price(1:N,1),ro,g,eff,x,N,k1,init_s(i),inFlow(:,i),alpha_par(i),b_par(i)),x0(1:2*N,i),A,b(:,i),Aeq,beq(:,i),LB(:,i),UB(:,i),[],options);      
    time_fmincon_parallel = toc(startTime);
    fprintf('Parallel FMINCON optimization takes %g seconds.\n',time_fmincon_parallel);
    end

Answer 1

[Walter Roberson]

Instead of running the fmincon calls in parallel, try running them in a loop, but using the option UseParallel to allow parallel estimation of the gradient.

Remember, it is common for Parallel processing to be slower than serial, depending on the amount of data to be transferred compared to the amount of work to be done per iteration, and taking into account that non-parallel workers can use the built-in parallelization of some operations on large "enough" matrices by calling into LaPACK / MKL.

[sensation]

Thanks for your comment Walter. I run it with parfor; useparaller and no parallel and here are results for reduced problem [10,1]:

parfor: 3022 sec

useparallel: 2195 sec

no parallel: 2629 sec

Still however, takes too long taking into account my matrix that is [4000,170] in fmincon. Do you have any other suggestion how to speed it up please?

Thanks and cheers!

[Walter Roberson]

Generalities:

• if you can rewrite the problem as a linear or quadratic problem, you probably should do that and use an appropriate solver, as there are analytic methods to deal with those
• nonlinear constraints are expensive; avoid them if you can. Fortunately for you, you do not have any
• linear equality constraints are slower than linear inequalities
• When you have a linear equality, examine your code. You might be able to reduce by one or more variables, rewriting the remaining variable internally in the objective function as (constant minus linear sum of the remaining variables). Reducing the number of variables also makes it faster to estimate gradient and hessian as the cost of estimating hessians grows with the square of the number of variables
• Sometimes there are regions where certain variables have no effect on the outcome. It can sometimes make sense to break up the problem into sub-problems with different objective functions that deal only with the regions where the variables are meaningful
• never load() inside your objective function; load first and pass the data in as part of the parameters to the objective function
• Give some thought to vectorizing your objective function. Your objective function will be executed many many times, so you want it to be efficient.
• if you have an efficient function for calculating gradients, then provide it in the options instead of having it estimate the gradients and hessian. But do keep in mind that having an analytic function to calculate the gradients or hessian does not always mean that you should calculate it analytically: getting all of the details right for exact solutions can be computationally expensive. The gradient and hessian estimation that is done is sort of like doing a truncated taylor expansion: because of the truncation it could end being much faster, but on the other hand the effective truncation can lead to some pretty misleading results. If your objective function is somewhat smooth at a macro level then gradient estimation can be enough, but if it has some pretty steep areas then you might need to provide a more analytic estimator.

[sensation]

Thanks for your insight Walter. Do you know how I can thus calculate hessian and grad of my function above knowing that:

      function R=revenue(price,ro,g,eff,x,N,k1,init_s,inFlow,alpha_par,b_par)
             HH=depth(alpha_par,b_par,init_s,inFlow,x,N); 
             R= -sum(price*((HH*ro*g*eff)*x(1:N))/k1);
      end

Thanks and cheers!

[Walter Roberson]

We would need to know depth() to be able to answer.

If you have the symbolic toolbox, then what can help is to create a symbolic vector the same of appropriate length and pass it into your objective function, hopefully getting back a symbolic expression for your objective function. Sometimes the result is just ugly messy, but other times the interaction with the inputs turns out to be fairly clean and you can use diff() and gradient() to come up with reasonable functions and then matlabFunction() to create appropriate function handles.

One trick with matlabFunction is to consider using the 'File' output option, which by default enables the optimization option. The symbolic engine then spends a bunch of time doing common sub-expression optimization to come out with a more efficient series of statements to do the calculation. But this process can take a lot of time; if the expression is longer and more complicated then it can take more than an hour, so it is only worth doing if the expression is not complicated or if the objective function will be executed a lot of times.

[sensation]

Thank you for your answer Walter. My other functions are:

       function HH=depth(alpha_par,b_par,init_s,inFlow,x,N)
             S=storage(init_s,inFlow,x,N);
             HH= (S/alpha_par).^(1/b_par)*0.8; 
       end

then S:

    function S=storage(init_s,inFlow,x,N)
         T=totalflow(x,N);

         S(1) = init_s(1) + inFlow(1)-T(1);  % x(1:N)+x(N+1:2*N)

        for ii = 2:N  
                S(ii) = S(ii-1) + inFlow(ii)-T(ii);

        end
   end

and T:

       function T=totalflow(x,N)
            T=x(1:N)+x(N+1:2*N);      
       end

with x being the one I want to optimize. I have symbolic toolbox, and beforehand I had all code written as sym but my simulations never ended. So thats why, I choose to go along without syms and it goes faster now. Do you suggest to use sym maybe only in the last function so I can use diff and grad as you suggested on that? Thank for any hints really. Best!

[Walter Roberson]

About how large is N?

[sensation]

[4000,200] =[N,M] Walter. Thanks!

[sensation]

Any idea Walter? Thanks!

[Walter Roberson]

If I recall correctly, with N even close to that large, asking matlabFunction to optimize the code takes far far too long, so I do not think you are going to be able to take advantage of that.

Answer 2

[Matt J]

This is a more optimal implementation of storage(),

    function S=storage(init_s,inFlow,x,N)

         D=inFlow-totalflow(x,N);

         D(1) = D(1) + ( init_s(1) + D(1) );  

         S=cumsum(D);

   end

[Matt J]

What are the dimensions of all the quantities in your objective function

R= -sum(price*((HH*ro*g*eff)*x(1:N))/k1);

In particular, are any of them scalars?

[sensation]

Hi Matt, they are all of the size (4000,200) meaning 4000 timesteps and 200 case studies. ro,g,eff,k1 are constants; price is [4000,1]; HH calculated [4000,200].

[Matt J]

ro,g,eff,k1 are constants

Are they scalar constants?

[Matt J]

HH calculated [4000,200].

It doesn't look possible for HH to be 4000x200. Your depth() function returns an N-vector, not an NxM matrix.

[sensation]

Are they scalar constants? Yes

[sensation]

in:

        function R=revenue(price,ro,g,eff,x,N,k1,init_s,inFlow,alpha_par,b_par)
             HH=depth(alpha_par,b_par,init_s,inFlow,x,N); 
             R= -sum(price*((HH*ro*g*eff)*x(1:N))/k1);
        end

so it is called in optimization for each M and all N and so on. At the end I meant if I derive HH values the vector would be [N,M] but ofcourse the calculations are done for each case study in optimization separately. Thanks!

[Matt J]

At the end I meant if I derive HH values the vector would be [N,M] but ofcourse the calculations are done for each case study in optimization separately.

OK, then for clarity please say again what the dimensions are of all of the variables in the calculation

   R= -sum(price*((HH*ro*g*eff)*x(1:N))/k1);

as they exist in the workspace of revenue() when it is executed. It looks like HH has to be a 1xN row vector, not an Nx1 column vector, based on the calculations I see in storage().

[sensation]

Exactly, they are all 1D vector for each case study Matt and then, they are passed to revenue as such.

[Matt J]

OK, so HH is a 1xN row , x(1:N) is an Nx1 column, What about price? Row or column?

And what about ro,g,eff,k1. Earlier, you said they are scalars, not vectors. Is that still true?

[sensation]

Price is as X (one column) Matt. The latter are all scalars yes. Thank you very much

[Matt J]

If that's the case, then your objective function only depends on 'price' through its sum, which you shoud pre-compute before calling fmincon. Your objective computation then simplifies to

sumprice=-(ro*g*eff/k1)*sum(price);

function R=revenue(x,N,sumprice,init_s,inFlow,alpha_par,b_par)

      HH=depth(alpha_par,b_par,init_s,inFlow,x,N); 

      R= sumprice*(HH*x(1:N));

 end

So now you have 2 ways of making the computation more efficient: the above modification of your objective function and the re-implementation of storage() in my initial post.

[sensation]

Thank you a lot Matt, I have implemented those suggestions and its a bit faster but still too slow:(. Do you know maybe how I can run it with quadprog instead of fmincon? how to calculate f and H? Thanks a lot! really appreciate your help!

[Matt J]

I have implemented those suggestions and its a bit faster

How fast is it now? It should have been a lot faster than what you were doing.

Do you know maybe how I can run it with quadprog instead of fmincon?

The problem doesn't look quadratic, except maybe when b_par=1.

[sensation]

With N being 4000 time steps, the optimization is quite long. It is a bit faster 1.2 times.

I was thinking to run fmincon with supplied hessian and gradient maybe. I run it in parallel with parfor but still not efficient.

If you have any idea to make it feasible, would appreciate it a lot Matt. Thanks!