2

ppval を使用して区分的多項式スプラインから一連の値を評価するループがあります。補間はループの中で最も時間のかかる部分であり、関数の効率を改善する方法を探しています。

より具体的には、有限差分スキームを使用して、摩擦溶接部の過渡温度分布を計算しています。これを行うには、時間ステップごとに材料特性を (温度と位置の関数として) 再計算する必要があります。レート制限要因は、これらの値の補間です。別の有限差分スキーム (時間領域での制限が少ない) を使用することもできますが、可能であれば、私が持っているものに固執したいと思います。

以下に MWE を含めました。

x=0:.1:10;
y=sin(x);
pp=spline(x,y);
tic
for n=1:10000
    x_int=10*rand(1000,1);
    y_int=ppval(pp,x_int);
end
toc

plot(x,y,x_int,y_int,'*') % plot for sanity of data

Elapsed time is 1.265442 seconds.

編集-値間の単純な線形補間に満足しているが、interp1関数はppvalよりも遅いことにおそらく言及する必要があります

x=0:.1:10;
y=sin(x);
tic
for n=1:10000
    x_int=10*rand(1000,1);
    y_int=interp1(x,y,x_int,'linear');
end
toc

plot(x,y,x_int,y_int,'*') % plot for sanity of data

Elapsed time is 1.957256 seconds.
4

2 に答える 2

2

JITの唯一の最も厄介な制限に直面しているため、これは遅いです。これは、SO の MATLAB タグに非常に多くの質問が寄せられる原因です。

MATLAB の JIT アクセラレータは、非組み込み関数を呼び出すループを高速化できません。

と は両方ppvalとinterp1も組み込まれていませんtype ppval(またはで確認してくださいedit interp1)。それらの実装は特に遅くはありませんが、ループに配置された場合は高速ではありません。

最近のバージョンの MATLAB では改善されている印象がありますが、「インライン」ループと「非インライン」ループの間にはまだ大きな違いがあります。彼らの JIT が非ビルトインに再帰するだけでこのタスクを自動化しない理由は、私にはまったくわかりません。

とにかく、これを修正するには、何が起こるかの本質をコピーppvalしてループ本体に貼り付ける必要があります。

% Example data
x = 0:.1:10;
y = sin(x);
pp = spline(x,y);


% Your original version
tic
for n = 1:10000
    x_int = 10*rand(1000,1);
    y_int = ppval(pp, x_int);
end
toc


% "inlined" version

tic

br = pp.breaks.';
cf = pp.coefs;

for n = 1:10000

    x_int = 10*rand(1000,1);

    [~, inds] = histc(x_int, [-inf; br(2:end-1); +inf]); 

    x_shf = x_int - br(inds);    
    zero  = ones(size(x_shf));
    one   = x_shf;
    two   = one .* x_shf;
    three = two .* x_shf;

    y_int = sum( [three two one zero] .* cf(inds,:), 2);
end
toc

プロファイラー:

ここに画像の説明を入力

私の安っぽいマシンでの結果：

Elapsed time is 2.764317 seconds.  % ppval
Elapsed time is 1.695324 seconds.  % "inlined" version

違いは実際には私が予想したよりも小さいですが、それは主にsum()-- このppval場合、通常、繰り返しごとに 1つのhistcサイトを評価するだけで済みます。 /ベクトル乗算x*y(BLAS) の代わりにsum(x.*y)(高速ですが、BLAS 高速ではありません)。

まあ、〜60％の削減は悪くありません:)

于 2013-09-06T09:27:51.183 に答える
2

interp1が よりも遅いのは少し意外ですがppval、そのソース コードをざっと見てみると、多くの特殊なケースをチェックする必要があり、すべてのポイントをループする必要があるようです。一定です。

I didn't check the timing, but I guess you can speed up the linear interpolation by a lot if you can guarantee that steps in x of your table are constant, and that the values to be interpolated are stricktly within the given range, so that you do not have to do any checking. In that case, linear interpolation can be converted to a simple lookup problem like so:

%data to be interpolated, on grid with constant step
x = 0:0.5:10;
y = sin(x);

x_int = 0:0.1:9.9;

%make sure it is interpolation, not extrapolation
assert(all(x(1) <= x_int & x_int < x(end)));

% compute mapping, this can be precomputed for constant grid
slope = (length(x) - 1) / (x(end) - x(1));
offset = 1 - slope*x(1); 

%map x_int to interval 1..lenght(i)
xmapped = offset + slope * x_int;
ind = floor(xmapped);
frac = xmapped - ind;
%interpolate by taking weighted sum of neighbouring points
y_int = y(ind) .* (1 - frac) + y(ind+1) .* frac;

% make plot to check correctness
plot(x, y, 'o-', x_int, y_int, '.')
于 2013-09-05T19:23:52.593 に答える