Using interp1 multiple times to find a root : interp1 on a matrix without loops. Using interp1 multiple times for different values of query points without loops.

Dear All, I have the following problem:

I have a vector of arguments x and a vector of values v (both are monotonic; the root may not necessarily exist). I would like to find a root of function v(x). I could flip the axes and use interp1 in point 0:

root = interp1(-v,x,0);

The problem is that I have to do it many times, i.e. for many values of parameters affecting the vector v, which is actually a vector v(alpha,beta,gamma). Obviosuly I could use nested loops:

for alpha = 1:a
for beta = 1:b
for gamma = 1:c

root(alpha,beta,gamma) = interp1(-v(alpha,beta,gamma),x,0); 

end
end
end

Unfortunately, this solution is very slow (a,b and c are large numbers). As far as I know interp1 does not allow to use a matrix as its first input so I do not know how to vectorize it.

How to solve this problem in a fast way? Either by speeding up the solution using vectors getting rid of loops or by using any other method than interp1. Interp1 has quite a large fixed cost of being called so it is very slow in the loop environment while it speeds up when vectorized (but as I know only the second argument (i.e. values ) can be used as a mtrix )

The other problem is:

Is it possible to use interp1 for multiple values of query points without using loops? something like this:

root = interp1(x,v(:,alpha,beta,gamma),x0(alpha,beta,gamma));

where x is a vector, v(:,alpha,beta,gamma) is a multidimensional matrix (array) and x0(alpha,beta,gamma) is a point. The given syntax is doing a 1 dimensional interpolation over the first dimension of v for every combination of alpha,beta and gamma for a given point x0, which varies for all combinations of parameters.

I would be grateful for any hints!

[Image Analyst]

How large are a, b, and c? And why do you need to make up this huge 3D matrix? Give us some background/context.

[Adam Pigon]

Actually there are more parameters/variables than just alpha, beta and gamma that were used just for the exposition of the issue. Unfortunately, the curse of dimensionality is quite a big problem in this case. Ideally the dimension should be like 50x50x50x2x3x60. Loops are simply prohibitively slow. I have already vectorized some of the code (incredible speed-up!) and the values of vector v used in the question above are just values of a function whose arguments are multidimensional matrices.

general background: this problem is a way (one of the steps) to circumvent the curse of dimensionality in the dynamic programming problem.

Answer 1

[Daniel kiracofe]

My first question here would be, are you pre-allocating root? i.e. do you have root = zeros(a,b,c) before your loop? If not, then probably what is taking up the time is memory allocation, not the actual interpolation.

If memory is not your problem, then you should be able to get some speedup by using a modified version of interp1(). i.e. interp1() is a pretty general function. it works for either increasing or decreasing vectors, it checks to make sure you don't have duplicated points, it allows multiple interpolation method, etc. If you know for sure that you only have, say, increasing vectors and never have duplicated points, you could take advantage of that knowledge to write a quicker routine. Or take advantage of such routines already written:

https://www.mathworks.com/matlabcentral/fileexchange/43325-quicker-1d-linear-interpolation--interp1qr

[Adam Pigon]

Thank you for your comments! In fact your suggestion of resigning from pre-allocation was very useful - it improved the performance by about 20-25% ! My question is then, why it was that useful? I have always read the all the matrices should be pre-allocated as matlab consumes a lot of time when creating them again and again in loops. Was I missing something?

I should apply the other suggestion about other interpolation routines - hopefully it will be useful as well.

I have not found any methods how to vectorize interp1 and I think this is not possible so using any hints like the ones above is the best solution.

[Daniel kiracofe]

Wait, you removed the pre-allocation and your code got 20% faster? I was suggesting that you should add pre-allocation. So it sounds like you are getting the exact opposite result of what I expected. Pre-allocation should be speeding up your code. If it's not, then something weird might be going on. Suggest you post all of your code (or at least a larger portion) so we can take a look.

Also, another possible speedup: pay attention to row-major order vs column major order when setting up your loops (https://en.wikipedia.org/wiki/Row-_and_column-major_order#Column-major_order). You might already be doing it the fastest way, but since we don't have your full code I'm not sure and I thought I would point it out. You might get a little speedup by changing the order in which you are traversing the array. Might not be huge, but if it's a big problem every bit helps.

Finally, if speed really is a big problem, you might just consider re-writing in fortran. Fortran is horribly painful to program in, compared to matlab, but it will be faster. How much faster is dependent on your exact code. Could be just a few percent, but could be an order of magnitude.

[Adam Pigon]

well... it seems that I haven't read your comment careful enough and I got rid of preallocation but I can confirm that it leads to this strange result - removing preallocation speeds up the code by these 15-20% (which is visible for larger number of points, for small sets of parameters it is negligible).

this is my code:

tic

for t = (T-1):-1:1  % loop over time periods t

    for q = 1:nx    % auxilliary loop for cases 1 and 2, necessary condition: na = nx

    % case 1,  
    AP1  = a(q);
    F1(q,t,:,:,:,:,:) = f(AP1,t,agrid,hgrid,xgrid,Agrid,Sgrid);

    % case 2,
    AP2  = xp(q)/(1+r-deltaA);    
    F2(q,t,:,:,:,:,:) = beta .* ( E_util_ap(AP2,t,agrid,hgrid,xgrid,Agrid,Sgrid) + ( r - deltaA - adj1(AP2,agrid) - E_adj2(AP2,t,agrid,hgrid,xgrid,Agrid,Sgrid) ) .* E_util_cp(AP2,t,agrid,hgrid,xgrid,Agrid,Sgrid) ) ./ ( 1 + adj1(AP2,agrid) );

    end             % end of the auxilliary loop

    % case 1, eta = 0, mu = 0, a' > 0, m' > 0 

    for aux1 = 1:na
        for aux2 = 1:nh
            for aux3 = 1:nx
                for aux4 = 1:nA
                    for aux5 = 1:nS
                        ap1(:,:,:,:,:) = interp1(-squeeze(F1(:,t,aux1,aux2,aux3,aux4,aux5)),a,0,'spline');
                    end
                end
            end
        end
    end

    % case 2, eta > 0, mu = 0, a' = x'/(1+r-deltaA), m' = 0
    for p = 1:nx
    eta2(:,:,p,:,:) = F2(p,t,:,:,p,:,:);
    end

    % case 3, eta = 0, mu > 0, a' = 0, m' = x'
    AP3 = 0;
    mu3(:,:,:,:,:) = - beta .* ( E_util_ap(AP3,t,agrid,hgrid,xgrid,Agrid,Sgrid) + ( r - deltaA - adj1(AP3,agrid) - E_adj2(AP3,t,agrid,hgrid,xgrid,Agrid,Sgrid) ) .* E_util_cp(AP3,t,agrid,hgrid,xgrid,Agrid,Sgrid) ) ;

end                   % end of the T loop   

toc

The t loop is not working yet but it it immaterial to this problem. The preallocation is like this:

ap1  = zeros(na,nh,nx,nA,nS);

Functions like f or E_util_cp have been vectorized and their arguments are multidimensional matrices like agrid, hgrid etc that were created by ndgrid command.

The functions are like this ('ap' is the scalar argument as I could not vectorize the interp1 function):

util_cp = @(ap,t,da,dhp,dxp,dA,dSp) ( xp(dxp) + A(dA) .* hp(dhp) .* Sp(dSp) - repmat( interp1( a, squeeze(efinal(1,:,:,:,:,1)), ap, 'spline'), na, 1,1,1,nS ) - repmat( interp1( a, squeeze(afinal(t+1,:,:,:,:,dSp)), ap, 'spline' ), na, 1,1,1,nS ) - repmat( interp1( a, squeeze(mfinal(t+1,:,:,:,:,dSp)), ap, 'spline' ), na, 1,1,1,nS ) - repmat( lambda/omega .* ( log( exp( omega .* ( interp1( a, squeeze(afinal(t+1,:,:,:,:,dSp)), ap, 'spline' ) - ap) ) + exp( - omega .* ( interp1( a, squeeze(afinal(t+1,:,:,:,:,dSp)), ap, 'spline' ) - ap) ) + 2 ) + log(0.25) ), na, 1,1,1,nS ) ).^( theta.*(1-sigma)-1 ) .* theta .* ( psi .* ap + epsilon ).^( (1-theta).*(1-sigma) );

[Daniel kiracofe]

Hmmm. There are some suspicious lines here. Look at this one

ap1(:,:,:,:,:) = interp1(-squeeze(F1(:,t,aux1,aux2,aux3,aux4,aux5)),a,0,'spline');

Are you sure you didn't mean to do something like this instead?

ap1(aux1,aux2,aux3,aux4,aux5) = interp1(-squeeze(F1(:,t,aux1,aux2,aux3,aux4,aux5)),a,0,'spline');

As written, this implies that you are overwriting the previous value of ap1 every time through the loop. i.e. you calculate a value, through it away, repeat. You would only keep the last one. If this is really what you are doing, then pre-allocation of ap1 shouldn't matter, because it's not changing size during each execution of the loop. I suspect you are trying to post a simplified version of the code, and just made a cut and paste error. If it is your real code, then you can save a ton of time by eliminating all 5 loops and just doing

ap1(:,:,:,:,:) = interp1(-squeeze(F1(:,t,na,nh,nx,nA,nS)),a,0,'spline')

On the other hand, preallocation of variables like F1, F2, eta2, should save time, as those are growing in the loop.