正者,正也。其心以为不然者,天门弗开矣。

作为《深度学习入门》的阅读笔记,本文简略概述深度学习的学习过程,损失函数,梯度法相关的内容,并辅以python实现

需要事先了解的

本文并不是从零开始的,需要实现了解

  • 蜻蜓点水python

    • 蜻蜓点水python_dlc
  • 简略感知机
  • 简易神经网络_推理和正向传播

不同的学习

首先,学习有不同的类型,人想,机器学习,深度学习;

然后,学习的基础是特征向量,机器学习中特征量还是人想出来的,深度学习则没有这一步骤,基本实现了全程机器管理
image-20251018220933182.png

损失函数

有道是能观测就能干涉,能干涉就能控制;那么应该观测什么变量才能知道学习的进度,从而干涉和控制学习的过程呢?这里直接给出答案:损失函数。损失函数理论上可以使用任何函数,这里使用均方误差交叉熵误差

这里有一个问题,为什么不使用精度(即准确度)来最为指标呢?首先说结论:在进行神经网络的学习时,不能将识别精度作为指标,并且对于大部分导数(后面会提到)为0的函数,都不能作为指标,会导致学习无法进行(无法更新参数),原因在了解学习的过程后自然能理解。

均方误差

直接公式,实现如下:
image-20251018223037143.png

  • y_k是前向传播的输出
  • t_k是监督数据的标签(一般是one-hot表示)
  • k是维度
import numpy as np
def mean_squared_error(y, t):
    return 0.5 * np.sum((y-t)**2)

print(mean_squared_error(np.array([0.1,0.7,0.2]),np.array([0,1,0]))) #判断正确
print(mean_squared_error(np.array([0.1,0.7,0.2]),np.array([1,0,0]))) #判断错误

#0.07000000000000002
#0.67

可以看出,如果判断正确,损失函数会是一个较小的数字,如果判断失败,损失函数会是一个较大的数字。

交叉熵误差

直接公式,实现如下:
image-20251018224310892.png

  • y_k是前向传播的输出
  • t_k是监督数据的标签(一般是one-hot表示)
  • k是维度
def cross_entropy_error(y, t):
    delta = 1e-7
    return -np.sum(t * np.log(y + delta))

print(cross_entropy_error(np.array([0.1,0.7,0.2]),np.array([0,1,0])))
print(cross_entropy_error(np.array([0.1,0.7,0.2]),np.array([1,0,0])))

#0.3566748010815999
#2.302584092994546
  • 其中delta是为了避免log的参数有0的情况

Mini-Batch学习

和推理的过程一样,学习也可以用批处理的方法完成,不过实现上有些区别;这里以交叉熵误差为例

第一个是损失函数公式的变化:
image-20251018230403347.png

可见其实就是取平均,实现如下:

def cross_entropy_error2(y, t):
    batch_size = 1
    if y.ndim > 1:
        batch_size = y.shape[0]
    return -np.sum(t * np.log(y + 1e-7)) / batch_size

#测试不同维度
print(cross_entropy_error2(np.array([0.1,0.7,0.2]),np.array([0,1,0])))
print(cross_entropy_error2(np.array([0.1,0.7,0.2]),np.array([1,0,0])))
print(cross_entropy_error2(np.array([[0.1,0.7,0.2],[0.1,0.7,0.2]]),np.array([[0,1,0],[0,1,0]])))
print(cross_entropy_error2(np.array([[0.1,0.7,0.2],[0.1,0.7,0.2]]),np.array([[1,0,0],[1,0,0]])))

第二个是训练数据其实还是很多的,想一次性全部学习是高耗能甚至有些不现实的,所以需要使用np.random.choice抽取一个batch来学习

# 这里的路径按本地实际来
from dfs.dataset.mnist import load_mnist
import sys,os
sys.path.append(os.pardir)

(x_train, t_train), (x_test, t_test) = load_mnist(flatten=True,
normalize=False,one_hot_label=True)

# 输出各个数据的形状
print(x_train.shape) # (60000, 784)
# 输出各个数据的形状
print(t_train.shape) # (60000,)
print(x_test.shape) # (10000, 784)
print(t_test.shape) # (10000,)

train_size = x_train.shape[0]
batch_size = 10
batch_mask = np.random.choice(train_size, batch_size)
print(batch_mask) # 输出抽选的索引
x_batch = x_train[batch_mask]
t_batch = t_train[batch_mask]

print(x_batch.shape) # (10, 784)
print(t_batch.shape) # (10,)


#(60000, 784)
#(60000, 10)
#(10000, 784)
#(10000, 10)
#[ 6681  4809 12952  4669  1444 20997 56685 39195 24206  8649]
#(10, 784)
#(10, 10)

第三如果标签不是one-hot表示,有那么一个特殊的实现方法,思路是虽然是求和,但只有正确项会有值

def cross_entropy_error2(y, t):
    if y.ndim == 1:
        t = t.reshape(1, t.size)
        y = y.reshape(1, y.size)
    batch_size = y.shape[0]
    print(y[np.arange(batch_size), t])
    return -np.sum(np.log(y[np.arange(batch_size), t] + 1e-7)) / batch_size

print(cross_entropy_error2(np.array([[0.1,0.7,0.2],[0.1,0.7,0.2]]),np.array([1,2])))
#[0.7 0.2]
#0.9830561067579126

另外有个概念会在后面用到这里先提出:

epoch是一个单位。一个epoch表示学习中所有训练数据均被使用过一次时的更新次数。比如,对于 10000笔训练数据,用大小为 100笔数据的mini-batch进行学习时,重复随机梯度下降法100次,所有的训练数据就都被“看过”了 。此时,100次就是一个epoch。

导数和梯度

有了损失函数,就可以观测学习的进度,但是如何在此基础上对其进行控制呢?要做到这点需要解决三个问题,函数的值是会通过自变量的变化而变化的,改哪个变量会使损失函数变小?如何改?要改多少?

直接说结论:修改的变量是权重和偏置,依靠导数计算出的梯度来修改自变量,改多少依赖于人工设定的学习率

导数

不准确的说,导数就是x周围很小范围的斜率(中心差分),既不是x左边的斜率也不是x右边的斜率(前向差分);数值微分和真实的解析解必然是有所差距的。

求导数的实现如下:

def numerical_diff(f,x,h=1e-4):
    return (f(x+h)-f(x-h))/(2*h)
def func_test(x):
    return 2*x+1

numerical_diff(func_test,np.array([2,3]))

偏导数和梯度

多个变量的函数的导数就是偏导数,当然一个偏导数只能对其中一个变量来做,另外的变量就当作常数对待,称为函数对该变量的偏导数,这是解析解的做法;对于数值解则更加的简单一点,固定其他变量,只变化一个变量来求导数就是函数针对这个变量的偏导数。

把所有变量的偏导数组成向量,就是梯度,梯度指向的方向是函数变化(增加)率最大的方向。所以对其取负数,就是损失函数下降最快的方向。

实现如下:

def numerical_gradient(f, x):
    h = 1e-4 # 0.0001
    grad = np.zeros_like(x)

    it = np.nditer(x, flags=['multi_index'], op_flags=['readwrite'])
    while not it.finished:
        idx = it.multi_index
        tmp_val = x[idx]
        x[idx] = float(tmp_val) + h
        fxh1 = f(x) # f(x+h)

        x[idx] = tmp_val - h
        fxh2 = f(x) # f(x-h)
        grad[idx] = (fxh1 - fxh2) / (2*h)

        x[idx] = tmp_val # 还原值
        it.iternext()

    return grad

def test_func(x):
    return x[0]*x[0]+x[1]*x[1]

numerical_gradient(test_func,np.array([3.0,4.0]))
#array([6., 8.])

梯度下降法

按上述说的,按照梯度反方向更新参数就是梯度下降法,公式如下:

image-20251019112529214.png

其中偏导数前面的常数就是学习率,人为指定,规定学习的程度,过大过小都不好,这里默认0.01,实现如下:

def gradient_descent(f,init_x,lr=0.01,step_num=100):
    x = init_x

    for i in range(step_num):
        grad = numerical_gradient(f,x)
        x -= lr*grad
    return x

def function_2(x):
    return x[0]*x[0]+x[1]*x[1]

init_x=np.array([-3.0,4.0])


gradient_descent(function_2,init_x,lr=0.1,step_num=100)
#array([-6.11110793e-10,  8.14814391e-10])

另外,像这种人为指定的参数,称为超参数,像是mini-batch的大小,输入的形状,学习的次数都是超参数

总结和实现

现在就可以继续上章,完成MINIST的学习步骤了,但是还有一些细节需要说明

神经网络的梯度法

在实际的神经网络中要求的是损失函数对于权重的梯度,也就是对权重矩阵中的每个成员求梯度,最后结果还是矩阵:

image-20251019120202053.png

对于实现而言,只要能保证矩阵的形状没有问题即可

学习算法的步骤

image-20251019120525911.png

步骤4则是重复1-3的步骤

而这种随机选取+计算梯度的方法称为:随机梯度下降法(SGD)

实现

先列出代码

# 这里的路径按本地实际来
import numpy as np
from dfs.dataset.mnist import load_mnist
import sys,os
sys.path.append(os.pardir)

(x_train, t_train), (x_test, t_test) = load_mnist(flatten=True,
normalize=False,one_hot_label=True)

# 输出各个数据的形状
#print(x_train.shape) # (60000, 784)
# 输出各个数据的形状
#print(t_train.shape) # (60000,10)
#print(x_test.shape) # (10000, 784)
#print(t_test.shape) # (10000,10)

train_size = x_train.shape[0]
batch_size = 10
batch_mask = np.random.choice(train_size, batch_size)
#print(batch_mask)
x_batch = x_train[batch_mask]
t_batch = t_train[batch_mask]

print(x_batch.shape) # (10, 784)
print(t_batch.shape) # (10,10)


def cross_entropy_error(y, t):
    batch_size = 1
    if y.ndim > 1:
        batch_size = y.shape[0]
    return -np.sum(t * np.log(y + 1e-7)) / batch_size

def softmax(a):
    c = np.max(a)
    exp_a = np.exp(a-c) #溢出对策
    sum_exp_a = np.sum(exp_a)
    y = exp_a /sum_exp_a
    return y

def sigmoid(x):
    return 1/(1+np.exp(-x))


class SigMoidLayer:
    def __init__(self,w,b,s="SigMoid"):
        self.layer_name = s
        self.w=w
        self.b = b
    def forward(self,x):
        a=np.dot(x,self.w)+self.b
        return sigmoid(a)

class OutputLayer2:
    def __init__(self,w,b,s="Ouput"):
        self.layer_name = s
        self.w=w
        self.b = b
    def forward(self,x):
        a=np.dot(x,self.w)+self.b
        return softmax(a)

def numerical_gradient(f, x):
    h = 1e-4 # 0.0001
    grad = np.zeros_like(x)

    it = np.nditer(x, flags=['multi_index'], op_flags=['readwrite'])
    while not it.finished:
        idx = it.multi_index
        tmp_val = x[idx]
        x[idx] = float(tmp_val) + h
        fxh1 = f(x) # f(x+h)

        x[idx] = tmp_val - h
        fxh2 = f(x) # f(x-h)
        grad[idx] = (fxh1 - fxh2) / (2*h)

        x[idx] = tmp_val # 还原值
        it.iternext()

    return grad

class TestMNIST:
    # 正常情况下,权重或者偏置都应该是学习的结果,这里直接赋值只是一种假设
    def __init__(self):
        self.rate=0.01
        self.grads_b = {}
        self.grads_w = {}
        self.network_w = {}
        self.network_b = {}
        self.network_w['L1'] = np.random.randn(784,50)
        self.network_b['L1'] = np.zeros(50)
        self.network_w['L2'] = np.random.randn(50,100)
        self.network_b['L2'] = np.zeros(100)

        # 输出层
        self.network_w['L3'] = np.random.randn(100,10)
        self.network_b['L3'] = np.zeros(10)

        self.layers=['L1','L2']

    def forward(self,x):
        tmp=x
        for l in self.layers:
            tmp = SigMoidLayer(self.network_w[l],self.network_b[l],l).forward(tmp)
            #print(l,tmp)
        tmp = OutputLayer2(self.network_w['L3'],self.network_b['L3'],'L3').forward(tmp)
        return tmp

    def loss(self,x,t):
        return cross_entropy_error(self.forward(x),t)

    def grad(self,x,t):
        loss_W = lambda W: self.loss(x, t)

        self.grads_w={}
        self.grads_b={}

        self.grads_w['L1']=numerical_gradient(loss_W,self.network_w['L1'])
        self.grads_b['L1']=numerical_gradient(loss_W,self.network_b['L1'])
        print('L1 Done.')

        self.grads_w['L2']=numerical_gradient(loss_W,self.network_w['L2'])
        self.grads_b['L2']=numerical_gradient(loss_W,self.network_b['L2'])
        print('L2 Done.')
        self.grads_w['L3']=numerical_gradient(loss_W,self.network_w['L3'])
        self.grads_b['L3']=numerical_gradient(loss_W,self.network_b['L3'])
        print('L3 Done.')

    def update(self):
        self.network_w['L1']-=self.rate*self.grads_w['L1']
        self.network_b['L1']-=self.rate*self.grads_b['L1']
        self.network_w['L2']-=self.rate*self.grads_w['L2']
        self.network_b['L2']-=self.rate*self.grads_b['L2']
        self.network_w['L3']-=self.rate*self.grads_w['L3']
        self.network_b['L3']-=self.rate*self.grads_b['L3']



nett = TestMNIST()
#测试代码
#nett.forward(np.random.randn(100,784))

# 学习过程
mini_batch_size = x_train.shape[0]
select_size = 100
loss_point =[]
train_time=10000

for i in range(train_time):
    # 选择数据
    mask = np.random.choice(mini_batch_size, select_size)
    x_batch = x_train[mask]
    t_batch = t_train[mask]
    print(x_batch.shape,t_batch.shape)
    print("Loop",i+1)
    # 损失函数
    loss_p=nett.loss(x_batch,t_batch)
    loss_point.append(loss_p)
    print("Loss:",loss_p)
    # 计算梯度
    nett.grad(x_batch,t_batch)
    # 更新参数
    nett.update()



#推理过程
batch_size=100
accuracy_cnt=0

for i in range(0,len(x_test),batch_size):
    x_batch = x_test[i:i+batch_size]
    t_batch = t_test[i:i+batch_size]

    y_batch = nett.forward(x_batch)
    #print(y_batch.shape)
    p=np.argmax(y_batch,axis=1)
    t=np.argmax(t_batch,axis=1)
    accuracy_cnt+=(np.sum(p==t)/batch_size)

print("Accuracy:" + str(float(accuracy_cnt) / (len(x_test)/batch_size)))


import matplotlib.pyplot as plt
 # 生成数据
xx = np.arange(0, train_time)
# 绘制图形
plt.plot(xx, loss_point)
plt.show()
  • 首先要注意的是矩阵的形状一定要“严丝合缝”,这是做矩阵计算的基础

    • np.zeros(50)是创建0矩阵的方法,第一个参数是shape
  • 然后当运行该代码时可以发现一点,就是随机梯度下降的处理速度十分缓慢,所以在之后会需要反向传播这一方法来加速梯度的计算

评价方法

评价一个神经网络的方法就是评价它的泛化能力,也就是对于陌生数据的识别精度:

  • 如果只对训练数据识别精度较好,则是发生了过拟合
  • 一般来说,当学习经过一个epoch,就可以对测试数据和训练数据进行推理,然后进行可视化对比,这里只给出伪代码:
    image-20251019192442837.png

image-20251019192500604.png

epoch是一个单位。一个epoch表示学习中所有训练数据均被使用过一次时的更新次数。比如,对于 10000笔训练数据,用大小为 100笔数据的mini-batch进行学习时,重复随机梯度下降法100次,所有的训练数据就都被“看过”了。此时,100次就是一个epoch。

不幸的结果

和原著对比,这里的效果差了很多(最重要的是运行时间巨长!!!!整整2个星期),原因会放在未来探索,可能的原因一方面是初始化参数的问题,一方面可能是层级结构的问题。准确度在:Accuracy:0.7011999999999997 而不是0.9以上,梯度下降的也不太对,如下图。

但不管怎么样,SGD是有效果的,接下来需要解决其速度太慢的问题。
image-20251110204818780.png

image-20251103191147939.png

化其万物而不知其禅之者,焉知其所终?焉知其所始?正而待之而已耳。

作为《深度学习入门》的阅读笔记,本文简略概述正向传播和激活函数的内容,并辅以python实现

需要事先了解的

本文并不是从零开始的,需要实现了解

  • 蜻蜓点水python

    • 蜻蜓点水python_dlc
  • 简略感知机

学习和推理

和求解机器学习问题的步骤(分成学习和推理两个阶段进行)一样,使用神经网络解决问题时,也需要首先使用训练数据(学习数据)进行权重参数的学习;进行推理时,使用刚才学习到的参数,对输入数据进行分类。

首先,通过机器学习解决问题通常分为两个阶段:学习推理

  • 对于感知机或者神经网络,学习的过程就是调整参数的过程,是需要事先执行的步骤;
  • 在完成学习后,使用确定的参数解决实际问题的过程,称为推理;在神经网络或感知机中,通过前向传播来实现;

从感知机到神经网络

回顾下异或问题,在感知机的实现中,输入先通过与与非门或门的处理,再使用其输出经过与门的处理就完成了异或的操作:

image-20251009221438817.png

这种多层感知机就是一种二层神经网络层数的算法和感知机的层数算法一致,只是表示法不同:

image-20251018113047122.png

先解释结构,x1 x2的位置称为输入层,y所在位置称为输出层,两层之间的不管有多少层都是中间层或者隐藏层

不同点在于门的运算过程都跑到上去了,回顾感知机的实现,所谓的实现就是不同权重的调整,所以边上既包含了参数也包含了固定的传播算法,节点或者称神经元上包含了输出输入的结果,但也包含了一些其他处理过程,如激活函数

一般而言,“朴素感知机”是指单层网络,指的是激活函数使用了阶跃函数(后述)的模型。“多层感知机”是指神经网络,即使用sigmoid函数(后述)等平滑的激活函数的多层网络。

传播算法和激活函数

首先,回顾下感知机中的权重和偏置和感知机的实现:

image-20251009215857997.png

简单来说感知机的实现就是该公式的实现,其中b是被称为偏置的参数,用于控制神经元被激活的容易程度;而w1和w2是表示各个信号的权重的参数,用于控制各个信号的重要性。 而b+w1x1+w2x2就是层和层之间的传播算法,从下图可以看出,偏置也可作为神经元为1而权重为b进行传播:

image-20251018114534693.png

解决了传播问题,再看一眼公式可以发现0和1的判断是依赖于b+w1x1+w2x2的输出的,而这0和1的判断规则也不一定是以公式中的方式出现的,是可以记作函数h(x)来表示的,这个函数就是激活函数,激活函数的输出才是神经元最后表示的输出。

简而言之:激活函数,是用来决定何时激活信号的函数

image-20251018115114656.png

阶跃函数

就从感知机使用的激活函数来看:

image-20251018115844401.png

该函数称为阶跃函数;实现和图像如下,注意到这并不是一个平滑的函数

def step_function(x):
    y = x > 0
    return y.astype(int)
x = np.array([-1.0,1.0,2.0])
print(x)
print(x>0)
y = step_function(x)
print(type(y))

#[-1.  1.  2.]
#[False  True  True]
#<class 'numpy.ndarray'>

import matplotlib.pyplot as plt
x = np.arange(-5.,5.,0.1)
y = step_function(x)
plt.plot(x,y)
plt.ylim(-0.1,1.1)
plt.show()

image-20251018120028301.png

Sigmoid函数

Sigmoid函数是神经网络经常使用一个函数,公式,实现和图像如下:

image-20251018121412781.png

def sigmoid(x):
    return 1/(1+np.exp(-x))

x = np.arange(-5.,5.,0.1)
y = sigmoid(x)
plt.plot(x,y)
plt.ylim(-0.1,1.1)
plt.show()

image-20251018141612162.png

函数对比

image-20251018121517952.png

上面说到,阶跃函数并不是平滑的函数,而sigmoid是,然后二者都是非线性函数,这里先给出结论:

  • 神经网络的激活函数必须使用非线性函数。

原因是由线性函数的性质造成的,如果把中间层合并在一起,所做的处理更加类似于一个超级复杂的复合函数,但是线性函数再怎么复合,都可以找到一个参数来等价的表示这些复合的操作,对于中间层来说就是可以用一层来代表N层,这样中间层的存在意义就没有了,神经网络也就无法发挥多层网络带来的优势。下图为简单的举例:

scanner_20251018_122457.jpg

ReLU函数

sigmoid函数很早就开始被使用了,而最近则主要使用ReLU(Rectified Linear Unit)函数。

公式,实现和图像如下:

image-20251018141008571.png

def relu(x):
    return np.maximum(0,x)
x = np.arange(-5.,5.,0.1)
y = relu(x)
plt.plot(x,y)

image-20251018141528425.png

神经网络和矩阵

矩阵乘法

回顾一下numpy的相关操作,shape表示数组形状dim表示数组维数;多维数组可以代表矩阵,横为行竖为列;

image-20251018150234147.png

然后简单了解下矩阵乘法

image-20251018150255344.png

代码如下:

A=np.array([[1,2,3],[4,5,6]])
B=np.array([[1,2],[3,4],[5,6]])
print(A.shape,B.shape)

C=np.dot(A,B)
print(C)
print(C.shape)

#(2, 3) (3, 2)
#[[22 28]
# [49 64]]
#(2, 2)

在矩阵的乘积运算中,对应维度的元素个数要保持一致 :

image-20251018151409921.png

神经网络和矩阵

从下图可见,之前所说层与层之间的传播算法,可以统统用矩阵乘法表示。而这一层接一层的传播就是所谓的(前)正向传播

image-20251018151655115.png

前向传播

前向传播用语言表达比较繁琐,直接看图和代码则比较方便,无非是一层接一层的矩阵计算和激活函数计算:
image-20251018153304137.png

image-20251018153324040.png

输出层的激活函数

激活函数已经了解,但是输出层的激活函数稍微有些不同,但也只是输出函数类型的问题。

首先 ,对于机器学习解决的问题有两种,一种是所谓分类问题,就好比数字识别,从图片判断是0-9的哪一个,还有一种是回归问题,好比从历史数据判断今天地铁的客流量,一个是选择题,一个是填空题。

在输出层,对于选择题使用softmax函数,对于回归问题使用恒等函数

顺带一提,输出层神经元的数量在分类问题中取决于分类的个数。

恒等函数

恒等函数图像如下,其实就是将输入原样输出,这里便不列出代码identity_function

image-20251018154912288.png

SoftMax函数

对于分类使用的SoftMax函数则有一些说法,先看公式和基本实现:

image-20251018155139122.png

首先这涉及了前一层的所有输入,类比sql语句,这便成为全连接,意思是每一个输出都完整的接受了上一个输入层中的所有数据;而计算时用所有输入做分母,其中一个做分子,这又有点概率的意思,事实上也是这样

最后,SoftMax函数一般只用于学习阶段,因为在推理过程中,一般而言只把输出值最大的神经元所对应的类别作为识别结果,使用softmax会造成额外的计算量

下面是代码实现:

def softmax(a):
    exp_a = np.exp(a)
    sum_exp_a=np.sum(exp_a)
    y=exp_a/sum_exp_a
    return y

a = np.array([0.3,2.9,4.0])
print(softmax(a))

#下两行会引发溢出异常
a = np.array([1010,1000,990])
print(softmax(a))

这样是用问题的,因为对于指数函数有溢出的风险,但如果按如下步骤进行改造则可以解决这个问题:

image-20251018185006687.png

其中 C'一般取a_n这个数组中的最大值取反,这样削减一下大小,就不会有溢出的问题了:

# softmax的改进变形
def softmax(a):
    c = np.max(a)
    exp_a = np.exp(a-c) #溢出对策
    sum_exp_a = np.sum(exp_a)
    y = exp_a /sum_exp_a

    return y

a = np.array([1010,1000,990])
print(softmax(a))

#[9.99954600e-01 4.53978686e-05 2.06106005e-09]

总结和阶段性实现

目前,一个支持前向传播的神经网络便可以实现,但是目前只能假设一个权重参数(因为没有学习阶段的支持),然后通过矩阵乘法和中间层和输出层的激活函数来实现前向传播。

所以目前以MNIST数据集为例来实现一个简易的前向传播

测试数据和训练数据

机器学习中,一般将数据分为训练数据(监督数据 )测试数据两部分来进行学习和实验等。

  • 首先,使用训练数据进行学习,寻找最优的参数;
  • 然后,使用测试数据评价训练得到的模型的实际能力

在处理测试数据过程中,如果发现效果很好,则说明该模型的泛化能力很好,即面对未知数据表现很好;反面的,如果只对训练数据表现好,则说明泛化能力差,这种现象称为过拟合

MNIST数据集的读取

读取该数据集需要一些工具函数的辅助,这里直接使用《深度学习入门》中提供的工具实现:

# 这里的路径按本地实际来
from dfs.dataset.mnist import load_mnist
import sys,os
sys.path.append(os.pardir)

(x_train, t_train), (x_test, t_test) = load_mnist(flatten=True,
normalize=False,one_hot_label=False)

# 输出各个数据的形状
print(x_train.shape) # (60000, 784)
# 输出各个数据的形状
print(t_train.shape) # (60000,)
print(x_test.shape) # (10000, 784)
print(t_test.shape) # (10000,)
  • train代表训练数据,test表示测试数据
  • x代表输入的图片,t代表最终的标签(在one_hot_label是false时会在显示标签0-9,在true时会返回一个10位数组,在标签数字对应的位置是1,其余为0)
from PIL import Image

def img_show(img):
    pil_img = Image.fromarray(np.uint8(img))
    pil_img.show()

img = x_train[0]
label = t_train[0]
print(label) # 5
print(img.shape) # (784,)
img = img.reshape(28, 28) # 把图像的形状变成原来的尺寸
print(img.shape) # (28, 28)
img_show(img)
  • 该数据集中,输入的是一堆图片和图片对应的标签,在这里将图片展开为一维数组(flatten=True),并使用原始像素信息,而不是将0-255的颜色值压缩到0-1的范围(normalize=false),而将数据压缩到0-1范围这种操作称为标准化或正交化
one-hot表示是仅正确解标签为1,其余皆为0的数组,就像[0,0,1,0,0,0,0,0,0,0]这样。当one_hot_label为False时,只是像 7、 2这样简单保存正确解标签;当 one_hot_label为 True时,标签则保存为one-hot表示。

实现

首先,由于并没有学习这一步骤,所以会给所有权重和偏置设置符合正态分布的默认值;

然后因为传播的实现就是矩阵乘法的实现,现在假设使用的是3层神经网络,则三个权重矩阵的形状要符合矩阵乘法的要求,所以假设矩阵形状如下:
image-20251018201706760.png

最后,由于一个图片一个图片处理过于缓慢,而矩阵乘法却能吃到SIMD的优化福利,所以会使用批处理的方法,即如下图所示,一次性打包处理100张图像,为此,可以把x的形状改为100 × 784,将100张图像打包作为输入数据:
image-20251018201945559.png

如此实现如下,首先是神经网络的实现和使用

import numpy as np
def softmax(a):
    c = np.max(a)
    exp_a = np.exp(a-c) #溢出对策
    sum_exp_a = np.sum(exp_a)
    y = exp_a /sum_exp_a
    return y

def sigmoid(x):
    return 1/(1+np.exp(-x))

class SigMoidLayer:
    def __init__(self,w,b,s="SigMoid"):
        self.layer_name = s
        self.w=w
        self.b = b
    def forward(self,x):
        a=np.dot(x,self.w)+self.b
        return sigmoid(a)

class OutputLayer2:
    def __init__(self,w,b,s="Ouput"):
        self.layer_name = s
        self.w=w
        self.b = b
    def forward(self,x):
        a=np.dot(x,self.w)+self.b
        return softmax(a)
class TestMNIST:
    # 正常情况下,权重或者偏置都应该是学习的结果,这里直接赋值只是一种假设
    def __init__(self):
        self.network_w = {}
        self.network_b = {}
        self.network_w['L1'] = np.random.randn(784,50)
        self.network_b['L1'] = 1
        self.network_w['L2'] = np.random.randn(50,100)
        self.network_b['L2'] = 1

        # 输出层
        self.network_w['L3'] = np.random.randn(100,10)
        self.network_b['L3'] = 1

        self.layers=['L1','L2']

    def forward(self,x):
        tmp=x
        for l in self.layers:
            tmp = SigMoidLayer(self.network_w[l],self.network_b[l],l).forward(tmp)
            print(l,tmp)
        tmp = OutputLayer2(self.network_w['L3'],self.network_b['L3'],'L3').forward(tmp)
        return tmp


nett = TestMNIST()
#测试代码
#nett.forward(np.random.randn(100,784))

batch_size=100

for i in range(0,len(x_test),batch_size):
    x_batch = x_test[i:i+batch_size]
    #print(x_batch.shape)
    y_batch = nett.forward(x_batch)
    p=np.argmax(y_batch,axis=1)
    print(p)
#[8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8
# 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8
# 8 8 8 8 8 8 6 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8]

可见结果并不好,都只会输出8,即使如此也要测一下准确率:

accuracy_cnt=0

for i in range(0,len(x_test),batch_size):
    x_batch = x_test[i:i+batch_size]
    y_batch = nett.forward(x_batch)
    p=np.argmax(y_batch,axis=1)
    accuracy_cnt += np.sum(p == t_test[i:i+batch_size])
print("Accuracy:" + str(float(accuracy_cnt) / len(x_test)))
#Accuracy:0.0965

太好了,基本都对不上。

通于一而万事毕,无心得而鬼神服。

本文以数组加法为例,简述CUDA的编程模型

基本概念

确保了解《CUDA的线程结构》中的内容再进行下去,因为对于线程结构本章会省略一部分描述

编程模型

首先,编程模型是写程序时可被控制的部分,通过编程模型可以控制计算设备的行为;其次,编程模型对于编译器和库函数来说是一层通信抽象;

再次,对于GPU来说,编程模型可分为

  • 核函数
  • 内存管理
  • 线程管理

最后,对于Cuda编程,可抽象出几个过程

  • 领域层(业务主导要实现的目标)
  • 逻辑层(基于逻辑组织线程代码)
  • 硬件层(思考如何将线程映射到机器上,以及思考如何提升性能)

接下来,会从以下方面说明Cuda的编程结构

  • 内存
  • 线程
  • 核函数

    • 启动
    • 编写
    • 验证
  • 错误处理

内存管理

首先应该理解异构环境下内存的样貌:通过PCIe总线分割和通信的两种内存:主机内存设备内存

一个异构环境,通常有多个CPU多个GPU,他们都通过PCIe总线相互通信,也是通过PCIe总线分隔开的。所以我们要区分一下两种设备的内存:

  • 主机:CPU及其内存
  • 设备:GPU及其内存

注意这两个内存从硬件到软件都是隔离的(CUDA6.0 以后支持统一寻址),我们目前先不研究统一寻址,我们现在还是用内存来回拷贝的方法来编写调试程序,以巩固大家对两个内存隔离这个事实的理解。

CUDA的内存管理是C风格的,列表如下:

标准C函数CUDA C 函数说明
malloccudaMalloc内存分配
memcpycudaMemcpy内存复制
memsetcudaMemset内存设置
freecudaFree释放内存

另外,CUDA中设备内存也是分层次的,分为单网格中公用的全局内存和单线程块中公用的共享内存,大致如图:

3.png

线程管理

在《CUDA的线程结构》中,曾经用草图展示过cuda线程的组织形式,这里用更好的图展示一下:

4.png

另外,不同线程块的线程物理隔离无法相互影响,同一个线程块中的线程可以实现以下两种操作:

  • 同步
  • 共享内存

最后再重申一下《CUDA的线程结构》中提过的,如何标识一个线程

  • blockIdx(线程块在线程网格内的位置索引)
  • threadIdx(线程在线程块内的位置索引)

这两个内置结构体基于 uint3 定义,包含三个无符号整数xyz的结构

还有blockIdx中三个字段的范围和threadIdx中三个字段的范围:

  • gridDim
  • blockDim

dim3类型(基于uint3定义的数据结构)的变量,也包含三个字段xyz

另外,线程块在3个轴上的个数和一个线程块的线程数是有限制的:

  • 线程块在xyz三轴上分别最多:2^31-1,65535,65535
  • 一个线程块线程总数最多1024,在三维上最多:1024,1024,64

核函数

之前hello world时已经用到了核函数,但是没有给一个明确定义,这里给出:

核函数就是在CUDA模型上诸多线程中运行的那段串行代码,这段代码在设备上运行,用NVCC编译,产生的机器码是GPU的机器码,所以我们写CUDA程序就是写核函数

  • 第一步我们要确保核函数能正确的运行产生正切的结果
  • 第二步优化CUDA程序的部分,无论是优化算法,还是调整内存结构,线程结构都是要调整核函数内的代码,来完成这些优化的。

接下来将以数组加法为例来实现和使用核函数,以及应用一些内存管理的函数

单机版数组加法

无需多言,直接上代码,明确我们要干什么:

/**
 * 单机版:两数组相加
 *
*/
#include<iostream>

const double a=1.23;
const double b=3.45;
const double c=4.68;
const int N = 100000000;
const double EPSILON = 1.0e-15;

// 函数1: 数组相加(数组默认是动态分配的)
void add(const double *x,const double *y,double *z,const int N) {
    for (int i=0;i<N;i++) {
        z[i] = x[i]+y[i];
    }
}

// 函数2: 检查结果
void check(const double *z,const int N) {
    bool has_error = false;
    for (int i=0;i<N;i++) {
        if (fabs(z[i]-c) > EPSILON) {
            has_error = true;
            break;
        }
    }
    printf("%s\n",has_error?"Has errors":"No errors");
}

int main() {

    //创建数组xyz
    double* x = new double[N];
    double* y = new double[N];
    double* z = new double[N];

    //初始化数组
    for (int i=0;i<N;i++) {
        x[i] = a;
        y[i] = b;
        z[i] = 0;
    }

    //计算数组
    add(x,y,z,N);
    //检查结果
    check(z,N);

    //删除数组指针
    delete[] x;
    delete[] y;
    delete[] z;

    return 0 ;
}

CUDA程序和内存管理

先看代码,比起hello world,这是一个完整的CUDA程序,它包含了下图所示的结构:

image-20251102212858586.png

/**
 * cuda版:两数组相加
 *
*/
#include<iostream>

const double a=1.23;
const double b=3.45;
const double c=4.68;
const int N = 100000000;
const int M = sizeof(double)*N;
const double EPSILON = 1.0e-15;

// 函数1: 数组相加(数组默认是动态分配的)
__global__ void add(const double *x,const double *y,double *z) {
    const int n  = blockDim.x*blockIdx.x + threadIdx.x;
    z[n] = x[n] + y[n];
}

// 函数2: 检查结果
void check(const double *z,const int N) {
    bool has_error = false;
    for (int i=0;i<N;i++) {
        if (fabs(z[i]-c) > EPSILON) {
            has_error = true;
            break;
        }
    }
    printf("%s\n",has_error?"Has errors":"No errors");
}

int main() {

    //主机: 创建数组xyz
    double* x = new double[N];
    double* y = new double[N];
    double* z = new double[N];

    //初始化数组
    for (int i=0;i<N;i++) {
        x[i] = a;
        y[i] = b;
        z[i] = 0;
    }

    //主机: 创建设备数组xd,yd

    double* xd;
    double* yd;

    cudaMalloc((void **)&xd,M);
    cudaMalloc((void **)&yd,M);

    //主机: 复制主机数据到设备
    cudaMemcpy(xd,x,M,cudaMemcpyHostToDevice);
    cudaMemcpy(yd,y,M,cudaMemcpyHostToDevice);

    //主机: 指定设备运行参数
    //const dim3 block_size(128,1,1);
    const int block_size=128;
    const dim3 grid_size(N/block_size,1,1);

    //设备: 计算数组
    add<<<grid_size,block_size>>>(xd,yd,yd);

    //主机: 复制设备结果到主机
    cudaMemcpy(z,yd,M,cudaMemcpyDeviceToHost);

    //检查结果
    check(z,N);


    //删除数组指针
    cudaFree(xd);
    cudaFree(yd);

    delete[] x;
    delete[] y;
    delete[] z;

    return 0 ;
}

先把核函数的部分放在一边,代码里有以下几点可以注意:

  • 隐式设备初始化:了解过opengl或d3d程序的可以明显感觉,CUDA是没有显式初始化设备(即GPU)的函数。在第一次调用一个和设备管理及版本查询功能无关的运行时 API函数时,设备将自动初始化。
  • 内存管理一节已经大致了解了内存管理的相关函数,这里就是应用的地方,需要注意的是cudaMalloc和cudaMemcpy

    • cudaMalloc的第一个参数是指针的指针
    • cudaMemcpy的最后一个参数规定了数据传输的方式(从名称就可以看出)

CUDA程序和核函数

然后来看核函数,可分为三个方面:

  • 启动核函数
  • 编写核函数
  • 验证核函数

首先,CPU是一个控制者,运行核函数,要从CPU发起,第一,核函数是异步的,第二,如《CUDA的线程结构》所说,三个尖括号<<<grid,block>>>内是对设备代码执行的线程结构的配置:

  • 网格的布局(线程块)
  • 线程块的布局(线程)

其次,同步函数有两种

  • cudaDeviceSynchronize显式同步
  • 涉及到设备端不执行完,主机没办法进行,比如内存拷贝函数,隐式同步

然后看到核函数的写法,cuda中不仅有核函数,还有一些其他类型,通过不同的标识表明,如核函数是__global__

限定符执行调用备注
global设备端执行可以从主机调用也可以从计算能力3以上的设备调用必须有一个void的返回类型
device设备端执行设备端调用
host主机端执行主机调用可以省略
有些函数可以同时定义为 devicehost ,这种函数可以同时被设备和主机端的代码调用,主机端代码调用函数很正常,设备端调用函数与C语言一致,但是要声明成设备端代码,告诉nvcc编译成设备机器码,同时声明主机端设备端函数,那么就要告诉编译器,生成两份不同设备的机器码。

其他方面核函数基本和普通函数一致,但有一些限制:

  1. 除了统一寻址,只能访问设备内存
  2. 必须有void返回类型
  3. 不支持可变数量的参数
  4. 不支持静态变量
  5. 显示异步行为
  6. 不可作为成员,但可以包装
  7. 在动态并行的前提下,可以递归调用

另外设备函数可用__noinline__ __forceinline__来作为内联函数和非内联函数

CUDA小技巧,当我们进行调试的时候可以把核函数配置成单线程的:

下面是一些例子:

  • 函数的使用
  • 斐波那契数列,本程序意义不大,旨在展示__device__ __host__的共用
  • 斐波那契数列,本程序意义不大,旨在展示递归的使用
/**
 * 函数的使用
*/
#include<iostream>

const double a=1.23;
const double b=3.45;
const double c=4.68;
const int N = 100000000;
const int M = sizeof(double)*N;
const double EPSILON = 1.0e-15;

// 函数1: 数组相加(数组默认是动态分配的)
//TODO: 注意到三种写法效率是一样的
__device__ double addfunc1(double x,double y) {
    return x+y;
}
__device__ double addfunc2(double x,double y,double *z) {
    *z = x+y;
}
__device__ double addfunc3(double x,double y,double &z) {
    z= x+y;
}
__global__ void add(const double *x,const double *y,double *z) {
    const int n  = blockDim.x*blockIdx.x + threadIdx.x;
    if (n<N) {
        z[n] = x[n] + y[n];
        //z[n] = addfunc1(x[n],y[n]);
        //addfunc2(x[n],z[n],&z[n]);
        //addfunc3(x[n],y[n],z[n]);
    }
}
/**
 * 斐波那契数列,本程序意义不大,旨在展示__device__ __host__的共用
 *
*/
#include<iostream>


const int N = 10;
const int M = sizeof(int)*N;

// 函数1: 斐波那契数列
__device__ __host__ void fibN(int* A, const int N) {
    if (N<0) return;
    if (A[N]==-1) {
        if (N>1) {
            fibN(A,N-1);
            fibN(A,N-2);
            A[N] = A[N-1] + A[N-2];
        }
        if (N<=1) {
            A[N]=N;
        }
    }
}

__global__ void fibNN(int* A, const int N) {
    fibN(A, N);
}


// 函数2: 检查结果
//0、1、1、2、3、5、8、13、21、34
void check(const int *A,const int N) {
    for (int i=0;i<N;i++) {
        std::cout<<A[i]<<" ";
    }
    std::cout<<std::endl;
}

int main() {

    //主机: 创建数组xyz
    int* x = new int[N];


    //初始化数组
    for (int i=0;i<N;i++) {
        x[i] = -1;
    }

    //主机: 创建设备数组xd,yd

    int* xd;
    cudaMalloc((void **)&xd,M);

    //主机: 复制主机数据到设备
    cudaMemcpy(xd,x,M,cudaMemcpyHostToDevice);

    //主机: 指定设备运行参数
    const dim3 block_size(128,1,1);
    const dim3 grid_size(128,1,1);

    //设备: 计算数组
    fibNN<<<grid_size,block_size>>>(xd,N-1);


    //主机: 复制设备结果到主机
    cudaMemcpy(x,xd,M,cudaMemcpyDeviceToHost);

    //检查结果
    check(x,N);


    for (int i=0;i<N;i++) {
        x[i] = -1;
    }
    fibN(x,N-1);
    check(x,N);

    //删除数组指针
    cudaFree(xd);
    delete[] x;


    return 0 ;
}
/**
 * cuda版:两数组相加
 *
*/
#include<iostream>

const double a=1.23;
const double b=3.45;
const double c=4.68;
const int N = 100000000;
const int M = sizeof(double)*N;
const double EPSILON = 1.0e-15;

// 函数1: 数组相加(数组默认是动态分配的)
__global__ void add(const double *x,const double *y,double *z) {
    const int n  = blockDim.x*blockIdx.x + threadIdx.x;
    if (n<N) {
        z[n] = x[n] + y[n];
    }
}

// 函数2: 检查结果
void check(const double *z,const int N) {
    bool has_error = false;
    for (int i=0;i<N;i++) {
        if (fabs(z[i]-c) > EPSILON) {
            has_error = true;
            break;
        }
    }
    printf("%s\n",has_error?"Has errors":"No errors");
}

int main() {

    //主机: 创建数组xyz
    double* x = new double[N];
    double* y = new double[N];
    double* z = new double[N];

    //初始化数组
    for (int i=0;i<N;i++) {
        x[i] = a;
        y[i] = b;
        z[i] = 0;
    }

    //主机: 创建设备数组xd,yd

    double* xd;
    double* yd;

    cudaMalloc((void **)&xd,M);
    cudaMalloc((void **)&yd,M);

    //主机: 复制主机数据到设备
    cudaMemcpy(xd,x,M,cudaMemcpyHostToDevice);
    cudaMemcpy(yd,y,M,cudaMemcpyHostToDevice);

    //主机: 指定设备运行参数
    //const dim3 block_size(128,1,1);
    const int block_size=128;
    const dim3 grid_size(N/block_size,1,1);

    std::cout << "Size: " << N/block_size <<"=="<<block_size <<std::endl;

    //设备: 计算数组
    //TODO: 注意到这里使用了会超过N的执行参数
    add<<<800000,block_size>>>(xd,yd,yd);

    //主机: 复制设备结果到主机
    cudaMemcpy(z,yd,M,cudaMemcpyDeviceToHost);

    //检查结果
    check(z,N);


    //删除数组指针
    cudaFree(xd);
    cudaFree(yd);

    delete[] x;
    delete[] y;
    delete[] z;

    return 0 ;
}
/**
 * 斐波那契数列,本程序意义不大,旨在展示递归的使用
 *
*/
#include<iostream>


const int N = 10;
const int M = sizeof(int)*N;

// 函数1: 斐波那契数列
__global__ void fibN(int* A, const int N) {
    if (N<0) return;
    if (A[N]==-1) {
        if (N>1) {
            fibN<<<1,1>>>(A,N-1);
            fibN<<<1,1>>>(A,N-2);
            A[N] = A[N-1] + A[N-2];
        }
        if (N<=1) {
            A[N]=N;
        }
    }
}



// 函数2: 检查结果
//0、1、1、2、3、5、8、13、21、34
void check(const int *A,const int N) {
    for (int i=0;i<N;i++) {
        std::cout<<A[i]<<" ";
    }
    std::cout<<std::endl;
}

int main() {

    //主机: 创建数组xyz
    int* x = new int[N];


    //初始化数组
    for (int i=0;i<N;i++) {
        x[i] = -1;
    }

    //主机: 创建设备数组xd,yd

    int* xd;
    cudaMalloc((void **)&xd,M);

    //主机: 复制主机数据到设备
    cudaMemcpy(xd,x,M,cudaMemcpyHostToDevice);

    //主机: 指定设备运行参数
    const dim3 block_size(128,1,1);
    const dim3 grid_size(128,1,1);

    //设备: 计算数组
    fibN<<<grid_size,block_size>>>(xd,N-1);


    //主机: 复制设备结果到主机
    cudaMemcpy(x,xd,M,cudaMemcpyDeviceToHost);

    //检查结果
    check(x,N);


    //删除数组指针
    cudaFree(xd);
    delete[] x;


    return 0 ;
}

//0 1 1 0 -2 -2 -2 -4 -2 -6 结果错误 但是可见可以这么用

验证和错误检测

所有编程都需要对错误进行处理,早期的编码错误,编译器会帮我们搞定,内存错误也能观察出来,但是有些逻辑错误很难发现,甚至到了上线运行时才会被发现,而且有些厉害的bug复现会很难,不总出现,但是很致命,而且CUDA基本都是异步执行的,当错误出现的时候,不一定是哪一条指令触发的,这一点非常头疼;这时候我们就需要对错误进行防御性处理了

错误检测宏如下:

#pragma once
#include <cstdio>

#define CHECK(call)                                   \
do                                                    \
{                                                     \
    const cudaError_t error_code = call;              \
    if (error_code != cudaSuccess)                    \
    {                                                 \
        printf("CUDA Error:\n");                      \
        printf("    File:       %s\n", __FILE__);     \
        printf("    Line:       %d\n", __LINE__);     \
        printf("    Error code: %d\n", error_code);   \
        printf("    Error text: %s\n",                \
            cudaGetErrorString(error_code));          \
        exit(1);                                      \
    }                                                 \
} while (0)

  • cudaGetErrorString可把error_code转为相应的提示信息

用这个宏函数包装大部分CUDA运行时API函数。有一个例外是cudaEventQuery()函数,因为它很有可能返回cudaErrorNotReady,但又不代表程序出错了

#include "error.cuh"
#include <cmath>
#include <cstdio>

const double EPSILON = 1.0e-15;
const double a = 1.23;
const double b = 2.34;
const double c = 3.57;
void __global__ add(const double *x, const double *y, double *z, const int N);
void check(const double *z, const int N);

int main(void)
{
    const int N = 100000000;
    const int M = sizeof(double) * N;
    double *h_x = (double*) malloc(M);
    double *h_y = (double*) malloc(M);
    double *h_z = (double*) malloc(M);

    for (int n = 0; n < N; ++n)
    {
        h_x[n] = a;
        h_y[n] = b;
    }

    double *d_x, *d_y, *d_z;
    CHECK(cudaMalloc((void **)&d_x, M));
    CHECK(cudaMalloc((void **)&d_y, M));
    CHECK(cudaMalloc((void **)&d_z, M));
    //TODO:此处错误复制方向搞反了
    CHECK(cudaMemcpy(d_x, h_x, M, cudaMemcpyDeviceToHost));
    CHECK(cudaMemcpy(d_y, h_y, M, cudaMemcpyDeviceToHost));

    const int block_size = 128;
    const int grid_size = (N + block_size - 1) / block_size;
    add<<<grid_size, block_size>>>(d_x, d_y, d_z, N);

    CHECK(cudaMemcpy(h_z, d_z, M, cudaMemcpyDeviceToHost));
    check(h_z, N);

    free(h_x);
    free(h_y);
    free(h_z);
    CHECK(cudaFree(d_x));
    CHECK(cudaFree(d_y));
    CHECK(cudaFree(d_z));
    return 0;
}

void __global__ add(const double *x, const double *y, double *z, const int N)
{
    const int n = blockDim.x * blockIdx.x + threadIdx.x;
    if (n < N)
    {
        z[n] = x[n] + y[n];
    }
}

void check(const double *z, const int N)
{
    bool has_error = false;
    for (int n = 0; n < N; ++n)
    {
        if (fabs(z[n] - c) > EPSILON)
        {
            has_error = true;
        }
    }
    printf("%s\n", has_error ? "Has errors" : "No errors");
}
用上述方法不能捕捉调用核函数的相关错误,因为核函数不返回任何值(回顾一下,核函数必须用void修饰)。有一个方法可以捕捉调用核函数可能发生的错误
#include "error.cuh"
#include <cmath>
#include <cstdio>

const double EPSILON = 1.0e-15;
const double a = 1.23;
const double b = 2.34;
const double c = 3.57;
void __global__ add(const double *x, const double *y, double *z, const int N);
void check(const double *z, const int N);

int main(void)
{
    const int N = 100000000;
    const int M = sizeof(double) * N;
    double *h_x = (double*) malloc(M);
    double *h_y = (double*) malloc(M);
    double *h_z = (double*) malloc(M);

    for (int n = 0; n < N; ++n)
    {
        h_x[n] = a;
        h_y[n] = b;
    }

    double *d_x, *d_y, *d_z;
    CHECK(cudaMalloc((void **)&d_x, M));
    CHECK(cudaMalloc((void **)&d_y, M));
    CHECK(cudaMalloc((void **)&d_z, M));
    CHECK(cudaMemcpy(d_x, h_x, M, cudaMemcpyHostToDevice));
    CHECK(cudaMemcpy(d_y, h_y, M, cudaMemcpyHostToDevice));


    //TODO:错误在线程块的大小太大
    const int block_size = 1280;
    const int grid_size = (N + block_size - 1) / block_size;
    add<<<grid_size, block_size>>>(d_x, d_y, d_z, N);
    CHECK(cudaGetLastError());
    CHECK(cudaDeviceSynchronize());

    CHECK(cudaMemcpy(h_z, d_z, M, cudaMemcpyDeviceToHost));
    check(h_z, N);

    dim3 a;
    free(h_x);
    free(h_y);
    free(h_z);
    CHECK(cudaFree(d_x));
    CHECK(cudaFree(d_y));
    CHECK(cudaFree(d_z));
    return 0;
}

void __global__ add(const double *x, const double *y, double *z, const int N)
{
    const int n = blockDim.x * blockIdx.x + threadIdx.x;
    if (n < N)
    {
        z[n] = x[n] + y[n];
    }
}

void check(const double *z, const int N)
{
    bool has_error = false;
    for (int n = 0; n < N; ++n)
    {
        if (fabs(z[n] - c) > EPSILON)
        {
            has_error = true;
        }
    }
    printf("%s\n", has_error ? "Has errors" : "No errors");
}

CHECK(cudaGetLastError());
CHECK(cudaDeviceSynchronize());

其中,第一条语句的作用是捕捉第二个语句之前的最后一个错误,第二条语句的作用是同步主机与设备。之所以要同步主机与设备,是因为核函数的调用是异步的,即主机发出调用核函数的命令后会立即执行后面的语句,不会等待核函数执行完毕。

万物一齐,而百兽率舞。

本文简略介绍CUDA编程的相关背景

异构计算

首先,异构计算就是CPU+X,这个X可以是GPU,也可以是NPU,FPGA,DSP等等,在CPU上运行的代码称为Host主机代码,X上运行的代码称为Device设备代码

其次,异构计算和同构计算的一个不同点在于,异构的处理要比同构复杂,因为涉及同构下自动处理的计算,控制,传输都要人工进行干预

GPU

CPU vs GPU

先看cpu的架构和gpu架构的不同:

1.png

先看几个共同的部分:

  • 控制器Control
  • 缓存器Cache
  • 随机存取存储器DRAM
  • 总线PCIe Bus
  • 算术逻辑单元ALU

不需要过多的计组原理知识,直观来看,GPU弱化了控制器的能力,大大堆砌了ALU的数量,且为一组ALU(SM,简单理解为线程束)分配了单独的Control和cache,所以直觉上,对了逻辑简单,数据量大的任务,GPU更高效

低并行逻辑复杂的程序适合用CPU
高并行逻辑简单的大数据计算适合GPU

CPU和GPU线程的区别:

  1. CPU线程是重量级实体,操作系统交替执行线程,线程上下文切换花销很大
  2. GPU线程是轻量级的,GPU应用一般包含成千上万的线程,多数在排队状态,线程之间切换基本没有开销。
  3. CPU的核被设计用来尽可能减少一个或两个线程运行时间的延迟,而GPU核则是大量线程,最大幅度提高吞吐量

GPU简介

对于Nvida,GPU分为4个系列:

  • Tegra
  • Geforce
  • Quadro
  • Tesla

针对不同的应用场景,比如Tegra用于嵌入式,Geforce是平时打游戏用到,Tesla主要用于计算,但目前看Geforce和Tesla貌似也没那么分得很清楚,只是某些系列上有容错的硬件部分

GPU也存在版本号,即其计算能力,但这和性能不是成正比的,只代表架构的更新,对于CUDA来说,新版的cuda只支持最近几个版本的GPU,参考:https://en.wikipedia.org/wiki/CUDA#GPUs_supported

真正评价GPU能力的分为两种,一种是计算性能指标,一种是容量指标:

  • 计算性能主要看FLOPS,一般单位在T级别(TFLOPS),称每秒浮点数运算次数,分单精度双精度,根据不同的架构,双精度一般是单精度的1/N
  • 容量指标上有:显存,带宽,核心数量等,特别的如果一个CUDA程序爆显存了,则不使用统一内存的情况下是无法运行程序的

CUDA

CUDA平台不是单单指软件或者硬件,而是建立在Nvidia GPU上的一整套平台,并扩展出多语言支持

对于CUDA API来说,分为driver和运行时两种API,二者互斥不能同时使用,性能没啥差距

  • driver API 低级,复杂,灵活,一般用于配合其他语言调用CUDA
  • 运行时简单,基于driver API,稍微高级一点

CUDA编程环境架构如下:

5.png

nvcc

nvcc之于cuda相当于普通c++之于g++。因为cuda程序中host和device的代码是混在一起的,nvcc会分离这两个部分分别处理:

6.png

其中,PTX相当于CUDA的汇编语言

nvdia-smi

可以通过nvidia-smi(Nvidia's system management interface)程序检查与设置设备。它包含在CUDA开发工具套装内。

CUDA的程序结构

一般CUDA程序分成下面这些步骤:

  1. 分配GPU内存
  2. 拷贝内存到设备
  3. 调用CUDA内核函数来执行计算
  4. 把计算完成数据拷贝回主机端
  5. 内存销毁

参考

  1. https://face2ai.com/program-blog/#GPU%E7%BC%96%E7%A8%8B%EF%BC%88CUDA%EF%BC%89

万物化作,萌区有状,盛衰之杀,变化之流也。

本文以helloworld为例,简述CUDA的线程结构

基本概念

即使是HelloWorld,Cuda程序也引入了一些不同的概念,确保了解《简略CUDA和GPU相关概念》中的内容在进行下去

从HelloWorld开始

下面是一个最简单的HelloWorld程序:

/**
 运行命令: cl ./hello_host_only.cpp
**/

#include <cstdio>
using namespace std;

int main() {
    printf("Host: Hello wolrd!");
    return 0;
}

在windows可在x64 Native Tools Command Prompt for VS 2022cl来编译它

只有Host代码的CUDA版本HelloWolrd

/**
 运行命令: nvcc -arch=sm_62 hello_host_only.cu
 sm_62的来源:
 https://docs.nvidia.com/cuda/archive/12.9.0/cuda-compiler-driver-nvcc/index.html
 https://en.wikipedia.org/wiki/CUDA#GPUs_supported
*/

#include <cstdio>
using namespace std;

int main() {
    printf("Host: Hello wolrd!");
    return 0;
}

代码完全没有变化,只有文件后缀名发生了变化

但这里性质发生了变化,代码分为了Host代码部分和Device代码部分由nvcc分别处理

nvcc指定运算架构

注意到注释部分,nvcc可以指定虚拟和真实架构的编号,可以通过链接查看本地设备可以使用的标识

如果需要匹配多个架构,可以使用多个-gencode -arch 虚拟架构 -code 真实架构来编译产生不同的输出,最后真实架构不能落后与虚拟架构

使用CUDA线程

从这里正式开始添加Device代码:

/**
 运行命令: nvcc  XX.cu
*/
#include<cuda.h>
#include<cuda_runtime.h>
#include<cstdio>

__global__ void hello_from_gpu()
{
    printf("Device: Hello World.");
}
int main()
{
    const int grid_size=1;
    const int block_size=1;
    hello_from_gpu<<<grid_size,block_size>>>();

    cudaDeviceSynchronize();
    return 0;
}

注意四点:

  • include包含了两个cuda相关头文件,这不是必须要添加的,nvcc会自动引入,其中包含了cstdlib的内容,但没有包含cstdio的内容,最后device程序是不认iostream的
  • __global___就是device代码,即核函数的标志
  • cudaDeviceSynchronize的功能是同步host和device,也就是说gpu上的代码和cpu上的代码是异步运行的
  • <<<grid_size,block_size>>>指定了用多少线程运行该核函数

使用多线程的CUDA程序

再进一步改进程序

/**
 运行命令: nvcc  XX.cu
*/
#include<cuda.h>
#include<cuda_runtime.h>
#include<cstdio>

__global__ void hello_from_gpu()
{
    const int bid = blockIdx.x;
    const int tid = threadIdx.x;
    printf("Device: Hello World from block %d and thread %d\n",bid,tid);
}
int main()
{
    const int grid_size=2;
    const int block_size=4;
    hello_from_gpu<<<grid_size,block_size>>>();

    cudaDeviceSynchronize();
    return 0;
}

CUDA将线程组织为网格,网格被划分为线程块,即网格由线程块组成,线程块由线程组成,三者之间并不是简单的包含关系(即一般论的包含关系,网格有N个线程块,每个线程块又有M个线程),但是目前示例中确实只是单纯的包含关系,可简单认为有网格两个线程块,每个线程块四个线程共八个线程

下图展示了多维线程的布局,这也是网格的真正分布:

scanner_20251011_221652.jpg

可以认为网格是一个三维立方体,被划分出的小格就是线程块,而小格又可以按三维继续分割,分割出的就是线程。

像上面简单的以2和4指定的线程块个数和线程个数,三维上可以分别认为是(2,1,1)和(4,1,1),只按x轴延伸的2个大方块,每个大方块又是由按x轴延展的4个小方块组成

多维线程的CUDA程序

既然有了三维建构,每个线程块和线程自然有其三维坐标,按照简单的数学逻辑也可以转换为1维坐标,这里直接把这些概念在代码中展示:

/**
 运行命令: nvcc  XX.cu
*/
#include<cuda.h>
#include<cuda_runtime.h>
#include<cstdio>

__global__ void hello_from_gpu()
{
    const int bid = blockIdx.x;
    const int tid_x = threadIdx.x;
    const int tid_y = threadIdx.y;
    const int tid_z = threadIdx.z;
    //1D坐标:
    const int bid_1D = blockIdx.x+blockIdx.y*gridDim.x+blockIdx.z*gridDim.x*gridDim.y;
    const int tid_1D = threadIdx.x + blockDim.x*threadIdx.y + blockDim.x*blockDim.y*threadIdx.z;
    printf("Device(%d)(%d): Hello World from block %d and thread (%d-%d-%d)\n",bid_1D,tid_1D,bid,tid_x,tid_y,tid_z);
}
int main()
{
    const int grid_size=2;
    const dim3 block_size(3,3,3);
    hello_from_gpu<<<grid_size,block_size>>>();

    cudaDeviceSynchronize();
    return 0;
}

注意到:

  • dim3类似于vec3,表3维向量
  • gridDim和blockDim是dim3固有变量,内容是当初规定的各个轴的上限
  • blockIdx和threadIdx也是dim3固有变量,包含了当前线程的坐标

另外,线程块在3个轴上的个数和一个线程块的线程数是有限制的:

  • 线程块在xyz三轴上分别最多:2^31-1,65535,65535
  • 一个线程块线程总数最多1024,在三维上最多:1024,1024,64

下面是超过了限制的程序:

/**
 运行命令: nvcc  XX.cu
*/
#include<cuda.h>
#include<cuda_runtime.h>
#include<cstdio>

__global__ void hello_from_gpu()
{
    const int bid = blockIdx.x;
    const int tid_x = threadIdx.x;
    const int tid_y = threadIdx.y;
    const int tid_z = threadIdx.z;
    //1D坐标:
    const int bid_1D = blockIdx.x+blockIdx.y*gridDim.x+blockIdx.z*gridDim.x*gridDim.y;
    const int tid_1D = threadIdx.x + blockDim.x*threadIdx.y + blockDim.x*blockDim.y*threadIdx.z;
    printf("Device(%d)(%d): Hello World from block %d and thread (%d-%d-%d)\n",bid_1D,tid_1D,bid,tid_x,tid_y,tid_z);
}
int main()
{
    const dim3 grid_size(1,65535,1); //大于规定的2^31-1,65535,65535范围虽然能通过编译,但没有输出
    const dim3 block_size(1024,1,1); //大于1024虽然能通过编译,但没有输出
    hello_from_gpu<<<grid_size,block_size>>>();

    cudaDeviceSynchronize();
    return 0;
}

线程束

还有一个没有在代码中体现出的概念:线程束,简单理解为同一线程块中相邻的N个线程,一般按架构为32个一组

image-20251011233201443.png