mul 的 Julia 实现！当 Y = A？答案

【问题标题】：Julia implementation of mul! when Y = A?mul 的 Julia 实现！当 Y = A？
【发布时间】：2022-10-23 18:18:48
【问题描述】：

矩阵乘法 Y = A * B 可以通过 mul!(Y, A, B) 来实现，以节省内存分配。但是穆尔！如果 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

在此先感谢您的帮助！

【问题讨论】：

标签： matrix julia matrix-multiplication memory-efficient

【解决方案1】：

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

【讨论】：

谢谢，很有帮助！作为对未来用户的说明，如果您有一个 for 循环并且想要实现它，您可以通过 for i = 1:(i_max) if i % 2 = 0 mul!(Y1, Y0, B_i) else mul! (Y0, Y1, B_i) end end 其中 B_i 是需要为相关迭代添加的任何矩阵

【解决方案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

【讨论】：