Julia 矩阵乘法比 numpy 的要慢答案

【问题标题】：Julia matrix multiplication is slower than numpy'sJulia 矩阵乘法比 numpy 的要慢
【发布时间】：2016-03-17 08:10:43
【问题描述】：

我正在尝试在 Julia 中进行一些矩阵乘法，以将其与 numpy 进行基准测试。

我的 Julia 代码如下：

function myFunc()
  A = randn(10000, 10000)
  B = randn(10000, 10000)
  return A*B
end

myFunc()

而python版本是：

A = np.random.rand(10000,10000)
B = np.random.rand(10000,10000)
A*B

Python 版本的执行时间不到 100 毫秒。 Julia 版本占用 13 秒！！看到他们在后台使用几乎相同的 BLAS 技术，Julia 版本似乎有什么问题？！

【问题讨论】：

标签： python numpy julia matrix-multiplication blas

【解决方案1】：

我不认为他们在做同样的事情。 numpy 表达式只执行逐个元素的乘法，而 Julia 表达式执行真正的矩阵乘法。

您可以通过使用较小的输入来看到差异。这是numpy 示例：

>>> A
array([1, 2, 3])
>>> B
array([[1],
       [2],
       [3]])
>>> A * B
array([[1, 2, 3],
       [2, 4, 6],
       [3, 6, 9]])
>>> B * A
array([[1, 2, 3],
       [2, 4, 6],
       [3, 6, 9]])

请注意，这里有 broadcasting，它“模拟”两个向量的外积，因此您可能会认为它是矩阵乘法。但它不可能，因为矩阵乘法不是可交换的，这里是(A * B) == (B * A)。看看当你在 Julia 中做同样的事情时会发生什么：

julia> A = [1, 2, 3]
3-element Array{Int64,1}:
 1
 2
 3

julia> B = [1 2 3]
1x3 Array{Int64,2}:
 1  2  3

julia> A * B
3x3 Array{Int64,2}:
 1  2  3
 2  4  6
 3  6  9

julia> B * A
1-element Array{Int64,1}:
 14

在这里，B * A 为您提供了适当的点积。如果您想进行真正的比较，请尝试numpy.dot。

如果您使用的是 Python 3.5 或更高版本，您还可以使用新的内置点积运算符！只需确保矩阵的形状对齐即可：

>>> A
array([[1, 2, 3]])
>>> B
array([[1],
       [2],
       [3]])
>>> A @ B
array([[14]])
>>> B @ A
array([[1, 2, 3],
       [2, 4, 6],
       [3, 6, 9]])

【讨论】：

感谢您启发我 Senderle。 numpy 中的以下代码比 Julia 版本花费的时间更长： A = np.mat(np.random.rand(10000,10000)); B = np.mat(np.random.rand(10000,10000)); C = A*B

【解决方案2】：

朴素矩阵乘法采用 N^3 次运算。您可以做一个简单的基准测试来查看这种增长：

function myFunc(N)
    A = rand(N, N)
    B = rand(N, N)

    A*B
end

myFunc(1)   # run once to compile

sizes = [floor(Int, x) for x in logspace(1, 3.5, 50)]

times = [@elapsed(myFunc(n)) for n in sizes]

using PyPlot

loglog(sizes, times, "o-")

为了更认真地做到这一点，我会在每种尺寸下平均运行几次。我得到类似下图的东西。事实上，外推到 N=10^4 在我的计算机上给出了大约 20 或 30 秒的时间。（同样，更严重的是，我会在对数图上拟合一条直线来进行外推。）

【讨论】：