Julia 中的并行梯度计算
Parallelising gradient calculation in Julia
前段时间我被说服放弃舒适的 matlab 编程并开始使用 Julia 编程。我长期以来一直在研究神经网络，我认为现在有了 Julia，我可以通过并行计算梯度来更快地完成工作。
不需要一次性在整个数据集上计算梯度；相反，可以拆分计算。例如，通过将数据集分成几部分，我们可以计算每个部分的部分梯度。然后通过将部分梯度相加来计算总梯度。
尽管原理很简单，但当我与 Julia 并行化时，性能下降了，即一个进程比两个进程快！我显然做错了什么......我已经咨询了论坛中提出的其他问题，但我仍然无法拼凑出答案。我认为我的问题在于有很多不必要的数据在移动，但我无法正确修复它。
为了避免发布乱七八糟的神经网络代码，我在下面发布了一个更简单的示例，它在线性回归的设置中复制了我的问题。
下面的代码块为线性回归问题创建了一些数据。代码解释了常量，但 X 是包含数据输入的矩阵。我们随机创建一个权重向量 w，当它与 X 相乘时创建一些目标 Y.
######################################
## CREATE LINEAR REGRESSION PROBLEM ##
######################################
# This code implements a simple linear regression problem
MAXITER = 100 # number of iterations for simple gradient descent
N = 10000 # number of data items
D = 50 # dimension of data items
X = randn(N, D) # create random matrix of data, data items appear row-wise
Wtrue = randn(D,1) # create arbitrary weight matrix to generate targets
Y = X*Wtrue # generate targets
下面的下一个代码块定义了用于衡量我们的回归（即负对数似然）和权重向量梯度的函数 w:
####################################
## DEFINE FUNCTIONS ##
####################################
@everywhere begin
#-------------------------------------------------------------------
function negative_loglikelihood(Y,X,W)
#-------------------------------------------------------------------
# number of data items
N = size(X,1)
# accumulate here log-likelihood
ll = 0
for nn=1:N
ll = ll - 0.5*sum((Y[nn,:] - X[nn,:]*W).^2)
end
return ll
end
#-------------------------------------------------------------------
function negative_loglikelihood_grad(Y,X,W, first_index,last_index)
#-------------------------------------------------------------------
# number of data items
N = size(X,1)
# accumulate here gradient contributions by each data item
grad = zeros(similar(W))
for nn=first_index:last_index
grad = grad + X[nn,:]' * (Y[nn,:] - X[nn,:]*W)
end
return grad
end
end
请注意，以上函数是有意未向量化的！我选择不矢量化，因为最终代码（神经网络案例）也不会接受任何矢量化（让我们不要深入了解这方面的更多细节）。
最后，下面的代码块显示了一个非常简单的梯度下降，它试图从给定数据 Y[=60 中恢复参数权重向量 w =] 和 X:
####################################
## SOLVE LINEAR REGRESSION ##
####################################
# start from random initial solution
W = randn(D,1)
# learning rate, set here to some arbitrary small constant
eta = 0.000001
# the following for-loop implements simple gradient descent
for iter=1:MAXITER
# get gradient
ref_array = Array(RemoteRef, nworkers())
# let each worker process part of matrix X
for index=1:length(workers())
# first index of subset of X that worker should work on
first_index = (index-1)*int(ceil(N/nworkers())) + 1
# last index of subset of X that worker should work on
last_index = min((index)*(int(ceil(N/nworkers()))), N)
ref_array[index] = @spawn negative_loglikelihood_grad(Y,X,W, first_index,last_index)
end
# gather the gradients calculated on parts of matrix X
grad = zeros(similar(W))
for index=1:length(workers())
grad = grad + fetch(ref_array[index])
end
# now that we have the gradient we can update parameters W
W = W + eta*grad;
# report progress, monitor optimisation
@printf("Iter %d neg_loglikel=%.4f\n",iter, negative_loglikelihood(Y,X,W))
end
正如所希望的那样，我在这里尝试以最简单的方式并行计算梯度。我的策略是在尽可能多的可用工作人员中中断梯度的计算。每个工人只需要在矩阵 X 的一部分上工作，该部分由 first_index 和 last_index 指定。因此，每个工人应该与 X[first_index:last_index,:] 一起工作。比如4个worker，N=10000，分工如下：
- 工人 1 => first_index = 1，last_index = 2500
- 工人 2 => first_index = 2501，last_index = 5000
- 工人 3 => first_index = 5001，last_index = 7500
- 工人 4 => first_index = 7501，last_index = 10000
不幸的是，如果我只有一个工人，整个代码运行得更快。如果通过 addprocs() 添加更多工人，代码 运行 会变慢。可以通过创建更多数据项来加剧这个问题，例如使用 N=20000。
随着更多的数据项，退化甚至更加明显。
在我具有 N=20000 和一个核心的特定计算环境中，代码 运行s 在 ~9 秒内完成。使用 N=20000 和 4 个内核需要大约 18 秒！
受此论坛中问题和答案的启发，我尝试了许多不同的方法，但遗憾的是都无济于事。我意识到并行化是天真的，数据移动一定是问题所在，但我不知道如何正确地做到这一点。似乎关于这个问题的文档也有点稀缺（Ivo Balbaert 的好书也是如此）。
非常感谢你的帮助，因为我已经被这个问题困扰了很长一段时间，我的工作确实需要它。对于任何想要 运行 代码的人，为了省去复制粘贴的麻烦，您可以获得代码 here.
感谢您花时间阅读这个冗长的问题！帮我把它变成一个模范答案，任何刚接触 Julia 的人都可以参考！
如果您想减少数据移动量，您应该强烈考虑使用 SharedArrays。您可以只预分配一个输出向量，并将其作为参数传递给每个工作人员。正如您所建议的，每个工人都设置了其中的一块。
我会说 GD 不适合使用任何建议的方法对其进行并行化：SharedArray 或 DistributedArray，或者自己实现数据块分布。
问题不在Julia，而在GD算法。
考虑代码：
主要流程：
for iter = 1:iterations #iterations: "the more the better"
δ = _gradient_descent_shared(X, y, θ)
θ = θ - α * (δ/N)
end
问题出在上面的for循环中，这是必须的。 _gradient_descent_shared再好，总的迭代次数扼杀了并行化的崇高概念。
阅读问题和上述建议后，我开始使用 SharedArray 实现 GD。请注意，我不是 SharedArrays 领域的专家。
主要流程部分（无正则化简单实现）：
run_gradient_descent(X::SharedArray, y::SharedArray, θ::SharedArray, α, iterations) = begin
N = length(y)
for iter = 1:iterations
δ = _gradient_descent_shared(X, y, θ)
θ = θ - α * (δ/N)
end
θ
end
_gradient_descent_shared(X::SharedArray, y::SharedArray, θ::SharedArray, op=(+)) = begin
if size(X,1) <= length(procs(X))
return _gradient_descent_serial(X, y, θ)
else
rrefs = map(p -> (@spawnat p _gradient_descent_serial(X, y, θ)), procs(X))
return mapreduce(r -> fetch(r), op, rrefs)
end
end
所有工人通用的代码：
#= Returns the range of indices of a chunk for every worker on which it can work.
The function splits data examples (N rows into chunks),
not the parts of the particular example (features dimensionality remains intact).=#
@everywhere function _worker_range(S::SharedArray)
idx = indexpids(S)
if idx == 0
return 1:size(S,1), 1:size(S,2)
end
nchunks = length(procs(S))
splits = [round(Int, s) for s in linspace(0,size(S,1),nchunks+1)]
splits[idx]+1:splits[idx+1], 1:size(S,2)
end
#Computations on the chunk of the all data.
@everywhere _gradient_descent_serial(X::SharedArray, y::SharedArray, θ::SharedArray) = begin
prange = _worker_range(X)
pX = sdata(X[prange[1], prange[2]])
py = sdata(y[prange[1],:])
tempδ = pX' * (pX * sdata(θ) .- py)
end
数据加载和训练。让我假设我们有：
- X::Array 大小为 (N,D) 的特征，其中 N - 示例数，特征的 D 维
- y::Array 大小的标签 (N,1)
主要代码可能如下所示：
X=[ones(size(X,1)) X] #adding the artificial coordinate
N, D = size(X)
MAXITER = 500
α = 0.01
initialθ = SharedArray(Float64, (D,1))
sX = convert(SharedArray, X)
sy = convert(SharedArray, y)
X = nothing
y = nothing
gc()
finalθ = run_gradient_descent(sX, sy, initialθ, α, MAXITER);
在实施这个和 运行（在我的 Intell Clore i7 的 8 核上）之后，我在训练多类（19 类 ) 训练数据（串行 GD 为 715 秒/共享 GD 为 665 秒）。
如果我的实现是正确的（请检查一下 - 我指望它）那么 GD 算法的并行化是不值得的。绝对可以在 1 核上使用随机 GD 获得更好的加速。
前段时间我被说服放弃舒适的 matlab 编程并开始使用 Julia 编程。我长期以来一直在研究神经网络，我认为现在有了 Julia，我可以通过并行计算梯度来更快地完成工作。
不需要一次性在整个数据集上计算梯度；相反，可以拆分计算。例如，通过将数据集分成几部分，我们可以计算每个部分的部分梯度。然后通过将部分梯度相加来计算总梯度。
尽管原理很简单，但当我与 Julia 并行化时，性能下降了，即一个进程比两个进程快！我显然做错了什么......我已经咨询了论坛中提出的其他问题，但我仍然无法拼凑出答案。我认为我的问题在于有很多不必要的数据在移动，但我无法正确修复它。
为了避免发布乱七八糟的神经网络代码，我在下面发布了一个更简单的示例，它在线性回归的设置中复制了我的问题。
下面的代码块为线性回归问题创建了一些数据。代码解释了常量，但 X 是包含数据输入的矩阵。我们随机创建一个权重向量 w，当它与 X 相乘时创建一些目标 Y.
######################################
## CREATE LINEAR REGRESSION PROBLEM ##
######################################
# This code implements a simple linear regression problem
MAXITER = 100 # number of iterations for simple gradient descent
N = 10000 # number of data items
D = 50 # dimension of data items
X = randn(N, D) # create random matrix of data, data items appear row-wise
Wtrue = randn(D,1) # create arbitrary weight matrix to generate targets
Y = X*Wtrue # generate targets
下面的下一个代码块定义了用于衡量我们的回归（即负对数似然）和权重向量梯度的函数 w:
####################################
## DEFINE FUNCTIONS ##
####################################
@everywhere begin
#-------------------------------------------------------------------
function negative_loglikelihood(Y,X,W)
#-------------------------------------------------------------------
# number of data items
N = size(X,1)
# accumulate here log-likelihood
ll = 0
for nn=1:N
ll = ll - 0.5*sum((Y[nn,:] - X[nn,:]*W).^2)
end
return ll
end
#-------------------------------------------------------------------
function negative_loglikelihood_grad(Y,X,W, first_index,last_index)
#-------------------------------------------------------------------
# number of data items
N = size(X,1)
# accumulate here gradient contributions by each data item
grad = zeros(similar(W))
for nn=first_index:last_index
grad = grad + X[nn,:]' * (Y[nn,:] - X[nn,:]*W)
end
return grad
end
end
请注意，以上函数是有意未向量化的！我选择不矢量化，因为最终代码（神经网络案例）也不会接受任何矢量化（让我们不要深入了解这方面的更多细节）。
最后，下面的代码块显示了一个非常简单的梯度下降，它试图从给定数据 Y[=60 中恢复参数权重向量 w =] 和 X:
####################################
## SOLVE LINEAR REGRESSION ##
####################################
# start from random initial solution
W = randn(D,1)
# learning rate, set here to some arbitrary small constant
eta = 0.000001
# the following for-loop implements simple gradient descent
for iter=1:MAXITER
# get gradient
ref_array = Array(RemoteRef, nworkers())
# let each worker process part of matrix X
for index=1:length(workers())
# first index of subset of X that worker should work on
first_index = (index-1)*int(ceil(N/nworkers())) + 1
# last index of subset of X that worker should work on
last_index = min((index)*(int(ceil(N/nworkers()))), N)
ref_array[index] = @spawn negative_loglikelihood_grad(Y,X,W, first_index,last_index)
end
# gather the gradients calculated on parts of matrix X
grad = zeros(similar(W))
for index=1:length(workers())
grad = grad + fetch(ref_array[index])
end
# now that we have the gradient we can update parameters W
W = W + eta*grad;
# report progress, monitor optimisation
@printf("Iter %d neg_loglikel=%.4f\n",iter, negative_loglikelihood(Y,X,W))
end
正如所希望的那样，我在这里尝试以最简单的方式并行计算梯度。我的策略是在尽可能多的可用工作人员中中断梯度的计算。每个工人只需要在矩阵 X 的一部分上工作，该部分由 first_index 和 last_index 指定。因此，每个工人应该与 X[first_index:last_index,:] 一起工作。比如4个worker，N=10000，分工如下：
- 工人 1 => first_index = 1，last_index = 2500
- 工人 2 => first_index = 2501，last_index = 5000
- 工人 3 => first_index = 5001，last_index = 7500
- 工人 4 => first_index = 7501，last_index = 10000
不幸的是，如果我只有一个工人，整个代码运行得更快。如果通过 addprocs() 添加更多工人，代码 运行 会变慢。可以通过创建更多数据项来加剧这个问题，例如使用 N=20000。
随着更多的数据项，退化甚至更加明显。
在我具有 N=20000 和一个核心的特定计算环境中，代码 运行s 在 ~9 秒内完成。使用 N=20000 和 4 个内核需要大约 18 秒！
受此论坛中问题和答案的启发，我尝试了许多不同的方法，但遗憾的是都无济于事。我意识到并行化是天真的，数据移动一定是问题所在，但我不知道如何正确地做到这一点。似乎关于这个问题的文档也有点稀缺（Ivo Balbaert 的好书也是如此）。
非常感谢你的帮助，因为我已经被这个问题困扰了很长一段时间，我的工作确实需要它。对于任何想要 运行 代码的人，为了省去复制粘贴的麻烦，您可以获得代码 here.
感谢您花时间阅读这个冗长的问题！帮我把它变成一个模范答案，任何刚接触 Julia 的人都可以参考！
如果您想减少数据移动量，您应该强烈考虑使用 SharedArrays。您可以只预分配一个输出向量，并将其作为参数传递给每个工作人员。正如您所建议的，每个工人都设置了其中的一块。
我会说 GD 不适合使用任何建议的方法对其进行并行化：SharedArray 或 DistributedArray，或者自己实现数据块分布。
问题不在Julia，而在GD算法。 考虑代码：
主要流程：
for iter = 1:iterations #iterations: "the more the better"
δ = _gradient_descent_shared(X, y, θ)
θ = θ - α * (δ/N)
end
问题出在上面的for循环中，这是必须的。 _gradient_descent_shared再好，总的迭代次数扼杀了并行化的崇高概念。
阅读问题和上述建议后，我开始使用 SharedArray 实现 GD。请注意，我不是 SharedArrays 领域的专家。
主要流程部分（无正则化简单实现）：
run_gradient_descent(X::SharedArray, y::SharedArray, θ::SharedArray, α, iterations) = begin
N = length(y)
for iter = 1:iterations
δ = _gradient_descent_shared(X, y, θ)
θ = θ - α * (δ/N)
end
θ
end
_gradient_descent_shared(X::SharedArray, y::SharedArray, θ::SharedArray, op=(+)) = begin
if size(X,1) <= length(procs(X))
return _gradient_descent_serial(X, y, θ)
else
rrefs = map(p -> (@spawnat p _gradient_descent_serial(X, y, θ)), procs(X))
return mapreduce(r -> fetch(r), op, rrefs)
end
end
所有工人通用的代码：
#= Returns the range of indices of a chunk for every worker on which it can work.
The function splits data examples (N rows into chunks),
not the parts of the particular example (features dimensionality remains intact).=#
@everywhere function _worker_range(S::SharedArray)
idx = indexpids(S)
if idx == 0
return 1:size(S,1), 1:size(S,2)
end
nchunks = length(procs(S))
splits = [round(Int, s) for s in linspace(0,size(S,1),nchunks+1)]
splits[idx]+1:splits[idx+1], 1:size(S,2)
end
#Computations on the chunk of the all data.
@everywhere _gradient_descent_serial(X::SharedArray, y::SharedArray, θ::SharedArray) = begin
prange = _worker_range(X)
pX = sdata(X[prange[1], prange[2]])
py = sdata(y[prange[1],:])
tempδ = pX' * (pX * sdata(θ) .- py)
end
数据加载和训练。让我假设我们有：
- X::Array 大小为 (N,D) 的特征，其中 N - 示例数，特征的 D 维
- y::Array 大小的标签 (N,1)
主要代码可能如下所示：
X=[ones(size(X,1)) X] #adding the artificial coordinate
N, D = size(X)
MAXITER = 500
α = 0.01
initialθ = SharedArray(Float64, (D,1))
sX = convert(SharedArray, X)
sy = convert(SharedArray, y)
X = nothing
y = nothing
gc()
finalθ = run_gradient_descent(sX, sy, initialθ, α, MAXITER);
在实施这个和 运行（在我的 Intell Clore i7 的 8 核上）之后，我在训练多类（19 类 ) 训练数据（串行 GD 为 715 秒/共享 GD 为 665 秒）。
如果我的实现是正确的（请检查一下 - 我指望它）那么 GD 算法的并行化是不值得的。绝对可以在 1 核上使用随机 GD 获得更好的加速。