标签归档:优化

CUDA, 软件抽象的幻影背后 之二

先更新到这儿,稍后再回来抛光查错。CUDA比较杂,我一写起来容易满嘴跑火车弄出错误,欢迎拍砖。

**********************************************************************

版权声明:原创作品,欢迎转载,但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明,否则将追究法律责任。

上一篇里说到,有两点对CUDA的计算能力影响甚大:数据并行,以及用多线程掩盖延迟。接下来我们要深入到其硬件实现,看一看这些机制是如何运作的。

通常人们经常说某GPU有几百甚至数千的CUDA核心,这很容易让人联想到多核CPU。不过事实上两种“核心”是不一样的概念,GPU的CUDA核心只相当于处理器中的执行单元,负责执行指令进行运算,并不包含控制单元。可以类比到CPU核心的是流多处理器(Streaming Multiprocessor,简写为SM. Kepler中叫做SMX,Maxwell中叫做SMM),通常一个GPU中有数个SM,而每个SM中包含几十或者上百个CUDA核心,以及数个warp scheduler(相当于控制单元)。如下图GM204中有16个SM,每个SM中有128个CUDA核心,4个warp scheduler。

GeForce_GTX_980_SM_Diagram-545x1024

图 1.  GM204的SM结构图

每个SM中有大量的寄存器资源,在GM204的例子中,有总共64k 32-bit寄存器,可以养活成千上万的线程。SM中另外一个重要资源是Shared Memory,没错,它正是软件抽象中Shared Memory的对应物。在GM204中,每个SM有96KB的Shared Memory.

到这里,SM在软件抽象里的对应也呼之欲出了,没错,正是Block。我们不妨先摆出这个对应:
Block <-> SM
Thread执行 <-> CUDA Cores
Thread数据 <-> Register/Local Memory

同一Grid下的不同Block会被分发到不同的SM上执行。SM上可能同时存在多个Block被执行,它们不一定来自同一个kernel函数。每个Thread中的局域变量被映射到SM的寄存器上,而Thread的执行则由CUDA核心来完成。

SM上可以同时存在多少个Block?这由硬件资源的消耗决定:每个SM会占用一定数量的寄存器和Shared Memory,因此SM上同时存活的Block数目不应当超过这些硬件资源的限制。由于SM上可以同时有来自不同kernel的Block存在,因此有时候即便SM上剩余资源不足以再容纳一个kernel A的Block,但却仍可能容纳下一个kernel B的Block.

接下来一个很重要的问题是Block如何被执行。我们可以看到,SM上的CUDA核心是有限的,它们代表了能够在物理上真正并行的线程数——软件抽象里,Block中所有的线程是并行执行的,这只是个逻辑上无懈可击的抽象,事实上我们不可能对一个任意大小的Block都给出一个同等大小的CUDA核心阵列,来真正并行的执行它们。
因而有了Warp这个概念:物理上,Block被划分成一块块分别映射到CUDA核心阵列上执行,每一块就叫做一个Warp.目前,CUDA中的Warp都是从threadIdx = 0开始,以threadIdx连续的32个线程为一组划分得到,即便最后剩下的线程不足32个,也将其作为一个Warp.CUDA kernel的配置中,我们经常把Block的size设置为32的整数倍,正是为了让它能够精确划分为整数个Warp(更深刻的原因和存储器访问性能有关,但这种情况下仍然和Warp的size脱不了干系)。
在GM204的SM结构图里我们可以看到,SM被划分成四个相同的块,每一块中有单独的Warp Scheduler,以及32个CUDA核心。Warp正是在这里被执行。
Warp的执行非常类似于SIMD. Warp中的活动线程由Warp Scheduler驱动,同步执行。我们可以看到,GM204中32个CUDA核心共享一个Warp Scheduler. 关于Warp执行中可能出现的复杂些的问题,留到下文另外说。

现在可以整理一下这个世界的图景了。SM上存活着几个Block,每个Block中的变量占据着自己的寄存器和Shared Memory,Block被划分成32个线程组成的Warp. 这样,大量的Warp生存在SM上,等待被调度到CUDA核心阵列去执行。

Warp Scheduler正如其名,是这个Warp世界里的调度者。当一个Warp执行中出现等待(存储器读写延迟等)后,Warp Scheduler就迅速切换到下一个可执行的Warp,对其发送指令直到这个Warp又一次出现等待,周而复始。这就是上一篇所说“用多线程掩盖延迟”在硬件图景下的模样。

CPU_GPU_COMPARE

图 2.  GPU用多个Warp掩盖延迟 / 与CPU计算模式的对比

本图引用自PPT “CUDA Overview” from Cliff Woolley, NVIDIA.

如图,GPU用多个Warp快速切换来掩盖延迟,而CPU用快速的寄存器来减小延迟。两者的重要区别是寄存器数目,CPU的寄存器快但少,因此Context Switch代价高;GPU寄存器多而慢,但寄存器数量保证了线程Context Switch非常快。

多少线程才能够掩盖掉常见的延迟呢?对于GPU,最常见的延迟大概要数寄存器写后读依赖,即一个局域变量被赋值后接着不久又被读取,这时候会产生大约24个时钟周期的延迟。为了掩盖掉这个延迟,我们需要至少24个Warp轮流执行,一个Warp遇到延迟后的空闲时间里执行其余23个Warp,从而保持硬件的忙碌。在Compute Capability 2.0,SM中有32个CUDA核心,平均每周期发射一条指令的情况下,我们需要24*32 = 768个线程来掩盖延迟。
保持硬件忙碌,用CUDA的术语来说,就是保持充分的Occupancy,这是CUDA程序优化的一个重要指标。

(未完待续)

***********************************************

一些后续补充。

网友邵:

SM上可以同时存在多少个Block,除了受到资源的限制之外,还受到设备上限的限制,每个SM有一个Device Limit,warps和blocks不能超过对应的上限。

Levenberg-Marquardt算法

匆匆一更。还没写完,对不起我太懒了T_T

版权声明:原创作品,欢迎转载,但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明,否则将追究法律责任。

Levenberg-Marquardt算法(下文简称LM算法)通常用于非线性最小二乘法的目标函数极小化。这是一个置信域方法(Trust-Region Method),为了防止步长太大而跳到非预期的局部极小值,这类算法自适应的调整步长。

假设f(x)为一个向量函数,每个分量都刻画了一个样本和模型预测之间的某种偏差,则非线性最小二乘目标函数为

(1)   \begin{eqnarray*} E&=&\transp{f}f \end{eqnarray*}

通常基于梯度的数值极小化算法,都是用一点的局域特征(如导数)建立该函数变化趋势的模型,以此来推测一个最可取的下降方向和步长。取一阶导数近似是线性模型,用切平面来近似函数在该点附近的变化趋势;取到二阶导数近似则是二阶模型,用一个抛物面来近似。当然,也可以取到更多阶,模型的几何直观也会越来越复杂,越来越精确的逼近真实的函数变化趋势,但是会付出更大的计算代价。当取到无穷阶时,模型在级数收敛区间内可以做到完全精确的刻画函数,此时这就是模型就是函数的泰勒展开了。

我们把f在任一点x附近展开f(x+\delta x)=f(x)+J(x)\delta x,这样Ex附近就是

(2)   \begin{eqnarray*} E(x+\delta x)&=&|f(x)+J(x)\delta x)|^2\\ &=&|f|^2+2\transp{f}J(x)\delta x + |J(x)\delta x|^2 \end{eqnarray*}

注意Jf而非E的导数,这里我们得到了E的二阶展开,即用一个抛物面来近似目标函数E在一点附近的形态。这个二阶模型显然是一个线性最小二乘问题——这正是我们所要的,把一个非线性最小二乘问题转换为一系列线性最小二乘问题来求解。

因此一个简单直接的想法就是,每一迭代步构造这样一个线性最小二乘问题,通过极小化它下降到下一个迭代点,循环直到收敛。但这个方案显然有问题。我们的二阶展开模型只能在一个小邻域上有效,当\delta x增大时,被我们忽略掉的高阶项就会越来越明显,最后让二阶模型彻底失效。这样,直接根据二阶模型计算出来的迭代点对目标函数而言可能并非我们想要的,它的函数值可能比上一点还要高(因而目标函数不降反增);也可能尽管比上一点下降了但是却跳到了另一个我们不想要的凸区域里(因而算法会收敛到离初始猜测解更远的局部极小值点上);更严重的情况,模型甚至可能是退化的,即一种极端形式的抛物面——就像一本被半卷起的书,此时我们在弯折的方向仍能找到极小值,但在其垂直方向已经无法找到极小值(或者说处处极小)。

为了解决以上问题,我们需要给模型加入一些约束,把迭代步长约束在一个合理的范围内,保证在这个范围内模型足够有效。同时,我们也希望在二阶模型退化的时候,约束能够拯救它,保证算法总能求出一个极小值。这个约束很简单,如下

(3)   \begin{eqnarray*} E(x+\delta x)&=&|f(x)+J(x)\delta x)|^2 + \lambda |C(x)\delta x|^2 \end{eqnarray*}

我们加入了一个二次约束项。单独看该项,它是一个中心在x点,开口向上的抛物面,或者说一个二次势阱。这个势阱像一个箍,把\delta x箍在x点附近的小邻域中。通过调节\lambda可以增强或减弱它对\delta x的约束:当\lambda充分大时,约束项成为主要贡献项,这时模型总有一个极小值,且步长较小;当\lambda较小时,约束项几乎可以忽略,算法以大步长快速往极小值接近。

再来看矩阵C(x). 这个矩阵通常被设置为对角的。此时,

(4)   \begin{eqnarray*} |C(x)\delta x|^2 &=& \transp{\delta x}\transp{C}C\delta x\\ &=&\sum_{i} \lambda_i \delta x_i^2 \end{eqnarray*}

\transp{C}C的本征值控制着每个方向上抛物面上升的快慢。可见,如果一个方向相应的本征值被设置的较小,那么该方向的约束项上升缓慢,约束较弱,步长跨度较大;否则反之。现实中使用的模型,目标函数对于各个参数的敏感度不同,有的参数略微调整就可以导致目标函数的强烈变化,有的则反之。C(x)这个参数矩阵使得我们可以针对各个不同的参数设置不同程度的步长约束,从而能够构造出一个稳定的算法。

(未完待补)