3

重心式による多項式補間を行うアルゴリズムを書くことになっている課題が与えられました。式は次のように述べています。

p（x）=（SIGMA_（j = 0 to n）w（j）* f（j）/（x --x（j）））/（SIGMA_（j = 0 to n）w（j）/（x --x（j）））

私はうまく機能するアルゴリズムを書きました、そして私は私が望む多項式出力を手に入れます。ただし、これには非常に長いループを使用する必要があり、グリッド数が多い場合は、多くの厄介なループ操作を実行する必要があります。したがって、これらすべてのループを回避するために、これをどのように改善できるかについて誰かが何かヒントを持っていれば、私は大いに感謝します。

アルゴリズムでは、補間することになっている与えられたポイントxをf表します。 wは、アルゴリズムを実行する前に計算された重心の重みを表します。そしてgrid、補間が行われるべきlinspaceは次のとおりです。

function p = barycentric_formula(x,f,w,grid)

%Assert x-vectors and f-vectors have same length.
if length(x) ~= length(f)
    sprintf('Not equal amounts of x- and y-values. Function is terminated.')
    return;
end

n = length(x);
m = length(grid);
p = zeros(1,m);

% Loops for finding polynomial values at grid points.  All values are
% calculated by the barycentric formula.
for i = 1:m
    var = 0;
    sum1 = 0;
    sum2 = 0;
    for j = 1:n
        if grid(i) == x(j)
            p(i) = f(j);
            var = 1;
        else
            sum1 = sum1 + (w(j)*f(j))/(grid(i) - x(j));
            sum2 = sum2 + (w(j)/(grid(i) - x(j)));
        end
    end
    if var == 0
        p(i) = sum1/sum2;
    end    
end
4

1 に答える 1

4

これは、matlabの「ベクトル化」の古典的なケースです。ループを削除するだけです。とても簡単です。まず、このコードを見てください：

function p = bf2(x, f, w, grid)

m = length(grid);
p = zeros(1,m);

for i = 1:m
    var = grid(i)==x;
    if any(var)
        p(i) = f(var);
    else
        sum1 = sum((w.*f)./(grid(i) - x));
        sum2 = sum(w./(grid(i) - x));
        p(i) = sum1/sum2;
    end
end
end

上の内側のループを削除しましたj。ここで行ったのは、実際には、インデックスを削除し(j)、算術演算子をtoから/toに./変更*することだけでし.*た。同じですが、操作が要素ごとに実行されることを示すために前にドットが付いています。これは、通常の行列演算子とは対照的に、配列演算子と呼ばれます。xまた、グリッドポイントが該当する特殊なケースの処理は、元の実装で行ったものと非常に似ていますが、のvarようなベクトルを使用するだけであることに注意してくださいx(var)==grid(i)。

これで、最も外側のループを削除することもできます。これはもう少し注意が必要です。MATLABでこれを行うには、2つの主要なアプローチがあります。私はそれをより簡単な方法で行います。これは効率が低下する可能性がありますが、より明確に読むことができます-使用repmat：

function p = bf3(x, f, w, grid)

% Find grid points that coincide with x.
% The below compares all grid values with all x values
% and returns a matrix of 0/1. 1 is in the (row,col)
% for which grid(row)==x(col)

var  = bsxfun(@eq, grid', x);

% find the logical indexes of those x entries
varx = sum(var, 1)~=0;

% and of those grid entries
varp = sum(var, 2)~=0;

% Outer-most loop removal - use repmat to
% replicate the vectors into matrices.
% Thus, instead of having a loop over j
% you have matrices of values that would be
% referenced in the loop

ww = repmat(w, numel(grid), 1);
ff = repmat(f, numel(grid), 1);
xx = repmat(x, numel(grid), 1);
gg = repmat(grid', 1, numel(x));

% perform the calculations element-wise on the matrices
sum1 = sum((ww.*ff)./(gg - xx),2);
sum2 = sum(ww./(gg - xx),2);
p    = sum1./sum2;

% fix the case where grid==x and return
p(varp) = f(varx);

end

bsxfun完全にベクトル化されたバージョンは、ではなくで実装できますrepmat。行列が明示的に形成されていないため、これは潜在的に少し速くなる可能性があります。ただし、システムサイズが小さい場合、速度差は大きくない場合があります。

また、ループが1つある最初のソリューションも、パフォーマンス的にはそれほど悪くありません。それらをテストして、何が良いかを確認することをお勧めします。たぶん、完全にベクトル化する価値はありませんか？最初のコードはもう少し読みやすく見えます。

于 2012-10-24T18:37:50.980 に答える