Speeding up (vectorizing) barycentric interpolation

2 次查看（过去 30 天）

Adam Pigon 2020-4-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/515399-speeding-up-vectorizing-barycentric-interpolation

评论： darova 2020-4-12

Dear All, I have the following problem:

is it possible to speed up the barycentric interpolation given below? I presume that vectorization could be particularly effective but I do not know how to do it. As the code has to be run many times (and the size of matrices involved is rather large) any increase in efficiency would be helpful. Could GPU computing, after vectorization, be of any help here?

na  = 100;
nh  = 100;
nm  = 100;
nap = na+2;
nhp = nh+2;
H   = randn(nap,nhp) ;
M   = randn(nap,nhp) ;
MP  = randn(nap,nhp) ;
ap_trial = NaN(nh,nm) ;
hp_trial = NaN(nh,nm) ;
mp_trial = NaN(nh,nm) ;
[h_grid,m_grid]   = ndgrid(1:nh,1:nm);
[ap_grid,hp_grid] = ndgrid(1:nap,1:nhp);
tic
for r0 = 0:1
    for r1 = 0:(nap-2)
        for r2 = 0:(nhp-2) % loops picking up triangles constructed by adjacent points in matrices H and M 
            h1 = H(1+r1+r0,1+r2+r0); 
            h2 = H(1+r1,2+r2); 
            h3 = H(2+r1,1+r2); 
            m1 = M(1+r1+r0,1+r2+r0); 
            m2 = M(1+r1,2+r2); 
            m3 = M(2+r1,1+r2);
            
            % barycentric interpolation (weights)
            w1 = ( (m2-m3)*( h_grid -h3) + (h3-h2)*( m_grid -m3) ) / ( (m2-m3)*(h1-h3) + (h3-h2)*(m1-m3) );
            w2 = ( (m3-m1)*( h_grid -h3) + (h1-h3)*( m_grid -m3) ) / ( (m2-m3)*(h1-h3) + (h3-h2)*(m1-m3) );
            w3 = 1 - w1 - w2;
                 
            % preventing extrapolation
            w1(w1 < 0) = NaN;
            w2(w2 < 0) = NaN;
            w3(w3 < 0) = NaN;
                  
	    % barycentric interpolation (results)
            ap = w1 * ap_grid(1+r1+r0,1+r2+r0) + w2 * ap_grid(1+r1,2+r2) + w3 * ap_grid(2+r1,1+r2); 
            hp = w1 * hp_grid(1+r1+r0,1+r2+r0) + w2 * hp_grid(1+r1,2+r2) + w3 * hp_grid(2+r1,1+r2); 
            mp = w1 * MP(1+r1+r0,1+r2+r0)      + w2 * MP(1+r1,2+r2)      + w3 * MP(2+r1,1+r2);
            
            ap_trial( isnan(ap) == 0 ) = ap( isnan(ap) == 0 );  
            hp_trial( isnan(hp) == 0 ) = hp( isnan(hp) == 0 );  
            mp_trial( isnan(mp) == 0 ) = mp( isnan(mp) == 0 ); 
        end
    end
end
toc

14 个评论
显示 12更早的评论隐藏 12更早的评论

Adam Pigon 2020-4-8

I see your point now. Unforunately, the example above using randn matrices may be a bit misleading as randn numbers are concentrated around 0. Therefore, we are operating with just a couple of indices like 1:5.

In reality:

Matrices ap_grid, hp_grid, h_grid and m_grid are exogenously given and they are regular, strictly monotonic, although not necessarily equidistant, grids.

Matrices H and M are endogenous - they are a function of many arguments like ap_grid and hp_grid matrices, some parameters, time iterations etc. and it is hard to say what their limits are. They are not regular grids. If they were, interpolation would be trivial (scatteredInterpolation is not an option though).

The goal of the entire code is to do the following:

Take a given triangle out of exogenous ap_grid and hp_grid matrices.
Take a corresponding triangle out of endogenous matrices H and M.
Ascertain which points of exogenous matrices h_grid and m_gris are inside of the corresponding H and M triangle.
Write down interpolated values of ap and hp for all points on the h_grid and m_grid (that's why we need overwriting of the ap_trial matrix).

We cannot say ex ante which points of the h_grid and m_grid are inside of a given H and M triangle. Maybe all of them are, maybe none. Practically, there are 1 or 2 or, most often, none.

I thought it could be an ample source of optimization. The question is how to effectively find only these parts of matrices that have non-NaN values and operate just with them? Would it be even possible to achieve a speed gain then?

In short, your idea is briliant but instead of indices 1:5 there should be endogenous indices x:y that change every loop iteration.

Adam Pigon 2020-4-12

Thank you for your idea! I have tested this improvement and my conclusions are as follows:

Matlab is pretty slow in indexing matrices (see this thread for interesting discussion: https://www.mathworks.com/matlabcentral/answers/54522-why-is-indexing-vectors-matrices-in-matlab-very-inefficient ), which means that my code needed some streamlining to achieve better performance. Including additional indices, i.e. of these points that are within the triangle, that have to be computed and implemented is costly. Therefore, dimensions nh and nm have to be large enough to achieve a good trade off beteween an additional cost of indexing and an efficiency gain from operating with smaller matrices that exploit only a fraction of matrices of size (nh x nm). In a nuthshell, if nm and nh are large enough, the efficiency gain can be sizable. An increase in nm is actually costless as nm does not show up in any of the loop sizes and therefore it doesn't put any additional indexing burden.

Thanks!

darova 2020-4-12

Thanks for the link

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

请先登录，再回答此问题。

类别

MATLAB Language Fundamentals Matrices and Arrays Matrix Indexing

在 Help Center 和 File Exchange 中查找有关 Matrix Indexing 的更多信息

产品

MATLAB

版本

R2018a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by

Speeding up (vectorizing) barycentric interpolation

14 个评论
显示 12更早的评论隐藏 12更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Speeding up (vectorizing) barycentric interpolation

14 个评论 显示 12更早的评论隐藏 12更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

14 个评论
显示 12更早的评论隐藏 12更早的评论