Julia实现mul！当Y=A时

Question

矩阵乘法Y=A*B可以用mul！(Y，A，B(以节省内存分配。但是mul！如果Y=A，则不能使用。是否有类似的有效方法来计算Y*=B？或者，如果不是，做矩阵乘法的最有效方法是什么Y*=B

小型工作示例：

n = 10
A = rand(n,n)
B = rand(n,n)
Y = zeros(n,n)
#mul! removes allocations 
@allocated Y = A * B   #896
@allocated mul!(Y, A, B) #0    
#mul! can't be applied in this case
@allocated Y *= B #896
#desired function performance
@allocated mul_2!(Y, B)  #0

提前感谢您的帮助！

Answer 1

我认为你不能有效地实现这一点(由于矩阵乘法中的阻塞是如何工作的(。您最好只保留另一个适当大小的矩阵作为缓冲区。

Answer 2

尽管没有任何分配很难做到，正如其他答案所解释的那样。通过将矩阵乘法分解为一系列独立的向量矩阵乘法，可以减少对O(n(而不是O(n^2(的分配。所以mul_2!看起来像这样：

function mul_2!(Y, B)
tmprow = Matrix{eltype(Y)}(undef,1,size(Y,2))
foreach(eachrow(Y)) do r
mul!(tmprow, r', B) 
r .= tmprow'
end
end

以及OP中调用的输出：

julia> @allocated mul_2!(Y, B)  # 144 instead of 0 ¯_(ツ)_/¯
144

为了将分配一直降到零，可以将临时向量作为参数传递给mul_2!并进行预分配。这意味着如果多次重复操作，则在计算过程中不会发生分配。

大致如下：

tmprow = Matrix{eltype(Y)}(undef,1,size(Y,2))
for i in 1:lots_of_iterations
mul_2_notmp!(Y,B,tmprow)
end

就速度而言，分配是缓慢的，所以没有分配更好。但确切的权衡应该取决于矩阵大小，因为标准矩阵乘法库已经经过了大量的优化。