Julia 矩阵乘法比 numpy 的慢
Julia matrix multiplication is slower than numpy's
我正在尝试在 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 版本似乎有什么问题?!
我不认为他们在做同样的事情。 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]])
请注意,这里我们有 广播,它 "simulates" 两个向量的外积,因此您可能认为它是矩阵乘法。但它不可能,因为矩阵乘法是不可交换的,这里 (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]])
朴素矩阵乘法采用 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 秒。 (再一次,更严重的是,我会用一条直线拟合对数对数图来进行外推。)
我正在尝试在 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 版本似乎有什么问题?!
我不认为他们在做同样的事情。 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]])
请注意,这里我们有 广播,它 "simulates" 两个向量的外积,因此您可能认为它是矩阵乘法。但它不可能,因为矩阵乘法是不可交换的,这里 (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]])
朴素矩阵乘法采用 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-")
为了更认真地做到这一点,我会对每个尺寸的几次运行进行平均。
我得到类似下图的结果。