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

请注意，上述函数是故意未矢量化的！我选择不进行矢量化，因为最终代码（神经网络情况）也不会接受任何矢量化（让我们不要深入了解这方面的更多细节）。

最后，下面的代码块显示了一个非常简单的梯度下降，试图恢复参数权重向量w从给定的数据Y and 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 的一部分进行工作，该部分由第一个索引 and 最后索引。因此，每个工人都应该与X[first_index:last_index,:]。例如，对于 4 个工人，N = 10000，工作应划分如下：

• 工人 1 => 第一个索引 = 1，最后一个索引 = 2500
• 工人 2 => 第一个索引 = 2501，最后一个索引 = 5000
• 工人 3 => 第一个索引 = 5001，最后一个索引 = 7500
• 工人 4 => 第一个索引 = 7501，最后一个索引 = 10000

不幸的是，如果我只有一名工作人员，则整个代码的运行速度会更快。如果通过添加更多工人addprocs()，代码运行速度较慢。人们可以通过创建更多数据项来加剧这一问题，例如使用N=20000。 随着数据项的增加，性能下降更加明显。 在我的特定计算环境中N=20000和一个核心，代码运行时间约为 9 秒。和N=200004 核则需要约 18 秒！

受到这个论坛中的问题和答案的启发，我尝试了很多不同的方法，但不幸的是没有成功。我意识到并行化很幼稚，数据移动一定是问题所在，但我不知道如何正确地做到这一点。似乎关于这个问题的文档也有点稀缺（Ivo Balbaert 的好书也是如此）。

我非常感谢您的帮助，因为我已经被困在这个问题上很长一段时间了，我的工作真的需要它。对于任何想要运行代码的人，为了省去复制粘贴的麻烦，您可以获得代码here https://www.dropbox.com/s/e22apyk5t3ucymz/simpleCode_1st_parallel_attempt.jl?dl=0.

感谢您花时间阅读这个很长的问题！帮助我将其变成一个模型答案，任何刚接触 Julia 的人都可以参考！

我想说，GD 并不是使用任何建议的方法进行并行化的良好候选者：要么SharedArray or DistributedArray，或者自己实现数据块的分布。

问题不在于Julia，而在于GD算法。 考虑代码：

主要流程：

for iter = 1:iterations #iterations: "the more the better"
    δ = _gradient_descent_shared(X, y, θ) 
    θ = θ -  α * (δ/N)   
end

问题出在上面的for循环中，这是必须的。不管多好_gradient_descent_shared也就是说，迭代总数扼杀了并行化的崇高概念。

阅读问题和上述建议后，我开始使用SharedArray。请注意，我不是 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);

实施此操作并运行（在我的 Intel Clore i7 的 8 核上）后，我在训练多类（19 类）训练数据上比串行 GD（1 核）获得了非常轻微的加速（串行 GD 为 715 秒/665 秒为共享GD）。

如果我的实现是正确的（请检查一下 - 我指望它），那么 GD 算法的并行化就不值得了。当然，在 1 核上使用随机 GD 可能会获得更好的加速。