2

我正在研究一个以 1xn 向量x作为输入并返回 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

一个深奥的答案：您可以从 1:n 开始计算每个 a、b、c，排除无关因素，然后沿着 c 维度进行所有计算。

[a, b, c] = ndgrid(1:n, 1:n, 1:n);

La = x(c)' < x(b) + (x(a) - x(b)) .* ((b - c)./(b-a));
La(b >= a | c <= b | c >= a) = true;

L = all(La, 3);

尽管 jit 可能会很好地处理 for 循环，因为它们做的很少。

编辑：仍然使用所有内存，但数学较少

[A, B, C] = ndgrid(1:n, 1:n, 1:n);

valid = B < A & C > B & C < A;
a = A(valid); b = B(valid); c = C(valid);

La = true(size(A));
La(valid) = x(c)' < x(b) + (x(a) - x(b)) .* ((b - c)./(b-a));
L = all(La, 3);

Edit2：替换最后一行以添加没有元素的 c 为 true 的子句

L = all(La,3) | ~any(valid,3);
于 2013-08-29T03:42:29.883 回答
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, b) 中从 a 到 c 的最大斜率，并将其用于 (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 回答