use*_*310 7 parallel-processing gradient linear-regression julia
不久前我被说服放弃了我舒适的matlab编程并开始在Julia编程.我一直在用神经网络工作很长时间,我认为,现在有了Julia,我可以通过并行计算梯度来更快地完成任务.
无需一次性对整个数据集计算梯度; 相反,人们可以拆分计算.例如,通过将数据集分成几部分,我们可以计算每个部分的部分梯度.然后通过将部分梯度相加来计算总梯度.
虽然原理很简单,但当我与Julia并行时,我会遇到性能下降,即一个进程比两个进程更快!我显然做错了什么......我已经咨询过论坛中提出的其他问题,但我仍然无法拼凑出答案.我认为我的问题在于有很多不必要的数据正在发生,但我无法正确修复它.
为了避免发布凌乱的神经网络代码,我发布了一个更简单的例子,它在线性回归的设置中复制了我的问题.
下面的代码块为线性回归问题创建了一些数据.代码解释了常量,但X是包含数据输入的矩阵.我们随机创建一个权重向量w ^当与乘以X创造一些目标ÿ.
######################################
## 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和X中恢复参数权重向量w:
####################################
## 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名工人和N = 10000,工作应划分如下:
不幸的是,如果我只有一个工作人员,那么整个代码的工作速 如果通过添加更多worker addprocs(),则代码运行得更慢.可以通过创建更多数据项来加剧此问题,例如使用N = 20000.随着更多的数据项,降级甚至更加明显.在我的特定计算环境中,N = 20000和一个核心,代码运行在~9秒.与N = 20000和4芯它需要〜18秒!
我尝试了很多不同的东西,这些东西都是在这个论坛的问题和答案的启发下,但遗憾的是无济于事.我意识到并行化是天真的,数据移动一定是问题,但我不知道如何正确地做到这一点.看来文档在这个问题上也有点稀缺(正如Ivo Balbaert的好书).
我很感激你的帮助,因为我已经坚持了很长一段时间,我真的需要它来完成我的工作.对于任何想要运行代码的人来说,为了省去复制粘贴的麻烦,你可以在这里获得代码.
感谢您抽出宝贵时间阅读这个非常冗长的问题!帮助我把它变成一个模型答案,任何Julia的新人都可以咨询!
我想说，GD 不是使用任何建议的方法进行并行化的良好候选者：或者SharedArray或DistributedArray，或者自己实现数据块的分布。
问题不在于 Julia，而在于 GD 算法。\n考虑代码：
\n\n主要流程：
\n\nfor iter = 1:iterations #iterations: "the more the better"\n \xce\xb4 = _gradient_descent_shared(X, y, \xce\xb8) \n \xce\xb8 = \xce\xb8 - \xce\xb1 * (\xce\xb4/N) \nend\n问题出在上面的 for 循环中，这是必须的。无论多好_gradient_descent_shared，迭代总数都会扼杀并行化的崇高概念。
阅读问题和上述建议后，我开始使用SharedArray. 请注意，我不是 SharedArrays 领域的专家。
主要流程部分（简单实现，无需正则化）：
\n\nrun_gradient_descent(X::SharedArray, y::SharedArray, \xce\xb8::SharedArray, \xce\xb1, iterations) = begin\n N = length(y)\n\n for iter = 1:iterations \n \xce\xb4 = _gradient_descent_shared(X, y, \xce\xb8) \n \xce\xb8 = \xce\xb8 - \xce\xb1 * (\xce\xb4/N) \n end \n \xce\xb8\nend\n\n_gradient_descent_shared(X::SharedArray, y::SharedArray, \xce\xb8::SharedArray, op=(+)) = begin \n if size(X,1) <= length(procs(X))\n return _gradient_descent_serial(X, y, \xce\xb8)\n else\n rrefs = map(p -> (@spawnat p _gradient_descent_serial(X, y, \xce\xb8)), procs(X))\n return mapreduce(r -> fetch(r), op, rrefs)\n end \nend\n所有工人通用的代码：
\n\n#= Returns the range of indices of a chunk for every worker on which it can work. \nThe function splits data examples (N rows into chunks), \nnot the parts of the particular example (features dimensionality remains intact).=# \n@everywhere function _worker_range(S::SharedArray)\n idx = indexpids(S)\n if idx == 0\n return 1:size(S,1), 1:size(S,2)\n end\n nchunks = length(procs(S))\n splits = [round(Int, s) for s in linspace(0,size(S,1),nchunks+1)]\n splits[idx]+1:splits[idx+1], 1:size(S,2)\nend\n\n#Computations on the chunk of the all data.\n@everywhere _gradient_descent_serial(X::SharedArray, y::SharedArray, \xce\xb8::SharedArray) = begin\n prange = _worker_range(X)\n\n pX = sdata(X[prange[1], prange[2]])\n py = sdata(y[prange[1],:])\n\n temp\xce\xb4 = pX\' * (pX * sdata(\xce\xb8) .- py) \nend\n数据加载和训练。让我假设我们有：
\n\n主要代码可能如下所示：
\n\nX=[ones(size(X,1)) X] #adding the artificial coordinate \nN, D = size(X)\nMAXITER = 500\n\xce\xb1 = 0.01\n\ninitial\xce\xb8 = SharedArray(Float64, (D,1))\nsX = convert(SharedArray, X)\nsy = convert(SharedArray, y) \nX = nothing\ny = nothing\ngc() \n\nfinal\xce\xb8 = run_gradient_descent(sX, sy, initial\xce\xb8, \xce\xb1, MAXITER);\n实施此操作并运行（在我的 Intel Clore i7 的 8 核上）后，我在训练多类（19 类）训练数据上比串行 GD（1 核）获得了非常轻微的加速（串行 GD 为 715 秒/对于串行 GD 为 665 秒）共享GD）。
\n\n如果我的实现是正确的（请检查一下 - 我指望它），那么 GD 算法的并行化就不值得了。当然，在 1 核上使用随机 GD 可能会获得更好的加速。
\n