一种对大量不同长度数组的 prod() 调用进行矢量化的方法？

Question

所以我的问题是，我想在没有 for 循环的情况下执行此操作。获取多个向量但长度不同的 prod()。

我正在处理与体素相交的光线。我通常有 1e6 个射线和 1e5 个体素，但这可能会有所不同。

intxRays 是与体素相交的光线列表。

gainList 是一个一维向量，每个元素都有一个值，对应于之前计算的特定射线体素交点（实际上是在你可爱的很多 here 的帮助下）。

rayIntxStart 和 rayIntxEnd 是索引的向量，在 gainlist 数组中，每条射线的对应值开始和结束（它们都是按顺序排列的）。

这是代码和一些示例以及预期的输出。

gainSum = zeros(1, 5);

% only interested in the intx uniques
intxSegCtr = 1;

% loop through all of the unique segments
for rayCtr = 1:max(intxRays)

    if rayCtr == intxRays(intxSegCtr)

        startInd = rayIntxStart(intxSegCtr);
        endInd = rayIntxEnd(intxSegCtr);

        % find which rows correspoond to those segements
        gainVals = gainList(startInd:endInd);
        gainProd = prod(gainVals);

        % get the product of the gains for those voxels
        gainSumIdx = intxRays(intxSegCtr);
        gainSum(gainSumIdx) = gainProd;

        % increment counter
        intxSegCtr = intxSegCtr + 1;

    end
end

五条射线和九个体素的示例数据。假设体素增益数组对于九个体素（在上一步中使用）看起来像这样（为简单起见）。

voxelGains = 10:10:90;

现在假设光线 1 和 3 不击中任何物体，光线 2 击中体素 1 和 2，光线 4 击中体素 2:7 并且射线 5 击中体素 6:9

intxRays = [2, 4, 5];

gainList = [10, 20, 20, 30, 40, 50, 60, 70, 60 70, 80, 90];

rayIntxStart = [1, 3, 9];

rayIntxEnd = [2, 8, 12];

对于这些数字，上面的代码会给出结果：

gainSum = [0, 200, 0, 5.0400e+09, 3.024e+07];

我希望这一切都有意义。

当我开发它时，我使用的光线和体素数要小得多，而且效果很好。不过，当我向上移动时，我的代码中的主要瓶颈是这个循环。实际上，仅 gainVals 和 gainProd 分配就占了我运行时间的 80% 和 15%。

这是我能找到的唯一可行的方法，填充等由于涉及的尺寸而不起作用。

有没有办法在没有这个循环的情况下获得我想要的值？

非常感谢！

Answer 1

好的，这是一个非常小的性能提升，但它可能会有所帮助。为了在没有循环的情况下测试矩阵方式，需要更大的数据样本。

这些是 3 个灵魂，您的原创，优化和优化的方式作为 oneliner。如果这已经为你做了什么，你能试试吗？

clear all
% orignial loop through all Rays
intxRays = [2, 4, 5];
gainList = [10, 20, 20, 30, 40, 50, 60, 70, 60 70, 80, 90];
rayIntxStart = [1, 3, 9];
rayIntxEnd = [2, 8, 12];
gainSum = zeros(1, 5);
tic
% only interested in the intx uniques
intxSegCtr = 1;

% loop through all of the unique segments
for rayCtr = 1:max(intxRays)

    if rayCtr == intxRays(intxSegCtr)

        startInd = rayIntxStart(intxSegCtr);
        endInd = rayIntxEnd(intxSegCtr);

        % find which rows correspoond to those segements
        gainVals = gainList(startInd:endInd);
        gainProd = prod(gainVals);

        % get the product of the gains for those voxels
        gainSumIdx = intxRays(intxSegCtr);
        gainSum(gainSumIdx) = gainProd;

        % increment counter
        intxSegCtr = intxSegCtr + 1;

    end
end
toc

clear all
%loop insted of every single one to max just through the intxRays
intxRays = [2, 4, 5];
gainList = [10, 20, 20, 30, 40, 50, 60, 70, 60 70, 80, 90];
rayIntxStart = [1, 3, 9];
rayIntxEnd = [2, 8, 12];
gainSum = zeros(1, 5);
tic
for rayCtr=1:length(intxRays)
    %no if as you just go through them
    %intxRays(rayCtr) is the corresponding element

     startInd = rayIntxStart(rayCtr);
     endInd = rayIntxEnd(rayCtr);
     % find which rows correspoond to those segements
     gainVals = gainList(startInd:endInd);
     gainProd = prod(gainVals);     

    % get the product of the gains for those voxels and set them to the ray
    gainSum(intxRays(rayCtr)) = gainProd;
end
%disp(gainSum);
toc

clear all
%same as above, but down to 1 line so no additional values are generated
intxRays = [2, 4, 5];
gainList = [10, 20, 20, 30, 40, 50, 60, 70, 60 70, 80, 90];
rayIntxStart = [1, 3, 9];
rayIntxEnd = [2, 8, 12];
gainSum = zeros(1, 5);
tic
for rayCtr=1:length(intxRays)
gainSum(intxRays(rayCtr))=prod(gainList(rayIntxStart(rayCtr):rayIntxEnd(rayCtr)));
end
toc