2

x入力として 1xn ベクトルを取り、nxn マトリックスを返す関数に取り組んでいますL。
ループをベクトル化して速度を上げたいのですが、私を困惑させる問題があります。ループ インデックスはループ インデックスbに依存しますa。どんな助けでも大歓迎です。

x = x(:); 
n = length(x);
L = zeros(n, n);
for a = 1 : n,
    for b = 1 : a-1,
        c = b+1 : a-1;
        if all(x(c)' < x(b) + (x(a) - x(b)) * ((b - c)/(b-a))),
            L(a,b) = 1;
        end
    end
end
4

3 に答える 3

1

簡単なテストから、下の三角形のみで何かをしているように見えます。これind2subのような醜いトリックを使用してベクトル化できるかもしれませんarrayfun

tril_lin_idx = find(tril(ones(n), -1));
[A, B] = ind2sub([n,n], tril_lin_idx);
C = arrayfun(@(a,b) b+1 : a-1, A, B, 'uniformoutput', false); %cell array
f = @(a,b,c) all(x(c{:})' < x(b) + (x(a) - x(b)) * ((b - c{:})/(b-a)));
L = zeros(n, n);
L(tril_lin_idx) = arrayfun(f, A, B, C);

x私はそれを持っておらず、期待される結果がわからないので、テストできません。私は通常、ベクトル化されたソリューションが好きですが、これはおそらく少しやりすぎです:)。私はあなたの明示的な for ループに固執します。これはより明確であり、Matlab のどの JIT が簡単に高速化できるはずです。if を に置き換えることができL(a,b) = all(...)ます。

編集1

~ n^3のスペースを無駄にしないように更新されたバージョンC:

tril_lin_idx = find(tril(ones(n), -1));
[A, B] = ind2sub([n,n], tril_lin_idx);
c = @(a,b) b+1 : a-1;
f = @(a,b) all(x(c(a, b))' < x(b) + (x(a) - x(b)) * ((b - c(a, b))/(b-a)));
L = zeros(n, n);
L(tril_lin_idx) = arrayfun(f, A, B);

編集2

ind2sub を使用せずb、より複雑な方法でa. 速度を上げるためにインライン化cしましたが、特に関数ハンドルの呼び出しはコストがかかるようです。

[A,B] = ndgrid(1:n);
v = B<A; % which elements to evaluate
f = @(a,b) all(x(b+1:a-1)' < x(b) + (x(a) - x(b)) * ((b - (b+1:a-1))/(b-a)));
L = false(n);
L(v) = arrayfun(f, A(v), B(v));
于 2013-08-27T21:45:47.430 に答える
1

あなたの問題を正しく理解してL(a, b) == 1いれば、a < c < b の任意の c について、(c, x(c)) が (a, x(a)) と (b, x(b)) を結ぶ線の「下」にある場合、 右？

これはベクトル化ではありませんが、別のアプローチを見つけました。新しい b ごとにすべての c を a < c < b と比較するのではなく、a から c への最大勾配を (a, b) に保存し、それを (a, b + 1) に使用しました。（一方向のみで試しましたが、両方向でも可能だと思います。）

x = x(:);
n = length(x);
L = zeros(n);

for a = 1:(n - 1)
  L(a, a + 1) = 1;
  maxSlope = x(a + 1) - x(a);
  for b = (a + 2):n
    currSlope = (x(b) - x(a)) / (b - a);
    if currSlope > maxSlope
      maxSlope = currSlope;
      L(a, b) = 1;
    end
  end
end

あなたのデータはわかりませんが、ランダムなデータを使用すると、結果は元のコードと同じになります（転置あり）。

于 2013-08-28T10:23:07.743 に答える