向量化一个短循环，使 Octave 匹配 Matlab 速度

Posted 2023-02-16

【中文标题】向量化一个短循环，使 Octave 匹配 Matlab 速度

【英文标题】：Vectorizing a short loop such that Octave matches Matlab speed

【发布时间】：2012-03-14 00:09:13

【问题描述】：

我的代码在 Octave 中运行大约需要 0.003 秒，在 Matlab 中运行大约需要 0.0007 秒。由于 Octave 没有 JIT 编译，我想 Matlab 正在做我应该自己做的幕后优化。

numerators = zeros(1, 64);
for c = 1 : C
  numerators(c) = py(c) * prod(diag(pxc(:, x)));
end

py 是一个1xC 行向量。 px 是一个包含C 元素的数组，每个元素都是一个DxV 矩阵。 x 是一个 Dx1 列向量，其值在 [1-V] 上是离散的。

prod(diag(...))oddity 是 pxc(d, x(d)) 与所有 d 相乘的矢量化方式：

p = 1;
for d = 1 : D
  p *= pxc(d, x(d))
endfor

cellfun 可能有效，但我陷入了传递参数的细​​节中。 （如果它可以完成，就这么说，我会自己弄清楚如何）。另一种选择可能是对px 使用 3-D 矩阵，但是，我尝试了这个并且已经足够新手了，我无法让任何工作。

【问题讨论】：

您能否添加一些有关此过程目标的信息，以便我们更好地了解如何重新制定它？

时间太短是不准确的。在这种情况下最好重复代码 1000 次。

tmpearce -- 如果您了解朴素贝叶斯，我正在尝试计算所有类的 p(y=k|x) 的分子，因此，p(y=k)*PROD_d p( x_d|y=k) 对于所有 k。因此，py 是 y 上的多项分布（具有 64 个类），而 px 是 x 维度上的类条件 (pxc) 分布，因此，x 具有 D 维度和 V 个可能值，pxc是一个 DxV 矩阵。

Richie -- 上面的代码对每个数据点运行 100,000 次。八度确实慢了约 2-3 倍。

【参考方案1】：

我会根据您提供的信息试一试。不过我没有时间信息。

使 px 成为 3D 矩阵 DxVxC（而不是 1xC 单元阵列 DxV 2D 矩阵） px_new=px(:,x,:) 重新组织矩阵，将您感兴趣的值放入前两个维度的主对角线 在对角线上创建一个逻辑索引 (mask=eye(D,D);)； repmat 所以它是 C 在第三维中的大小。 索引到px_new 和reshape，这样就有C 列。 prod 列，留下一个1xC 向量

将此（按元素）与py 相乘以获得您的输出

px = nan(3,4,5);    %# create test 3D matrix  
px(:, :, 1) = [1 2 3 4; 4 5 6 4; 7 8 9 4]; 
px(:,:,2)=px(:,:,1)*1.5;  
px(:,:,3)=px(:,:,2)*1.5;  
px(:,:,4)=px(:,:,3)*1.5;      
px(:,:,5)=px(:,:,4)*1.5;  
x = [4 2 3];            %# 1xD vector discrete on 1-V  
px_new=px(:,x,:);       %# reorganize into DxDxC  

idx=logical(repmat(eye(size(pd_new,1))),[1,1,size(pd_new,3)])); %# logical index  
P = prod(reshape(pd_new(idx),[],size(pd_new,3))); %#P is now 1xC vector

此代码已在 http://www.online-utility.org/math/math_calculator.jsp 测试

编辑：我最初采取了一些不必要的步骤。我已将其更新为更简洁。