是否有比这更好的矢量化技术？

Question

我试图看看是否有其他方法可以更有效地编码此代码示例.这里,y是1xM矩阵(例如,1x1000),z是NxM矩阵(例如,5x1000).

mean(ones(N,1)*y.^3 .* z,2)

这段代码工作正常,但我担心如果N增加很多,ones(N,1)*y.^3可能会浪费太多,让一切都变慢.

思考？

Answer 1

对于一个小的矩阵来说,这并不可怕.很多时候你可以在这样的问题上使用bsxfun.这里的矩阵实在太小了,无法获得任何东西.

>> N = 5;M =1000;
>> y = rand(1,M);
>> z = rand(N,M);
>> mean(ones(N,1)*y.^3 .* z,2)
ans =
      0.12412
      0.11669
      0.12102
      0.11976
      0.12196

>> mean(bsxfun(@times,y.^3,z),2)
ans =
      0.12412
      0.11669
      0.12102
      0.11976
      0.12196

>> z*y.'.^3/M
ans =
      0.12412
      0.11669
      0.12102
      0.11976
      0.12196

如您所见,所有三种解决方案都返回相同的结果.一切都同样有效.

现在我将比较所需的时间.

>> timeit(@() mean(ones(N,1)*y.^3 .* z,2))
ans =
   0.00023018

>> timeit(@() mean(bsxfun(@times,y.^3,z),2))
ans =
   0.00026829

>> timeit(@() z*y.'.^3/M)
ans =
   0.00016594

正如我所说,你没有获得太多收益.事实上,bsxfun根本没有获得,甚至有点慢.但是如果你把表达重写成我提出的第三种形式,你可以获得一点.不多,但有点.

编辑:如果N很大,那么时间会有所改变.

>> N = 2000;M = 1000;
>> y = rand(1,M);
>> z = rand(N,M);
>> timeit(@() mean(ones(N,1)*y.^3 .* z,2))
ans =
     0.034664

>> timeit(@() mean(bsxfun(@times,y.^3,z),2))
ans =
     0.012234

>> timeit(@() z*y.'.^3/M)
ans =
    0.0017674

不同之处在于第一个解决方案显式创建了扩展的y.^ 3矩阵.这是低效的.

bsxfun解决方案更好,因为它从未明确地形成扩展y.^ 3矩阵.但它仍然形成一个N乘M的产品矩阵.所以这个解决方案仍然必须抓住并填充大量内存.

你应该理解为什么矩阵向量乘法在所有情况下都是最好的.没有形成大的矩阵.由于*只是一个点积(因此是产品的总和),因此它必须更有效.然后我在事后划分M以创建所需的均值.

比上一次略有改进......

>> timeit(@() z*(y.*y.*y).'/M)
ans =
    0.0015793

这比电力运行略有增加.

和timeit？这来自文件交换,这是一个非常有用的实用程序,由Steve Eddins编写的时间代码片段.