逻辑回归的向量化实现样例
From Ufldl
我们想用批量梯度上升法对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按列而不是按行存放输入训练样本,同样,
是行向量而不是列向量。)
以下是梯度运算代码的一种实现,非常恐怖,速度极慢:
<span class="co1"><span style="color:#228b22;">% 代码1</span></span> grad = <span class="kw2"><span style="color:#00ff;">zeros</span></span><span class="br0"><span style="color:#0880;">(</span></span>n+<span class="nu0"><span style="color:#3333ff;">1</span></span>,<span class="nu0"><span style="color:#3333ff;">1</span></span><span class="br0"><span style="color:#0880;">)</span></span>; <span class="kw1"><span style="color:#00ff;">for</span></span> <span class="kw2"><span class="re0"><span style="color:#3333ff;">i</span></span></span>=<span class="nu0"><span style="color:#3333ff;">1</span></span>:m, h = sigmoid<span class="br0"><span style="color:#0880;">(</span></span>theta'*x<span class="br0"><span style="color:#0880;">(</span></span>:,<span class="kw2"><span class="re0"><span style="color:#3333ff;">i</span></span></span><span style="color:#0880;"><span class="br0">)</span><span class="br0">)</span></span>; temp = y<span class="br0"><span style="color:#0880;">(</span></span><span class="kw2"><span class="re0"><span style="color:#3333ff;">i</span></span></span><span class="br0"><span style="color:#0880;">)</span></span> - h; <span class="kw1"><span style="color:#00ff;">for</span></span> <span class="kw2"><span class="re0"><span style="color:#3333ff;">j</span></span></span>=<span class="nu0"><span style="color:#3333ff;">1</span></span>:n+<span class="nu0"><span style="color:#3333ff;">1</span></span>, grad<span class="br0"><span style="color:#0880;">(</span></span><span class="kw2"><span class="re0"><span style="color:#3333ff;">j</span></span></span><span class="br0"><span style="color:#0880;">)</span></span> = grad<span class="br0"><span style="color:#0880;">(</span></span><span class="kw2"><span class="re0"><span style="color:#3333ff;">j</span></span></span><span class="br0"><span style="color:#0880;">)</span></span> + temp * x<span class="br0"><span style="color:#0880;">(</span></span><span class="kw2"><span class="re0"><span style="color:#3333ff;">j</span></span></span>,<span class="kw2"><span class="re0"><span style="color:#3333ff;">i</span></span></span><span class="br0"><span style="color:#0880;">)</span></span>; <span class="kw1"><span style="color:#00ff;">end</span></span>; <span class="kw1"><span style="color:#00ff;">end</span></span>;
嵌套的for循环语句使这段代码的运行非常缓慢。以下是更典型的实现方式,它对算法进行部分向量化,带来更优的执行效率:
<span class="co1"><span style="color:#228b22;">% 代码2</span></span> grad = <span class="kw2"><span style="color:#00ff;">zeros</span></span><span class="br0"><span style="color:#0880;">(</span></span>n+<span class="nu0"><span style="color:#3333ff;">1</span></span>,<span class="nu0"><span style="color:#3333ff;">1</span></span><span class="br0"><span style="color:#0880;">)</span></span>; <span class="kw1"><span style="color:#00ff;">for</span></span> <span class="kw2"><span class="re0"><span style="color:#3333ff;">i</span></span></span>=<span class="nu0"><span style="color:#3333ff;">1</span></span>:m, grad = grad + <span class="br0"><span style="color:#0880;">(</span></span>y<span class="br0"><span style="color:#0880;">(</span></span><span class="kw2"><span class="re0"><span style="color:#3333ff;">i</span></span></span><span class="br0"><span style="color:#0880;">)</span></span> - sigmoid<span class="br0"><span style="color:#0880;">(</span></span>theta'*x<span class="br0"><span style="color:#0880;">(</span></span>:,<span class="kw2"><span class="re0"><span style="color:#3333ff;">i</span></span></span><span style="color:#0880;"><span class="br0">)</span><span class="br0">)</span><span class="br0">)</span></span>* x<span class="br0"><span style="color:#0880;">(</span></span>:,<span class="kw2"><span class="re0"><span style="color:#3333ff;">i</span></span></span><span class="br0"><span style="color:#0880;">)</span></span>; <span class="kw1"><span style="color:#00ff;">end</span></span>;
但是,或许可以向量化得更彻底些。如果去除for循环,我们就可以显著地改善代码执行效率。特别的,假定b是一个列向量,A是一个矩阵,我们用以下两种方式来计算A*b:
<span class="co1"><span style="color:#228b22;">% 矩阵-向量乘法运算的低效代码</span></span> grad = <span class="kw2"><span style="color:#00ff;">zeros</span></span><span class="br0"><span style="color:#0880;">(</span></span>n+<span class="nu0"><span style="color:#3333ff;">1</span></span>,<span class="nu0"><span style="color:#3333ff;">1</span></span><span class="br0"><span style="color:#0880;">)</span></span>; <span class="kw1"><span style="color:#00ff;">for</span></span> <span class="kw2"><span class="re0"><span style="color:#3333ff;">i</span></span></span>=<span class="nu0"><span style="color:#3333ff;">1</span></span>:m, grad = grad + b<span class="br0"><span style="color:#0880;">(</span></span><span class="kw2"><span class="re0"><span style="color:#3333ff;">i</span></span></span><span class="br0"><span style="color:#0880;">)</span></span> * A<span class="br0"><span style="color:#0880;">(</span></span>:,<span class="kw2"><span class="re0"><span style="color:#3333ff;">i</span></span></span><span class="br0"><span style="color:#0880;">)</span></span>; <span class="co1"><span style="color:#228b22;">% 通常写法为A(:,i)*b(i)</span></span> <span class="kw1"><span style="color:#00ff;">end</span></span>; <span class="co1"><span style="color:#228b22;">% 矩阵-向量乘法运算的高效代码</span></span> grad = A*b;
我们看到,代码2是用了低效的for循环语句执行梯度上升(译者注:原文是下降)运算,将b(i)看成(y(i) - sigmoid(theta'*x(:,i))),A看成x,我们就可以使用以下高效率的代码:
<span class="co1"><span style="color:#228b22;">% 代码3</span></span> grad = x * <span class="br0"><span style="color:#0880;">(</span></span>y- sigmoid<span class="br0"><span style="color:#0880;">(</span></span>theta'*x<span style="color:#0880;"><span class="br0">)</span><span class="br0">)</span></span>;
这里我们假定Matlab/Octave的sigmoid(z)函数接受一个向量形式的输入z,依次对输入向量的每个元素施行sigmoid函数,最后返回运算结果,因此sigmoid(z)的输出结果是一个与z有相同维度的向量。
当训练数据集很大时,最终的实现(译者注:代码3)充分发挥了Matlab/Octave高度优化的数值线性代数库的优势来进行矩阵-向量操作,因此,比起之前代码要高效得多。
想采用向量化实现并非易事,通常需要周密的思考。但当你熟练掌握向量化操作后,你会发现,这里面有固定的设计模式(对应少量的向量化技巧),可以灵活运用到很多不同的代码片段中。