1、矢量化编程
当使用学习算法时,一段更快的代码通常意味着项目进展更快。例如,如果你的学习算法需要花费20分钟运行完成,这意味着你每个小时能“尝试”3个新主意。但是假如你的程序需要20个小时来运行,这意味着你一天只能“尝试”一个新主意,因为你需要花费这么长时间来等待程序的反馈。对于后者,假如你可以提升代码的效率让其只需要运行10个小时,那么你的效率差不多提升一倍。
矢量化编程是提高算法速度的一种有效方法。为了提升特定数值运算操作(如矩阵相乘、矩阵相加、矩阵-向量乘法等)的速度,数值计算和并行计算的研究人员已经努力了几十年。矢量化编程的思想就是尽量使用这些被高度优化的数值运算操作来实现我们的学习算法。
例如,假设 和
为向量,需要计算
,那么可以按以下方式实现(使用Matlab):
z = 0; for i=1:(n+1), z = z + theta(i) * x(i); end;
或者可以更加简单的写为:
z = theta' * x;
第二段程序代码不仅简单,而且运行速度更快。
通常,一个编写Matlab/Octave程序的诀窍是:
-
- 代码中尽可能避免显式的for循环。
上面的第一段代码使用了一个显式的for循环。通过不使用for循环实现相同功能,可以显著提升运行速度。对Matlab/Octave代码进行矢量化的工作很大一部分集中在避免使用for循环上,因为这可以使得Matlab/Octave更多地利用代码中的并行性,同时其解释器的计算开销更小。
关于编写代码的策略,开始时你会觉得矢量化代码更难编写、阅读和调试,但你需要在编码和调试的便捷性与运行时间之间做个权衡。因此,刚开始编写程序的时候,你可能会选择不使用太多矢量化技巧来实现你的算法,并验证它是否正确(可能只在一个小问题上验证)。在确定它正确后,你可以每次只矢量化一小段代码,并在这段代码之后暂停,以验证矢量化后的代码计算结果和之前是否相同。最后,你会有望得到一份正确的、经过调试的、矢量化且有效率的代码。
一旦对矢量化常见的方法和技巧熟悉后,你将会发现对代码进行矢量化通常并不太费劲。矢量化可以使你的代码运行的更快,而且在某些情况下,还简化了你的代码。
2、逻辑回归的向量化实现样例
我们想用批量梯度上升法对logistic回归分析模型进行训练,其模型如下:
让我们遵从公开课程视频与CS229教学讲义的符号规范,设 ,于是
,
,
为截距。假设我们有m个训练样本{(
,
) ,...,(
,
)},而批量梯度上升法的更新法则是:
,这里的
是对数似然函数,
是其导函数。
[注:下文的符号规范与<公开课程视频>或<教学讲义CS229:机器学习>中的相同,详细内容可以参见公开课程视频或教学讲义#1 http://cs229.stanford.edu/]
于是,我们需要如下计算梯度:
我们用Matlab/Octave风格变量x表示输入数据构成的样本矩阵,x(:,i)代表第 i个训练样本,x(j,i)就代表
(译者注:第i个训练样本向量的第j个元素)。同样,用Matlab/Octave风格变量y表示由训练样本集合的全体类别标号所构成的行向量,则该向量的第i个元素y(i)就代表上式中的
。(注意这里跟公开课程视频及CS229的符号规范不同,矩阵x按列而不是按行存放输入训练样本,同样,
是行向量而不是列向量。)
以下是梯度运算代码的一种实现,非常恐怖,速度极慢:
% 代码1 grad = zeros(n+1,1); for i=1:m, h = sigmoid(theta'*x(:,i)); temp = y(i) - h; for j=1:n+1, grad(j) = grad(j) + temp * x(j,i); end; end;
嵌套的for循环语句使这段代码的运行非常缓慢。以下是更典型的实现方式,它对算法进行部分向量化,带来更优的执行效率:
% 代码2 grad = zeros(n+1,1); for i=1:m, grad = grad + (y(i) - sigmoid(theta'*x(:,i)))* x(:,i); end;
但是,或许可以向量化得更彻底些。如果去除for循环,我们就可以显著地改善代码执行效率。特别的,假定b是一个列向量,A是一个矩阵,我们用以下两种方式来计算A*b:
% 矩阵-向量乘法运算的低效代码 grad = zeros(n+1,1); for i=1:m, grad = grad + b(i) * A(:,i); % 通常写法为A(:,i)*b(i) end; % 矩阵-向量乘法运算的高效代码 grad = A*b;
我们看到,代码2是用了低效的for循环语句执行梯度上升(译者注:原文是下降)运算,将b(i)看成(y(i) - sigmoid(theta'*x(:,i))),A看成x,我们就可以使用以下高效率的代码:
% 代码3 grad = x * (y- sigmoid(theta'*x));
这里我们假定Matlab/Octave的sigmoid(z)函数接受一个向量形式的输入z,依次对输入向量的每个元素施行sigmoid函数,最后返回运算结果,因此sigmoid(z)的输出结果是一个与z有相同维度的向量。
当训练数据集很大时,最终的实现(译者注:代码3)充分发挥了Matlab/Octave高度优化的数值线性代数库的优势来进行矩阵-向量操作,因此,比起之前代码要高效得多。
想采用向量化实现并非易事,通常需要周密的思考。但当你熟练掌握向量化操作后,你会发现,这里面有固定的设计模式(对应少量的向量化技巧),可以灵活运用到很多不同的代码片段中。
3、神经网络向量化
在本节,我们将引入神经网络的向量化版本。在前面关于神经网络介绍的章节中,我们已经给出了一个部分向量化的实现,它在一次输入一个训练样本时是非常有效率的。下边我们看看如何实现同时处理多个训练样本的算法。具体来讲,我们将把正向传播、反向传播这两个步骤以及稀疏特征集学习扩展为多训练样本版本。
Contents[hide] |
正向传播
考虑一个三层网络(一个输入层、一个隐含层、以及一个输出层),并且假定x是包含一个单一训练样本 的列向量。则向量化的正向传播步骤如下:
这对于单一训练样本而言是非常有效的一种实现,但是当我们需要处理m个训练样本时,则需要把如上步骤放入一个for循环中。
更具体点来说,参照逻辑回归向量化的例子,我们用Matlab/Octave风格变量x表示包含输入训练样本的矩阵,x(:,i)代表第个训练样本。则x正向传播步骤可如下实现:
% 非向量化实现 for i=1:m, z2 = W1 * x(:,i) + b1; a2 = f(z2); z3 = W2 * a2 + b2; h(:,i) = f(z3); end;
这个for循环能否去掉呢?对于很多算法而言,我们使用向量来表示计算过程中的中间结果。例如在前面的非向量化实现中,z2,a2,z3都是列向量,分别用来计算隐层和输出层的激励结果。为了充分利用并行化和高效矩阵运算的优势,我们希望算法能同时处理多个训练样本。让我们先暂时忽略前面公式中的b1和b2(把它们设置为0),那么可以实现如下:
% 向量化实现 (忽略 b1, b2) z2 = W1 * x; a2 = f(z2); z3 = W2 * a2; h = f(z3)
在这个实现中,z2,a2,z3都是矩阵,每个训练样本对应矩阵的一列。在对多个训练样本实现向量化时常用的设计模式是,虽然前面每个样本对应一个列向量(比如z2),但我们可把这些列向量堆叠成一个矩阵以充分享受矩阵运算带来的好处。这样,在这个例子中,a2就成了一个s2 X m的矩阵(s2是网络第二层中的神经元数,m是训练样本个数)。矩阵a2的物理含义是,当第i个训练样本x(:i)输入到网络中时,它的第i列就表示这个输入信号对隐神经元 (网络第二层)的激励结果。
在上面的实现中,我们假定激活函数f(z)接受矩阵形式的输入z,并对输入矩阵按列分别施以激活函数。需要注意的是,你在实现f(z)的时候要尽量多用Matlab/Octave的矩阵操作,并尽量避免使用for循环。假定激活函数采用Sigmoid函数,则实现代码如下所示:
% 低效的、非向量化的激活函数实现 function output = unvectorized_f(z) output = zeros(size(z)) for i=1:size(z,1), for j=1:size(z,2), output(i,j) = 1/(1+exp(-z(i,j))); end; end; end % 高效的、向量化激活函数实现 function output = vectorized_f(z) output = 1./(1+exp(-z)); % "./" 在Matlab或Octave中表示对矩阵的每个元素分别进行除法操作 end
最后,我们上面的正向传播向量化实现中忽略了b1和b2,现在要把他们包含进来,为此我们需要用到Matlab/Octave的内建函数repmat:
% 正向传播的向量化实现 z2 = W1 * x + repmat(b1,1,m); a2 = f(z2); z3 = W2 * a2 + repmat(b2,1,m); h = f(z3)
repmat(b1,1,m)的运算效果是,它把列向量b1拷贝m份,然后堆叠成如下矩阵:
这就构成一个s2 X m的矩阵。它和W1 * x相加,就等于是把W1 * x矩阵(译者注:这里x是训练矩阵而非向量, 所以W1 * x代表两个矩阵相乘,结果还是一个矩阵)的每一列加上b1。如果不熟悉的话,可以参考Matlab/Octave的帮助文档获取更多信息(输入“help repmat”)。rampat作为Matlab/Octave的内建函数,运行起来是相当高效的,远远快过我们自己用for循环实现的效果。
反向传播
现在我们来描述反向传播向量化的思路。在阅读这一节之前,强烈建议各位仔细阅读前面介绍的正向传播的例子代码,确保你已经完全理解。下边我们只会给出反向传播向量化实现的大致纲要,而由你来完成具体细节的推导(见向量化练习)。
对于监督学习,我们有一个包含m个带类别标号样本的训练集。 (对于自编码网络,我们只需令y(i) = x(i)即可, 但这里考虑的是更一般的情况。)
假定网络的输出有s3维,因而每个样本的类别标号向量就记为。在我们的Matlab/Octave数据结构实现中,把这些输出按列合在一起形成一个Matlab/Octave风格变量y,其中第i列y(:,i)就是y(i)。
现在我们要计算梯度项和
。对于梯度中的第一项,就像过去在反向传播算法中所描述的那样,对于每个训练样本(x,y),我们可以这样来计算:
在这里表示对两个向量按对应元素相乘的运算(译者注:其结果还是一个向量)。为了描述简单起见,我们这里暂时忽略对参数b(l)的求导,不过在你真正实现反向传播时,还是需要计算关于它们的导数的。
假定我们已经实现了向量化的正向传播方法,如前面那样计算了矩阵形式的变量z2, a2, z3和h,那么反向传播的非向量化版本可如下实现:
gradW1 = zeros(size(W1)); gradW2 = zeros(size(W2)); for i=1:m, delta3 = -(y(:,i) - h(:,i)) .* fprime(z3(:,i)); delta2 = W2'*delta3(:,i) .* fprime(z2(:,i)); gradW2 = gradW2 + delta3*a2(:,i)'; gradW1 = gradW1 + delta2*a1(:,i)'; end;
在这个实现中,有一个for循环。而我们想要一个能同时处理所有样本、且去除这个for循环的向量化版本。
为做到这一点,我们先把向量delta3和delta2替换为矩阵,其中每列对应一个训练样本。我们还要实现一个函数fprime(z),该函数接受矩阵形式的输入z,并且对矩阵的按元素分别执行。这样,上面for循环中的4行Matlab代码中每行都可单独向量化,以一行新的(向量化的)Matlab代码替换它(不再需要外层的for循环)。
在向量化练习中,我们要求你自己去推导出这个算法的向量化版本。如果你已经能从上面的描述中了解如何去做,那么我们强烈建议你去实践一下。虽然我们已经为你准备了反向传播的向量化实现提示,但还是鼓励你在不看提示的情况下自己去推导一下。
稀疏自编码网络
稀疏自编码网络中包含一个额外的稀疏惩罚项,目的是限制神经元的平均激活率,使其接近某个(预设的)目标激活率ρ。其实在对单个训练样本上执行反向传播时,我们已经考虑了如何计算这个稀疏惩罚项,如下所示:
在非向量化的实现中,计算代码如下:
% 稀疏惩罚Delta sparsity_delta = - rho ./ rho_hat + (1 - rho) ./ (1 - rho_hat); for i=1:m, ... delta2 = (W2'*delta3(:,i) + beta*sparsity_delta).* fprime(z2(:,i)); ... end;
但在上面的代码中,仍旧含有一个需要在整个训练集上运行的for循环,这里delta2是一个列向量。
作为对照,回想一下在向量化的情况下,delta2现在应该是一个有m列的矩阵,分别对应着m个训练样本。还要注意,稀疏惩罚项sparsity_delta对所有的训练样本一视同仁。这意味着要向量化实现上面的计算,只需在构造delta2时,往矩阵的每一列上分别加上相同的值即可。因此,要向量化上面的代码,我们只需简单的用repmat命令把sparsity_delta加到delta2的每一列上即可(译者注:这里原文描述得不是很清楚,看似应加到上面代码中delta2行等号右边第一项,即W2'*delta3上)。
中英文对照
- 矢量化 vectorization
-
- 逻辑回归 Logistic Regression
- 批量梯度上升法 batch gradient ascent
- 截距 intercept term
- 对数似然函数 the log likelihood
- 导函数 derivative
- 梯度 gradient
-
- 向量化 vectorization
- 正向传播 forward propagation
- 反向传播 backpropagation
- 训练样本 training examples
- 激活函数 activation function
- 稀疏自编码网络 sparse autoencoder
- 稀疏惩罚 sparsity penalty
- 平均激活率 average firing rate
-
-
- 链接:
-
http://deeplearning.stanford.edu/wiki/index.php/%E7%A5%9E%E7%BB%8F%E7%BD%91%E7%BB%9C%E5%90%91%E9%87%8F%E5%8C%96