Daniel Warfield使用Midjourney的《Thread Master》。除非另有说明，否则作者的所有图像。文章最初在Intuitively and Exhaustively Explained上提供。

在本文中，我们将使用CUDA在GPU上训练AI模型，基本上是从头开始实现AI，假设几乎没有先验知识。

首先，我们将探索现代计算机的一些核心组件，然后我们将深入研究GPU，以描述它是什么，它是如何工作的，以及为什么它对人工智能有用。然后，我们将完成对CUDA的介绍。我们将描述什么是CUDA，并解释它如何允许我们编程同时利用CPU和GPU的应用程序。一旦我们了解了CUDA编程的工作原理，我们将使用CUDA在分类任务上构建、训练和测试神经网络。

这对谁有用？任何想深入了解人工智能的人。

这个帖子有多高级？鉴于高级主题，对于那些有一些机器学习经验的人来说，这篇文章可能更平易近人。不过，如果你没有机器学习经验，你肯定会通过阅读这篇文章学到很多东西。只需逐节取，然后用谷歌搜索很多。

先决条件：基本软件开发技能。接触C++可能会有所帮助，但不是必需的。熟悉衍生物的一般概念也可能是有用的。

归因：本文的灵感主要来自两个教程。这个关于CUDA的精彩YouTube系列，以及这个关于在CUDA中实现神经网络的教程。

为什么CUDA值得学习

在我们深入研究这篇大文章之前，我想解决一个可能出现在许多数据科学家脑海中的问题。“当我只能使用PyTorch时，我为什么要打扰CUDA呢？”

PyTorch是一个机器学习库，允许在GPU上创建、训练和运行高级AI模型。在很大程度上，过去几年人工智能的许多快速进步都归功于PyTorch的易用性和强大功能。

但是，在人工智能进步的前沿，PyTorch并不总是足够的。例如，Flash Attention通过重新设计PyTorch在CUDA中实现的注意力，将注意力速度提高了10倍。这篇文章不是关于Flash Attention（我计划很快报道），而是触及了一个更大的趋势；随着人工智能作为一门学科的成熟，效率正变得越来越重要。

研究是明确的——更大的模型更好。人工智能模型的训练和服务效率越高，它们就越大，就越有用。虽然PyTorch在推动最先进的技术水平方面至关重要，但低级CUDA仍然是高效实施尖端人工智能的重要工具。

此外，PyTorch使用CUDA与GPU接口，因此即使您使用PyTorch，您可能也在引擎盖下使用CUDA。您可以在CUDA中创建自定义PyTorch功能，因此如果您已经使用PyTorch，CUDA可以成为学习的有用技能。

现代计算机的组成

现代计算机的一个例子。本示例中使用的3D模型由来源Daniel Cardona命名为“Dream Computer Setup”。

与这个一般模型有偏差（即Groq）但99%的计算机由以下元素组成：

典型计算机的主要组件。

为了本文（使用CUDA在GPU上实现AI）的目的，我们主要感兴趣的是了解主板、CPU、RAM和GPU。

主板是计算机的支柱。它本质上是一个大电路板，允许计算机的组件相互通信。

主板，将所有计算机组件连接在一起。

主板中间是几乎每台现代计算机的跳动心脏：CPU或“中央处理单元”。CPU负责执行运行程序所需的计算，因此是计算机最重要的组件之一。

CPU，执行运行计算机所需的主要计算。

因为CPU产生热量，所以CPU通常被散热器覆盖，这就是为什么我们无法在这个渲染中看到它。通常，CPU看起来像一个金属正方形，底部有一堆针脚。

英特尔至强CPU源的图像，由英特尔的使用条款授予。

非常接近CPU的是RAM，它代表“随机访问内存”。RAM是CPU的工作内存，旨在允许快速访问相关信息。

RAM，存储与CPU相关的信息。

许多计算机的另一个流行组件，也是本文的主题，是GPU，或“图形处理单元”。我们将在以下部分中确切了解GPU的作用，但目前只需知道GPU旨在帮助CPU进行某些类型的计算就足够了。

GPU，旨在加速某些计算。

玩家喜欢他们的GPU，所以这个是正面和中央的；使用带状电缆将卡与主板连接。

在本例中，带状电缆将GPU连接到主板。

像许多消费硬件空间一样，这纯粹是美学。在企业设置中，GPU将尽可能接近其他组件，因此它可能会直接安装到PCI-E端口。

PCI-E或“Peripheral Component Interconnect Express”是许多主板上的一组端口，允许CPU与外部设备通信。PCI-E最常见的用途是连接GPU，但PCI-E是一个灵活的接口；允许将存储设备、专用卡和其他设备连接到计算机。

PCI-E总线

在许多方面，PCI-E总线上的组件是计算机核心的“插件”。CPU和RAM对计算机的运行至关重要，而像GPU这样的设备就像CPU可以激活的工具，可以执行某些事情。即使在最低的编程水平上，这种概念化仍然相关。因此，CPU和RAM通常被称为“主机”，而GPU被称为“设备”。

CPU和RAM是计算机的基本组件，有时统称为主机。GPU可以被CPU调用来协助一些任务，有时也被称为设备。

在结构上，如果我们移除GPU的外壳并检查其内部工作，我们发现GPU与主机（CPU和RAM）有些相似。

如果我们从GPU上拆下风扇组件和散热片，我们可以看到GPU的底层结构。ramyissa的名为“RTX 3060 TI Zotac Twin Edge OC（GPU模型）”的模型。来源。

在这一点上，澄清我们的定义是有意义的。虽然整个显卡通常被称为GPU，但实际上GPU一词是指显卡中的处理单元。

虽然显卡的整个组件通常被称为GPU，但实际上显卡的主要处理单元是GPU。

除了GPU，显卡还有其他主要组件，如vRAM，它本质上相当于CPU的RAM的显卡。

显卡内的VRAM。

因此，我们得出本节的主要结论：设备（显卡）有一个处理单元和内存。主机（CPU和RAM）有一个处理单元和内存。设备和主机都有独立工作的必要资源。构建同时使用设备和主机的低级应用程序需要对两个协同工作的实体进行编程。

可能会出现一个自然的问题：“如果设备和主机如此相似，只是一个计算模块和一些内存，为什么我们需要它们”？

为什么GPU在人工智能中很重要

要了解GPU存在的原因，首先更彻底地了解CPU是有帮助的。我将通过这一节进行闪电战，但如果您有兴趣深入挖掘，我会在另一篇文章中更彻底地介绍这些概念：

Groq和人工智能的硬件——直观而详尽的解释

分析用于运行人工智能的主要计算机硬件，以及一个新的重击手。

towardsdatascience.com

CPU由“核心”组成。核心是设计用于快速完成基本数学运算的电路。虽然简单，但核心可以做的数学是基本的，结合许多这些简单的计算可以产生几乎任何可以想象的数学运算。

现代CPU有几个这样的内核，因此CPU可以同时处理多件事物。CPU还具有不同类型的内存（称为缓存），旨在加快CPU的速度。

现代CPU主要组件的框图。拥有多个内核允许现代计算机同时处理几件事。

CPU的整个想法是运行一个程序，该程序可以尽快被概念化为操作列表。虽然CPU可以并行做事情（称为多线程），但CPU的主要重点是尽可能快地进行背靠背计算。因此，CPU的每个核心都旨在尽快完成每个操作。

这就是CPU设计人员试图实现的目标。这个想法是通过让CPU尽快在该程序中运行任务来使程序更快。来自我关于人工智能硬件的文章。

这种方法对许多应用程序至关重要。一些计算只需要一个接一个地进行，更快完成这些计算的唯一方法是尽可能快地完成每一步。然而，一些应用程序可以并行化。

CPU通过快速完成单个任务来使事情变得快速。GPU通过并行执行任务来快速完成任务。来自我关于人工智能硬件的文章。

GPU的想法不是优化单个计算，而是优化批量运行许多计算。一般来说，在给定的计算中，GPU比CPU慢得多，因此只有当有足够的计算可以并行运行时，GPU才会真正发光。然而，当应用于正确的问题时，GPU可以比CPU快100倍甚至1000倍。

CPU可以非常快速地完成单个计算，而GPU可以并行完成更多的计算，但速度较慢。归根结底，GPU可以完成更多的计算，但你必须等待它们通过。来自我关于人工智能硬件的文章

在专注于并行计算时，GPU与CPU有两个核心偏差：

1. GPU的内核数量比CPU多，但这些内核的能力较差。这意味着GPU可以一次进行更多计算，但这些计算通常较慢。
2. GPU的设置允许SIMD，即“单指令多数据”。基本上，单个控制电路控制GPU中的多个内核。

CPU和GPU的概念图。两者之间有很多相似之处，但对单个操作速度与并行性的关注在性能方面有重大影响。来自我关于AI硬件的文章

最初，GPU是为视频游戏渲染图形而创建的，但此后人气激增，因此NVIDIA的估值激增。

NVIDIA股价，来源。

这一爆炸性成功的原因主要是因为一个偶然的巧合。事实证明，除了图形，你还可以用GPU做许多其他有用的事情。例如，人工智能模型是GPU的完美用例，因为人工智能基本上是通过一堆简单且或多或少独立的计算完成的。

人工智能模型的简化概念图。即使在这个简单的模型中，计算输出所需的31个顺序操作也可以分为9个并行步骤，理论上使模型的运行速度快2.4倍。在GPU上运行时，真正的人工智能模型可以看到100倍、1000倍甚至更高的速度。来自我关于Groq的文章。

好的，所以我们了解计算机的核心组件，也理解为什么GPU在人工智能中有用。让我们开始了解如何让GPU做事情。

CUDA简介

完整的代码可以在这里找到。

CUDA，或“计算统一设备架构”，是NVIDIA的并行计算平台。CUDA本质上是一套用于构建在CPU上运行的应用程序的工具，可以与GPU接口以进行并行数学。

运行CUDA最受欢迎的语言可能是C++，所以这就是我们将要使用的。不过别担心，我是一名主要使用Python的数据科学家，所以本文中的C++应该相当平易近人（希望不是更资深的C++开发人员的判断来源）。

我将使用Google Colab作为我的开发环境。有一个方便的Jupyter扩展，可以编译和运行CUDA代码，就像它是一个普通的代码块一样。这允许Jupyter连接到nvcc,NVIDIA的CUDA编译器。幸运的是，GoogleColab已经设置了nvcc，所以我们只需要安装和加载扩展。

！pip安装nvcc4jupyter

%load_ext nvcc4jupyter

使用nvcc4jupyter，我们可以像普通python块一样运行CUDA代码块，方法是在单元格顶部添加%%cuda：

%%cuda

//we'll cover what all this stuff means soon
#include 

__global__ void hello(){
    printf("Hello from block: %u, thread: %u\n", blockIdx.x, threadIdx.x);
}

__host__ int main(){
    hello<<<2, 2>>>();
    cudaDeviceSynchronize();
}

使用%%cuda运行单元格块的结果

文章读到一半时，我们将使用一种略微不同的方法来启动CUDA，这使我们能够处理更复杂的项目，但现在我们可以使用这种简单的方法来探索CUDA的工作原理。

CUDA内核和启动配置

CUDA的第一个基本想法是内核。基本上，CUDA内核是一个跨GPU并行运行多次的函数。在CUDA中，我们在函数定义之前用__global__关键字定义内核。我们还可以使用__host__关键字定义在CPU上运行的代码。

%%cuda
#include 

//this runs on the GPU because it's __global__
__global__ void hello(){
    printf("Hello from block: %u, thread: %u\n", blockIdx.x, threadIdx.x);
}

//this runs on the CPU because it's __host__
__host__ int main(){
    hello<<<2, 2>>>();
    cudaDeviceSynchronize();
}

因为这是C++，首先执行的是main函数，它恰好在CPU（又名主机）上。然后，在CPU上，被调用的第一行代码是hello<<<2, 2>>>();。这被称为启动配置，并在GPU上启动CUDA内核hello。

在CUDA中，并行作业被组织成“线程”，其中线程是可以并行工作的程序的一部分。您可以在GPU上拼凑多个线程，它们都可以一起工作来完成一些任务。

多线程允许您将大型任务划分为许多较小的并行任务。

这些线程存在于所谓的“线程块”中。在现代GPU上，每个线程块可以有1024个线程。同一线程块中的线程在vRAM（相当于RAM的GPU）中共享相同的地址空间，因此它们能够在相同的数据上协同工作。

线程块允许线程在同一共享内存上协同工作

虽然每个块只能有1024个线程，但您可以在GPU上排队大量线程块。

大量线程块可以排队并准备在GPU上执行

大多数GPU都有资源一次执行多个线程块。一个GPU可能一次只能执行两个线程块，另一个可能能够执行其中八个线程块。这意味着相同的CUDA代码可以利用大型和小型GPU。

一个小GPU和大GPU处理同一任务的例子。

我们定义在启动配置中将使用多少线程和线程块。在三个雪佛龙<<< >>>，我们指定线程块的数量作为第一个参数，将每个块的线程数量指定为第二个参数。

%%cuda
#include 

__global__ void hello(){
    printf("Hello from block: %u, thread: %u\n", blockIdx.x, threadIdx.x);
}

__host__ int main(){
    //Launching out kernel across three thread blocks,
    //each consisting of two threads
    hello<<<3, 2>>>();
    cudaDeviceSynchronize();
}

注意：由于这些是并行执行的，因此没有先得出结论的。

由于内核在GPU上运行，为了将打印结果重新打印到CPU上，以便在jupyter中显示，我正在使用cudaDeviceSynchronize等待GPU上的所有线程完成，这允许打印输出从设备传播到主机。

当您运行内核时，CUDA会自动创建一些方便的变量。两个非常重要的是blockIdx和threadIdx。您可以使用这些变量来了解当前正在运行哪个块和哪个线程。

线程块可以最多以三维方式组织线程，块本身可以最多以三维方式组织线程，这意味着blockIdx和threadIdx具有x、y和z属性

%%cuda
#include 

__global__ void hello() {
    printf("Hello from block: (%u, %u, %u), thread: (%u, %u, %u)\n", 
           blockIdx.x, blockIdx.y, blockIdx.z, 
           threadIdx.x, threadIdx.y, threadIdx.z);
}

int main() {
    // Define the dimensions of the grid and blocks
    dim3 gridDim(2, 2, 2);   // 2x2x2 grid of blocks
    dim3 blockDim(2, 2, 2);  // 2x2x2 grid of threads per block

    // Launch the kernel
    hello<<>>();

    // Wait for GPU to finish before accessing on host
    cudaDeviceSynchronize();

    return 0;
}

换句话说，CUDA允许您构建多达三维的计算网格，并通过blockIdx和threadIdx通知线程它在该网格中的位置。我们将在未来的一节中进一步探索这个想法。

主机和设备之间的通信

要将数据发送到GPU，我们可以首先使用cudaMalloc在vRAM上保留一些位置，我们可以使用cudaMemcpy在RAM和vRAM之间复制数据，当我们在GPU上完成时，我们可以使用cudaFree让GPU知道我们不再需要vRAM中的数据。在数据发送到vRAM后，我们可以在该数据上运行内核，然后将结果复制回RAM，有效地使用GPU进行数学运算。

让我们探索这个代码，它向GPU发送两个值，使用GPU将这些数字相加，然后将结果重新放入RAM并打印结果。

%%cuda
#include 
#include 

using namespace std;

// Defining the kernel
__global__ void addIntsCUDA(int *a, int *b) {
    a[0] += b[0];
}

// Running main on host, which triggers the kernel
int main() {
    // Host values
    int a = 1, b = 2;

    //printing expression
    cout << a << " + " << b <<" = ";

    // Device pointers (GPU)
    int *d_a, *d_b;

    // Allocating memory on the device (GPU)
    cudaMalloc(&d_a, sizeof(int));
    cudaMalloc(&d_b, sizeof(int));

    // Copying values from the host (CPU RAM) to the device (GPU)
    cudaMemcpy(d_a, &a, sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, &b, sizeof(int), cudaMemcpyHostToDevice);

    // Calling the kernel to add the two values at the two pointer locations.
    addIntsCUDA<<<1, 1>>>(d_a, d_b);

    // The addition function overwrites the a pointer with the sum. Thus
    // this copies the result.
    cudaMemcpy(&a, d_a, sizeof(int), cudaMemcpyDeviceToHost);

    //printing result
    cout << a << endl;

    //freeing memory.
    cudaFree(d_a);
    cudaFree(d_b);

    return 0;
}

上述代码的结果。a和b的值被复制到vRAM上，GPU计算总和并将结果写入vRAM，然后将总和复制回RAM。

对于那些对C++不适应的人来说，像*和&之类的东西可能会让你大吃一现。这些都与指针有关。基本上，你可以把RAM和vRAM都视为一个大的值阵列。当代码int a = 1, b = 2;在CPU上触发时，RAM中的两个点被分配来存储这两个值。

当a和b被创建时，RAM中为每个变量保留了某些位置。

这些值中的每一个都有一些地址。当我们调用&a或&b我们得到的是这些价值观的地址，而不是价值观本身。

通过使用&获取值的地址

如果您想将地址作为变量进行跟踪，您可以使用*创建称为指针的东西。指针是内存中指向内存中一些其他数据的数据。例如，如果我们运行代码int *pointer_to_a = &a，我们将在内存中创建一个点，该点在内存中保存对应于a点。

创建指向一个指针的结果，用于演示目的。

这种指向内存中数据的想法在低级计算中至关重要，因此在CUDA中至关重要。因为我们处理的是两个独立的系统（设备和主机），我们通常会创建指向两个内存空间之间的指针。

在这个例子中，在我们创建整数后，我们创建了两个指针

// Device pointers (GPU)
int *d_a, *d_b;

这些是CPU RAM上的指针，最初没有值，但它们是带有d_的名称，因为这些指针将存储GPU vRAM上的数据位置。

此时RAM的概念图。RAM上定义了两个空指针。

在此之后，我们调用cudaMalloc，它旨在在GPU的vRAM上分配一些内存点。

// 在设备上分配内存（GPU）
cudaMalloc(&d_a, sizeof(int));
cudaMalloc(&d_b, sizeof(int));

运行cudaMalloc对vRAM的影响

然而，为了在内存中实际使用这些位置，主机（CPU/RAM）需要了解它们。这就是为什么我们将&d_a和&d_b传递到cudaMalloc，我们告诉cudaMalloc在RAM上将新的vRAM分配存储在哪里。

运行cudaMalloc对vRAM和RAM的影响

因此，通过运行cudaMalloc，我们现在在vRAM上分配了点，以及RAM上的指针，这些指针告诉我们这些点在vRAM上的位置。

我们可以使用以下代码将a和b的值复制到vRAM

// 将值从主机（CPU RAM）复制到设备（GPU）
cudaMemcpy(d_a, &a, sizeof(int), cudaMemcpyHostToDevice);
cudaMemcpy（d_b，&b，sizeof（int），cudaMemcpyHostToDevice）；

此代码告诉cudaMemcpy我们希望数据存储在哪里（d_a和d_b）和我们想要复制的数据的位置（&a和&b），以及我们是否从主机复制到设备，反之亦然。这导致a和b的值被复制到设备中。

然后，当我们启动内核时，我们可以传递内核d_a和d_b，以便内核知道vRAM上的哪些值需要修改。

// Calling the kernel to add the two values at the two pointer locations.
addIntsCUDA<<<1, 1>>>(d_a, d_b);

addsIntsCUDA在GPU上运行，GPU获取指针，并将b值添加到a

// Defining the kernel
__global__ void addIntsCUDA(int *a, int *b) {
    a[0] += b[0];
}

运行内核的结果

然后，d_a指向的值被重新添加到RAM上，然后打印其值，从而打印最终结果。

// 加法函数用总和覆盖指针。因此
// 这复制了结果。
cudaMemcpy（&a，d_a，sizeof（int），cudaMemcpyDeviceToHost）；

//打印结果
cout << a << endl;

此时的记忆状态。

GPU不知道CPU即将结束其程序，因此GPU上的值将一直存在，直到应用程序完成。对于这个小问题，这并不重要，但如果你必须添加1和2十亿次，你可能不希望十亿1和2占用vRAM的空间。因此，我们调用cudaFree来释放内存。

//freeing memory.
cudaFree(d_a);
cudaFree(d_b);

释放内存中的两个点后的vRAM。

显然，这不是将两个数字加在一起的最有效方法，但我希望它能让你了解数据如何在RAM和vRAM之间来回传递。在我们从头开始构建人工智能模型之前，我认为回顾一个稍微高级的示例可能会有所帮助，这样我们就可以真正感受到对CUDA的某种程度的掌握。

在GPU上并行化CPU程序

在前面的几节中，我们讨论了如何启动内核，以及如何在主机和设备之间传输数据。让我们用这些概念来并行化为CPU编写的程序。

下面是一个蛮力实现，它试图为一组点中的每个3D点找到最接近的其他点。最终结果应该是一个列表，其中列表中i点中的值应该是i点最接近的点的索引。因此，如果有三个点，点1最接近第3点，点2最接近第3点，第3点最接近第2点，输出将看起来像[3, 3, 2]。以下是CPU上的实现：

%%cuda
#include 
#include 
#include 
#include 
#include 

using namespace std;

//brute force approach to finding which point
void findClosestCPU(float3* points, int* indices, int count) {
    // Base case, if there's 1 point don't do anything
    if(count <=1) return;
    // Loop through every point
    for (int curPoint = 0; curPoint < count; curPoint++) {
        // set as close to the largest float possible
        float distToClosest = 3.4028238f ;
        // See how far it is from every other point
        for (int i = 0; i < count; i++) {
            // Don't check distance to itself
            if(i == curPoint) continue;
            float dist_sqr = (points[curPoint].x - points[i].x) *
                (points[curPoint].x - points[i].x) +
                (points[curPoint].y - points[i].y) *
                (points[curPoint].y - points[i].y) +
                (points[curPoint].z - points[i].z) *
                (points[curPoint].z - points[i].z);
            if(dist_sqr < distToClosest) {
                distToClosest = dist_sqr;
                indices[curPoint] = i;
            }
        }
    }
}

int main(){

    //defining parameters
    const int count = 10000;
    int* indexOfClosest = new int[count];
    float3* points = new float3[count];

    //defining random points
    for (int i = 0; i < count; i++){
        points[i].x = (float)(((rand()%10000))-5000);
        points[i].y = (float)(((rand()%10000))-5000);
        points[i].z = (float)(((rand()%10000))-5000);
    }

    long fastest = 1000000000;

    cout << "running brute force nearest neighbor on the CPU..."<

在CPU上运行几次迭代

为了在GPU上并行化此代码，我们需要在vRAM上获取点，然后在GPU上作为内核启动findClosest。

这是GPU上的并行化版本，让我们来浏览一下：

%%cuda
#include 
#include 
#include 
#include 
#include 

using namespace std;

// Brute force implementation, parallelized on the GPU
__global__ void findClosestGPU(float3* points, int* indices, int count) {
    if (count <= 1) return;
    int idx = threadIdx.x + blockIdx.x * blockDim.x;
    if (idx < count) {
        float3 thisPoint = points[idx];
        float smallestSoFar = 3.40282e38f;

        for (int i = 0; i < count; i++) {
            if (i == idx) continue;

            float dist_sqr = (thisPoint.x - points[i].x) *
                             (thisPoint.x - points[i].x) +
                             (thisPoint.y - points[i].y) *
                             (thisPoint.y - points[i].y) +
                             (thisPoint.z - points[i].z) *
                             (thisPoint.z - points[i].z);
            if (dist_sqr < smallestSoFar) {
                smallestSoFar = dist_sqr;
                indices[idx] = i;
            }
        }
    }
}

int main() {
    // Defining parameters
    const int count = 10000;
    int* h_indexOfClosest = new int[count];
    float3* h_points = new float3[count];

    // Defining random points
    for (int i = 0; i < count; i++) {
        h_points[i].x = (float)(((rand() % 10000)) - 5000);
        h_points[i].y = (float)(((rand() % 10000)) - 5000);
        h_points[i].z = (float)(((rand() % 10000)) - 5000);
    }

    // Device pointers
    int* d_indexOfClosest;
    float3* d_points;

    // Allocating memory on the device
    cudaMalloc(&d_indexOfClosest, sizeof(int) * count);
    cudaMalloc(&d_points, sizeof(float3) * count);

    // Copying values from the host to the device
    cudaMemcpy(d_points, h_points, sizeof(float3) * count, cudaMemcpyHostToDevice);

    int threads_per_block = 64;
    cout << "Running brute force nearest neighbor on the GPU..." << endl;
    for (int i = 1; i <= 10; i++) {
        long start = clock();

        findClosestGPU<<<(count / threads_per_block) + 1, threads_per_block>>>(d_points, d_indexOfClosest, count);
        cudaDeviceSynchronize();

        // Copying results from the device to the host
        cudaMemcpy(h_indexOfClosest, d_indexOfClosest, sizeof(int) * count, cudaMemcpyDeviceToHost);

        double duration = (clock() - start) / (double)CLOCKS_PER_SEC;
        cout << "Test " << i << " took " << duration << " seconds" << endl;
    }

    // Freeing device memory
    cudaFree(d_indexOfClosest);
    cudaFree(d_points);

    // Freeing host memory
    delete[] h_indexOfClosest;
    delete[] h_points;

    return 0;
}

在GPU上运行几次迭代。请注意，时间比CPU快几个数量级。

它开始时与CPU方法相同，我们创建了一堆我们需要整理的随机点。

//defining parameters
const int count = 10000;
int* d_indexOfClosest = new int[count];
float3* d_points = new float3[count];

//defining random points
for (int i = 0; i < count; i++){
    d_points[i].x = (float)(((rand()%10000))-5000);
    d_points[i].y = (float)(((rand()%10000))-5000);
    d_points[i].z = (float)(((rand()%10000))-5000);
}

然后，我们在vRAM上为点分配一些空间，以及存储结果的地方。

// 设备指针
int* d_indexOfClosest；
float3* d_points；

// 在设备上分配内存
cudaMalloc（&d_indexOfClosest，sizeof（int）*计数）；
cudaMalloc（&d_points，sizeof（float3）*计数）；

我们将点从RAM复制到vRAM

// Copying values from the host to the device
cudaMemcpy(d_points, h_points, sizeof(float3) * count, cudaMemcpyHostToDevice);

然后我们运行内核

findClosestGPU<<<(count / threads_per_block) + 1, threads_per_block>>>(d_points, d_indexOfClosest, count);

内核本身只需要进行一些修改即可将其从CPU转换为GPU。

// Brute force implementation, parallelized on the GPU
__global__ void findClosestGPU(float3* points, int* indices, int count) {
    if (count <= 1) return;
    int idx = threadIdx.x + blockIdx.x * blockDim.x;
    if (idx < count) {
        float3 thisPoint = points[idx];
        float smallestSoFar = 3.40282e38f;

        for (int i = 0; i < count; i++) {
            if (i == idx) continue;

            float dist_sqr = (thisPoint.x - points[i].x) *
                             (thisPoint.x - points[i].x) +
                             (thisPoint.y - points[i].y) *
                             (thisPoint.y - points[i].y) +
                             (thisPoint.z - points[i].z) *
                             (thisPoint.z - points[i].z);
            if (dist_sqr < smallestSoFar) {
                smallestSoFar = dist_sqr;
                indices[idx] = i;
            }
        }
    }
}

首先，我们可以为每个线程分配一个单独的点，指定为idx。我们可以通过从所有线程块中计算当前正在执行的线程来做到这一点。

int idx = threadIdx.x + blockIdx.x * blockDim.x;

我们可以使用这个唯一的ID将每个线程分配给一个特定的点。因此，线程0负责找到最接近点0的点，线程50找到最接近点50点等。

通常有太多的线程，因为可能没有足够的点来填充最后一个线程块，所以我们只评估有效的idx值

int idx = threadIdx.x + blockIdx.x * blockDim.x;
if (idx < count) {
  ...

其余的几乎完全相同，除了指针现在指向vRAM中的空格，而不是RAM。

一旦内核完成执行，所有值都更新，我们可以将它们复制回主机，然后用它们做任何我们想做的事情。

// 将结果从设备复制到主机
cudaMemcpy（h_indexOfClosest，d_indexOfClosest，sizeof（int）* count，cudaMemcpyDeviceToHost）；

就这样，我们在GPU上做了一些诚实到好的工作。我们几乎已经准备好从头开始构建人工智能模型，但首先，我想简要介绍一个超级重要的话题。

英伟达分析器

在CUDA中实施复杂的并行程序时，了解哪些类型的操作对性能的影响最大是有用的。

我们可以%% writefile findClosestGPU.cu将我们的CUDA代码块写为文本文件，而不是在代码块顶部调用%%cuda。然后，我们可以使用nvprof来分析用于执行CUDA代码的计算。

！nvprof ./findClosestGPU.out

nvprof在我们CUDA代码上的结果

可以看出，用内存做事情占用了该应用程序（cudaMalloc和cudaDeviceSynchronize）超过99.19%的时间。在本教程中，我们不会进行优化，但一般来说，在进行CUDA编程时，大部分时间都花在优化内存和设备间通信上，而不是计算（这就是Flash注意力如何将尖端人工智能的加速10倍）。

在本教程的其余部分，我们实际上不会使用分析器，但它是实际构建自己的CUDA程序的关键工具。现在，让我们在CUDA中从头开始构建人工智能。

加入IAEE

定义“从零开始的人工智能”

当然，我不会出去开采硅胶，这样我就可以自己手工制作的电脑芯片了。我们将依靠许多预先存在的技术来在本教程中实现人工智能。从铁杆低级开发人员的角度来看，CUDA实际上是一个高级编程工具。然而，与大多数数据科学家在日常生活中所做的相比，我们基本上是从基岩开始工作。

为了在CUDA中“从头开始”使人工智能，我们将实施一些事情：

• 我们需要从定义一些实用程序开始。我们将构建一个数据结构，帮助我们跟踪神经网络中数据的形状，我们将构建CUDA中错误的抽象，以帮助我们进行调试，我们将构建矩阵的抽象，我们将构建二进制交叉熵。
• 一旦我们所有的实用程序都设置好了，我们将构建一些构成模型本身的类。模型层的通用类，以及线性层、乙状体激活和ReLu激活的实现。
• 在那之后，我们将构建、训练和测试我们的模型。

如果这听起来像很多工作，那就是。但是，如果您是一名数据科学家，我想您会发现，在如此低的水平上查看这些实现在概念上是非常有帮助的，如果您对数据科学完全陌生，并且不知道什么是“ReLu”激活，您将学到很多东西。

在我们开始制作这些东西之前，我想简要讨论一下我们将如何编码此功能。

完整的代码可以在这里找到

这个项目中C++的结构

之前，我们在Google Colab的单元格块顶部使用%%cuda来运行CUDA代码。这是一个简单的方法，但这意味着我们需要在一个单元格块中运行所有代码。展望未来，我们将使用单元格块使用%%writefile魔术编写文本文件，然后使用nvcc显式编译代码。

每当我们定义一些功能时，我们都会首先定义一个头文件。以下是头文件的任意示例：

%%writefile someClass.hh

// 这是如此，如果someClass在多个次中多次导入
//文档，它实际上只被导入一次。
#pragma 一次

类 ClassWithFunctionality {
    // 定义私人物品供内部使用
私人：
    // 定义私人数据
    int someValue；
    另一个价值；

    // 定义私有函数
    void privateFunction1();
    void privateFunction2();

    // 定义在对象之外可以访问的东西
公众：
    // 定义公共数据
    int somePublicValue；
    int someOtherPublicValue；

    // 定义公共函数
ClassWithFunctionality（int构造函数输入）；
    void doSomething1();
    void doSomething2();
};

如果你对此不熟悉，不要太担心。这是C++的东西。我们定义了一个标题文件，以基本上阻止某些东西将具有哪些功能。然后，我们可以在相应的CUDA文件中实现该内容：

%%writefile someClass.cu

#include "someClass.hh"

// defining constructor
ClassWithFunctionality::ClassWithFunctionality(int constructorInput)
    : someValue(constructorInput), anotherValue(2), somePublicValue(3), someOtherPublicValue(4)
{}

void ClassWithFunctionality::doSomething1() {
    return;
}

void ClassWithFunctionality::doSomething2() {
    return;
}

void ClassWithFunctionality::privateFunction1() {
    return;
}

void ClassWithFunctionality::privateFunction2() {
    return;
}

如果你来自Python世界，你可能会认为构造函数语法有点奇怪。这使用一种称为“初始化列表”的语法。基本上，当构建类的实例时，这行代码设置了所有内容的值。

someValue(constructorInput), anotherValue(2), somePublicValue(3), someOtherPublicValue(4)

现在我们同时拥有了标头和CUDA文件，我们实际上无法运行它，因为没有main函数。这很有用，因为我们可以将其导入具有main函数的代码中（每个程序只能有一个），但这意味着我们需要快速制作一个，以便我们可以测试我们的功能。

%%writefile main.cu
#include 
#include "someClass.hh"

// testing SomeClass
int main(void) {
    ClassWithFunctionality example(3);
    std::cout << "it works!" << std::endl;
    return 0;
}

然后，要实际运行此代码，我们可以编译代码，运行它，将输出保存到一个文件中，然后使用以下代码块打印出文件。

！nvcc someClass.cu main.cu -o outputFile.out
!./outputFile.out

现在我们已经解决了这个问题，我们可以定义我们的第一个功能。

实用程序1：形状

我们将实施的人工智能大量利用了2D形状的矩阵。首先，我们可以定义一个名为Shape的数据结构，我们可以用它来跟踪2D大小。

%%writefile shape.hh
#pragma once

struct Shape {
 size_t x, y;

 Shape(size_t x = 1, size_t y = 1);
};

%%writefile shape.cu
#include "shape.hh"

Shape::Shape(size_t x, size_t y) :
 x(x), y(y)
{ }

%%writefile main.cu
#include "shape.hh"
#include 
#include 

using namespace std;

//testing
int main( void ) {
    Shape shape = Shape(100, 200);
    cout << "shape x: " << shape.x << ", shape y: " << shape.y << endl;
}

！nvcc shape.cu main.cu -o shape.out
!./shape.out

实用程序2：NNException

如果GPU存在问题，该问题可能需要一段时间才能传播回CPU。这可能会使调试CUDA程序变得困难，因为错误可能会在奇怪的时期抛出错误。为了缓解这种情况，我们可以使用cudaGetLastError()来检查GPU上的最新错误。NNException是一个围绕cudaGetLastError构建的轻量级包装器，允许我们在整个代码中检查GPU上是否有错误。

%%writefile nn_exception.hh
#pragma once

#include 
#include 

class NNException : std::exception {
private:
 const char* exception_message;

public:
 NNException(const char* exception_message) :
  exception_message(exception_message)
 { }

 virtual const char* what() const throw()
 {
  return exception_message;
 }

 static void throwIfDeviceErrorsOccurred(const char* exception_message) {
  cudaError_t error = cudaGetLastError();
  if (error != cudaSuccess) {
   std::cerr << error << ": " << exception_message;
   throw NNException(exception_message);
  }
 }
};

%%writefile main.cu
//With error handling

#include "nn_exception.hh"
#include 

int main() {
    // Allocate memory on the GPU
    float* d_data;
    cudaError_t error = cudaMalloc((void**)&d_data, 100 * sizeof(float));

    // Check for CUDA errors and throw an exception if any
    try {
        NNException::throwIfDeviceErrorsOccurred("Failed to allocate GPU memory");
    } catch (const NNException& e) {
        std::cerr << "Caught NNException: " << e.what() << std::endl;
        return -1; // Return an error code
    }

    // Free the GPU memory
    error = cudaFree(d_data);

    // Check for CUDA errors again
    try {
        NNException::throwIfDeviceErrorsOccurred("Failed to free GPU memory");
    } catch (const NNException& e) {
        std::cerr << "Caught NNException: " << e.what() << std::endl;
        return -1; // Return an error code
    }

    std::cout << "CUDA operations completed successfully" << std::endl;
    return 0; // Return success
}

！nvcc main.cu shape.cu -o nnexception.out
!./nnexception.out

实用程序3：矩阵

该类抽象了设备和主机之间的一些通信，允许在内存位置之间轻松传递值矩阵。它允许：

• 在GPU上为矩阵分配内存
• 在CPU上为矩阵分配的内存
• 在CPU和GPU上为矩阵分配内存
• 分配内存，如果还没有分配的话
• 将数据从CPU RAM复制到GPU VRAM
• 将数据从GPU VRAM复制到CPU RAM
• 覆盖以允许矩阵像数组一样索引

%%writefile矩阵.hh
#pragma once

#include "shape.hh"

#include 

class Matrix {
private:
 bool device_allocated;
 bool host_allocated;

 void allocateCudaMemory();
 void allocateHostMemory();

public:
 Shape shape;

 std::shared_ptr data_device;
 std::shared_ptr data_host;

 Matrix(size_t x_dim = 1, size_t y_dim = 1);
 Matrix(Shape shape);

 void allocateMemory();
 void allocateMemoryIfNotAllocated(Shape shape);

 void copyHostToDevice();
 void copyDeviceToHost();

 float& operator[](const int index);
 const float& operator[](const int index) const;
};

%%writefile matrix.cu
#include "matrix.hh"
#include "nn_exception.hh"

using namespace std;

Matrix::Matrix(size_t x_dim, size_t y_dim) :
 shape(x_dim, y_dim), data_device(nullptr), data_host(nullptr),
 device_allocated(false), host_allocated(false)
{ }

Matrix::Matrix(Shape shape) :
 Matrix(shape.x, shape.y)
{ }

void Matrix::allocateCudaMemory() {
 if (!device_allocated) {
  float* device_memory = nullptr;
  cudaMalloc(&device_memory, shape.x * shape.y * sizeof(float));
  NNException::throwIfDeviceErrorsOccurred("Cannot allocate CUDA memory for Tensor3D.");
  data_device = std::shared_ptr(device_memory,
            [&](float* ptr){ cudaFree(ptr); });
  device_allocated = true;
 }
}

void Matrix::allocateHostMemory() {
 if (!host_allocated) {
  data_host = std::shared_ptr(new float[shape.x * shape.y],
             [&](float* ptr){ delete[] ptr; });
  host_allocated = true;
 }
}

void Matrix::allocateMemory() {
 allocateCudaMemory();
 allocateHostMemory();
}

void Matrix::allocateMemoryIfNotAllocated(Shape shape) {
 if (!device_allocated && !host_allocated) {
  this->shape = shape;
  allocateMemory();
 }
}

void Matrix::copyHostToDevice() {
 if (device_allocated && host_allocated) {
  cudaMemcpy(data_device.get(), data_host.get(), shape.x * shape.y * sizeof(float), cudaMemcpyHostToDevice);
  NNException::throwIfDeviceErrorsOccurred("Cannot copy host data to CUDA device.");
 }
 else {
  throw NNException("Cannot copy host data to not allocated memory on device.");
 }
}

void Matrix::copyDeviceToHost() {
 if (device_allocated && host_allocated) {
  cudaMemcpy(data_host.get(), data_device.get(), shape.x * shape.y * sizeof(float), cudaMemcpyDeviceToHost);
  NNException::throwIfDeviceErrorsOccurred("Cannot copy device data to host.");
 }
 else {
  throw NNException("Cannot copy device data to not allocated memory on host.");
 }
}

float& Matrix::operator[](const int index) {
 return data_host.get()[index];
}

const float& Matrix::operator[](const int index) const {
 return data_host.get()[index];
}

%%writefile main.cu
#include 
#include "matrix.hh"
#include "nn_exception.hh"

int main() {
    // Create a Matrix object with dimensions 10x10
    Matrix matrix(10, 10);

    // Allocate memory on both host and device
    matrix.allocateMemory();
    std::cout << "Memory allocated on host and device." << std::endl;

    // Initialize host data
    for (size_t i = 0; i < 100; ++i) {
        matrix[i] = static_cast(i);
    }
    std::cout << "Host data initialized." << std::endl;

    // Copy data from host to device
    matrix.copyHostToDevice();
    std::cout << "Data copied from host to device." << std::endl;

    // Clear host data
    for (size_t i = 0; i < 100; ++i) {
        matrix[i] = 0.0f;
    }
    std::cout << "Host data cleared." << std::endl;

    // Copy data back from device to host
    matrix.copyDeviceToHost();
    std::cout << "Data copied from device to host." << std::endl;

    // Verify the data
    bool success = true;
    for (size_t i = 0; i < 100; ++i) {
        if (matrix[i] != static_cast(i)) {
            success = false;
            break;
        }
    }

    if (success) {
        std::cout << "Test passed: Data verification successful." << std::endl;
    } else {
        std::cout << "Test failed: Data verification unsuccessful." << std::endl;
    }

    return 0;
}

！nvcc main.cu matrix.cu shape.cu -o matrix.out
!./matrix.out

定义损失函数

现在我们达到了我们的第一个机器学习概念。

在本文结束时，我们将就一个令人难以置信的简单问题训练我们的人工智能模型。我们将在2D空间中创建一堆随机点，如果这些点在左下角或右上象限，它们应该被分配一个1的标签，如果它们在右下角或左上象限，它们应该有一个标签0

模型将试图解决的建模问题。

因此，我们将把我们的模型拟合到“二进制分类问题”。给定一点，我们希望模型给出两个值之一，1或0，给定一些二维向量。

为了训练我们的模型，我们需要实现一种叫做“损失函数”的东西。损失函数告诉模型在输出预测时有多错误。这个想法是根据模型的错误程度来更新模型，这样在那个特定的例子中就不会那么错误了。在查看了许多点并更新模型以减少错误后，模型应该开始表现良好。

损失函数的一个重要特征是不仅要弄清楚模型有多错误，还要弄清楚什么输出会更好。这是必要的，这样我们才能弄清楚如何更新模型参数，以产生更好的输出。

我不会为这个损失函数太努力地计算，我计划在另一篇文章中介绍它。基本上，有一个称为“二进制交叉熵”的函数，给定一个二进制预测，可以根据已知值计算出该预测的错误程度。

二进制交叉熵，其中y是预测（0-1），x是实际值（0-1）。来自PyTorch文档。

至关重要的是，这个函数是“可微的”，这意味着我们也可以计算增加预测概率是否会带来更好的预测。我计划在另一篇文章中更彻底地介绍这方面的数学，但现在我们可以将这些方程视为理所当然，并在GPU上实现它们。

CUDA中的二进制交叉熵

首先，我们可以看看我们的标题文件，以了解需要什么。注意，有时“损失”也称为“成本”，这就是为什么这个特定函数被称为BCECost。在机器学习中，“损失”和“成本”这两个术语是可以互换的，并代表相同的概念。

%%writefile bce_cost.hh
#pragma once
#include "matrix.hh"

class BCECost {
public:
 float cost(Matrix predictions, Matrix target);
 Matrix dCost(Matrix predictions, Matrix target, Matrix dY);
};

要实现二进制交叉熵进行建模，我们只需要实现两个函数：

• 一个是计算给定输出列表的预测列表的损失（或成本）
• 一个来计算损失的导数。因此，本质上是一个“这个预测应该更大，这个预测应该更小”的列表。

我们可以首先将这些功能中的每一个实现为可以在GPU上执行的内核，允许我们并行进行此计算。此内核计算损失：

__global__ void binaryCrossEntropyCost(float* predictions, float* target,
                                       int size, float* cost) {
    int index = blockIdx.x * blockDim.x + threadIdx.x;

    if (index < size) {
        // Clamp predictions to avoid log(0)
        float pred = predictions[index];
        pred = fmaxf(fminf(pred, 1.0f - 1e-7), 1e-7);

        float partial_cost = target[index] * logf(pred)
                + (1.0f - target[index]) * logf(1.0f - pred);
        atomicAdd(cost, - partial_cost / size);
    }
}

这个内核计算损失的导数：

__global__ void dBinaryCrossEntropyCost(float* predictions, float* target, float* dY,
                                        int size) {
    int index = blockIdx.x * blockDim.x + threadIdx.x;

    if (index < size) {
        // Clamp predictions to avoid division by zero
        float pred = predictions[index];
        pred = fmaxf(fminf(pred, 1.0f - 1e-7), 1e-7);

        dY[index] = -1.0 * (target[index] / pred - (1 - target[index]) / (1 - pred));
    }
}

每个内核都需要一个指针，分别用于存储输出、float* cost和float* dY。损失函数只需将二进制交叉熵函数应用于所有点，并在cost位置将它们相加。二进制交叉熵函数的导数在所有点上遍发生，并计算该点的导数，并存储该值是否应该在dY中更高或更低。

这些内核中的每一个实际上都由BCECost类的cost和dCost方法调用

//cost (or loss)
float BCECost::cost(Matrix predictions, Matrix target) {
 assert(predictions.shape.x == target.shape.x);

 float* cost;
 cudaMallocManaged(&cost, sizeof(float));
 *cost = 0.0f;

 dim3 block_size(256);
 dim3 num_of_blocks((predictions.shape.x + block_size.x - 1) / block_size.x);
 binaryCrossEntropyCost<<>>(predictions.data_device.get(),
                target.data_device.get(),
                predictions.shape.x, cost);
 cudaDeviceSynchronize();
 NNException::throwIfDeviceErrorsOccurred("Cannot compute binary cross entropy cost.");

 float cost_value = *cost;
 cudaFree(cost);

 return cost_value;
}

//derivative of cost (aka derivative of loss)
Matrix BCECost::dCost(Matrix predictions, Matrix target, Matrix dY) {
 assert(predictions.shape.x == target.shape.x);

 dim3 block_size(256);
 dim3 num_of_blocks((predictions.shape.x + block_size.x - 1) / block_size.x);
 dBinaryCrossEntropyCost<<>>(predictions.data_device.get(),
                 target.data_device.get(),
                 dY.data_device.get(),
                 predictions.shape.x);
 NNException::throwIfDeviceErrorsOccurred("Cannot compute derivative for binary cross entropy.");

 return dY;
}

除了一个关键差异外，这两个函数的行为通常相似。cost创建自己的值并返回它，而dCost期望它输出dY作为输入并修改它。我们将在后面的章节中探讨为什么会更多。现在，我们可以通过给它一些数据来玩弄我们的损失函数。

%%writefile main.cu
#include 
#include 
#include "matrix.hh"
#include "bce_cost.hh"
#include "nn_exception.hh"

// Helper function to initialize a Matrix with data
void initializeMatrix(Matrix& matrix, const std::vector& data) {
    for (size_t i = 0; i < data.size(); ++i) {
        matrix[i] = data[i];
    }
    matrix.copyHostToDevice();
}

int main() {
    // Define the size of the data
    const int size = 10;

    // Create predictions and target data
    std::vector predictions_data = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95};
    std::vector target_data = {0, 0, 1, 0, 1, 0, 1, 1, 1, 0};

    // Create Matrix objects for predictions and targets
    Matrix predictions(size, 1);
    Matrix target(size, 1);
    predictions.allocateMemory();
    target.allocateMemory();

    // Initialize matrices with data
    initializeMatrix(predictions, predictions_data);
    initializeMatrix(target, target_data);

    // Compute the binary cross-entropy cost
    BCECost bce_cost;
    float cost_value = bce_cost.cost(predictions, target);
    std::cout << "Binary Cross-Entropy Cost: " << cost_value << std::endl;

    // Compute the gradient of the binary cross-entropy cost
    Matrix dY(size, 1);
    dY.allocateMemory();
    Matrix dCost_matrix = bce_cost.dCost(predictions, target, dY);
    dCost_matrix.copyDeviceToHost();

    // Print the gradient values
    std::cout << "Gradient of Binary Cross-Entropy Cost: ";
    for (int i = 0; i < size; ++i) {
        std::cout << dCost_matrix[i] << " ";
    }
    std::cout << std::endl;

    return 0;
}

！nvcc main.cu matrix.cu shape.cu bce_cost.cu -o bce.out
!./bce.out

注意：梯度是派生词的花哨词。

在这个测试中，我们给出了损失函数这些预测值：

std::vector predictions_data = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95};

和这些期望值：

std::vector target_data = {0, 0, 1, 0, 1, 0, 1, 1, 1, 0};

当我们将这两个值都输入我们的导数损失函数时，我们得到了这些结果：

1.11, 1.25 -3.33, 1.66, -2.00, 2.50, -1.43, -1.25, -1.11, 20 

请注意，当预测值大于目标值时，导数的结果是正的。当预测值低于目标值时，导数为负数。此外，预测越错误，导数的幅度（正或负）就越大。

损失衍生品的这种一般质量至关重要，因为它允许我们根据小问题以小方式更新模型，根据大问题以大方式更新模型，并根据模型错误的方式以正确的方式更新模型。

定义模型

在上一节中，我们定义了我们的问题。给定二维空间中的一个点，我们的模型应该预测它是在左下角还是右上象限（输出为0），还是在右下角或左上象限（输出为1）

模型将试图解决的建模问题。

为了真正解决这个问题，我们将使用神经网络。具体来说，我们将使用这个神经网络：

//the neural network we'll ultimately define
NeuralNetwork nn;
nn.addLayer(new LinearLayer("linear_1", Shape(2, 30)));
nn.addLayer(new ReLUActivation("relu_1"));
nn.addLayer(new LinearLayer("linear_2", Shape(30, 1)));
nn.addLayer(new SigmoidActivation("sigmoid_output"));

在我们开始编码之前，让我们先了解一下每个层应该做什么。

线性层

这一层确实是神经网络得名的方式。线性网络由一系列节点和边缘组成，每个边缘都有一些权重，每个节点都有一些偏差。当值通过线性层时，值乘以权重并添加偏置。

ReLU激活

ReLU激活函数是一个常见的激活函数，旨在为线性层添加非线性。通过将ReLU等非线性函数注入线性网络，它允许线性网络学习更复杂的关系。

ReLU，这只是一条整流线。任何输入（x轴）都映射到输出（y轴）

Sigmoid激活

Sigmoid与ReLU相似，只是它经常被用作概率输出的激活函数。我们将使用此激活功能来创建最终输出。

Sigmoid，它只是一条整线。任何输入（x轴）都映射到输出（y轴）

定义一个图层

为了使我们的生活更轻松，我们将实现一个名为NNLayer的抽象类，所有层都将继承。

%%writefile nn_layer.hh
#pragma once

#include 
#include "matrix.hh"

class NNLayer {
protected:
 std::string name;

public:
 virtual ~NNLayer() = 0;

 virtual Matrix& forward(Matrix& A) = 0;
 virtual Matrix& backprop(Matrix& dZ, float learning_rate) = 0;

 std::string getName() { return this->name; };

};

inline NNLayer::~NNLayer() {}

每个NNLayer将实现两个关键功能：

• 前进：给定对网络的输入，该层将计算传递给下一层的值，最终产生预测。
• backprop：鉴于损失函数的导数所指示的是更好的输出，该函数告知前一层它应该如何更新其参数以产生更好的输出。

这两个步骤是现代机器学习的标志，几乎存在于每个尖端人工智能模型中。

培训过程的概念图。来自我关于LoRA的文章。

让我们实现我们的第一层，基于NNLayer

线性层

线性层是神经网络获得其名称的方式，可能是项目中最复杂的模块。基本上，你可以把它想象成一堆节点（受人类大脑中神经元的启发），它们相互插入，产生输出。

回想一下之前线性层的概念图。左边的小图是连续三个线性层。我们正在实现单个线性层的功能。

这些连接中的每一个都有一个权重，每个神经元都有一个偏差。通过更新模型的权重和偏差，模型可以学习过滤输入以获得一些输出。对于给定的神经元，每个输入值乘以相应的权重，所有这些值都相加，然后添加偏置项。

这正是用于线性网络前向传递的CUDA内核所做的。

//in this code W is the weight matrix, A is the input, Z is the output
//and b is the bias, then there are dimensions for the weight matrix
//and the input matrix.
__global__ void linearLayerForward( float* W, float* A, float* Z, float* b,
                                    int W_x_dim, int W_y_dim,
                                    int A_x_dim, int A_y_dim) {

    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    int Z_x_dim = A_x_dim;
    int Z_y_dim = W_y_dim;

    float Z_value = 0;

    if (row < Z_y_dim && col < Z_x_dim) {
        for (int i = 0; i < W_x_dim; i++) {
            Z_value += W[row * W_x_dim + i] * A[i * A_x_dim + col];
        }
        Z[row * Z_x_dim + col] = Z_value + b[row];
    }
}

线性层也需要一个用于反向传播的核心，所以如果我们知道线性层的输出应该如何变化，我们就可以弄清楚线性网络的输入应该如何变化。

//in this code W is the weight matrix, A is the input, dZ is the output
//of the linear layer should change, dA is how the input of the linear
//layer should change, then there are dimensions for the weight matrix
//and the input matrix.
__global__ void linearLayerBackprop(float* W, float* dZ, float *dA,
                                    int W_x_dim, int W_y_dim,
                                    int dZ_x_dim, int dZ_y_dim) {

    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;

    // W is treated as transposed
    int dA_x_dim = dZ_x_dim;
    int dA_y_dim = W_x_dim;

    float dA_value = 0.0f;

    if (row < dA_y_dim && col < dA_x_dim) {
        for (int i = 0; i < W_y_dim; i++) {
            dA_value += W[i * W_x_dim + row] * dZ[i * dZ_x_dim + col];
        }
        dA[row * dA_x_dim + col] = dA_value;
    }
}    

对此的计算非常直观。如果输出dZ应该变大，那么无论将输入与输出连接的重量如何，输入都应该变大dZ倍。输入和输出之间的线性关系，即线的斜率是权重，这就是为什么线性层被称为线性层。

输入的变化是所有连接的输出加权变化的总和，如for循环中例证。

for (int i = 0; i < W_y_dim; i++) {
  dA_value += W[i * W_x_dim + row] * dZ[i * dZ_x_dim + col];
}

为了实际更新模型的参数，我们可以创建两个额外的内核。一个更新模型的权重，另一个更新偏差。

//in this code W is the weight matrix, A is the input, dZ is the output
//of the linear layer should change, dA is how the input of the linear
//layer should change, then there are dimensions for the weight matrix
//and the input matrix.
__global__ void linearLayerUpdateWeights(  float* dZ, float* A, float* W,
                                           int dZ_x_dim, int dZ_y_dim,
                                           int A_x_dim, int A_y_dim,
                                           float learning_rate) {

    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;

    // A is treated as transposed
    int W_x_dim = A_y_dim;
    int W_y_dim = dZ_y_dim;

    float dW_value = 0.0f;

    if (row < W_y_dim && col < W_x_dim) {
        for (int i = 0; i < dZ_x_dim; i++) {
            dW_value += dZ[row * dZ_x_dim + i] * A[col * A_x_dim + i];
        }
        W[row * W_x_dim + col] = W[row * W_x_dim + col] - learning_rate * (dW_value / A_x_dim);
    }
}

//in this code W is the weight matrix, A is the input, dZ is the output
//of the linear layer should change, dA is how the input of the linear
//layer should change, then there are dimensions for the weight matrix
//and the input matrix.
__global__ void linearLayerUpdateBias(  float* dZ, float* b,
                                        int dZ_x_dim, int dZ_y_dim,
                                        int b_x_dim,
                                        float learning_rate) {
    int index = blockIdx.x * blockDim.x + threadIdx.x;

    if (index < dZ_x_dim * dZ_y_dim) {
        int dZ_x = index % dZ_x_dim;
        int dZ_y = index / dZ_x_dim;
        atomicAdd(&b[dZ_y], - learning_rate * (dZ[dZ_y * dZ_x_dim + dZ_x] / dZ_x_dim));
    }
}

这些在功能上都相同。如果损失的反向传播表明模型的输出应该更大，那么该输出的权重和偏差都试图使输出更大。如果输出应该更小，权重和偏差会更新，以便值更小。

这真的是更新权重的关键代码

for (int i = 0; i < dZ_x_dim; i++) {
    dW_value += dZ[row * dZ_x_dim + i] * A[col * A_x_dim + i];
}
W[row * W_x_dim + col] = W[row * W_x_dim + col] - learning_rate * (dW_value / A_x_dim);

您可能会想，“如果一个权重对应一个输出，为什么我们需要一起对一个权重进行迭和添加更改”？原因是我们正在分批。批处理的想法是，当您去训练模型时，您向其提供一批输入，并使用整个批处理的结果来通知参数更改。这使得模型的更新质量更高，计算效率更高。

如果批次大小为32，则根据这32个预测中的每个预测情况来更新每个权重。这就是为什么我们需要遍特dZ_x_dim因为这是批处理维度。

对于批次中的一个特定示例，我们通过将输出应变化的量乘以输入的大小来计算权重应该如何变化。

dW_value += dZ[row * dZ_x_dim + i] * A[col * A_x_dim + i];

为什么这是乘法的数学有点超出本文的范围，但它与dZ表示相对于线性层输出的损失梯度有关。有一种东西叫做连锁规则，这是一个整体，我将在以后的文章中围绕这个话题讨论更多数学。目前，我们可以认为这是理所当然的，并使用此表达式根据特定的批次示例来计算权重应该如何变化。

一旦我们了解了权重应该如何变化（dW_value），我们可以用以下表达式更新我们的权重：

W[row * W_x_dim + col] = W[row * W_x_dim + col] - learning_rate * (dW_value / A_x_dim);

这个表达中有两个重要概念。第一，我们使用学习率来更改我们更新模型的速率。这是机器学习中一个非常有用的参数，允许我们根据单个训练批次来决定模型本身的更新程度。第二个关键想法是，我们还将dW_value批次数（A_x_dim）。我们这样做是为了将所有批次的累积转化为平均值。

我们可以使用这些内核来实现线性层的主要功能，但需要实现一堆其他锅炉板才能使线性网络正常工作，从头文件可以看出。

%%writefile linear_layer.hh
#pragma once
#include "nn_layer.hh"


class LinearLayer : public NNLayer {
private:
 const float weights_init_threshold = 0.01;

 Matrix W;
 Matrix b;

 Matrix Z;
 Matrix A;
 Matrix dA;

 void initializeBiasWithZeros();
 void initializeWeightsRandomly();

 void computeAndStoreBackpropError(Matrix& dZ);
 void computeAndStoreLayerOutput(Matrix& A);
 void updateWeights(Matrix& dZ, float learning_rate);
 void updateBias(Matrix& dZ, float learning_rate);

public:
 LinearLayer(std::string name, Shape W_shape);
 ~LinearLayer();

 Matrix& forward(Matrix& A);
 Matrix& backprop(Matrix& dZ, float learning_rate = 0.01);

 int getXDim() const;
 int getYDim() const;

 Matrix& getWeightsMatrix();
    Matrix& getBiasVector();
};

我们可以通过此功能快速闪电战：

• initializeBiasWithZeros：此函数将线性网络中所有节点的偏置参数设置为零。
• 初始化WeightsRandomly：此函数将线性网络的所有权重参数设置为值在0和1之间的随机正态分布。
• computeAndStoreBackpropError：调用 linearLayerBackprop 内核。
• computeAndStoreLayerOutput：调用 linearLayerForward 内核。
• updateWeights：调用 linearLayerUpdateWeights kenrel。
• updateBias：调用 linearLayerUpdateBias 内核。

如果您对实现感到好奇，请在此处查看代码。

测试培训

我们还有很多工作要做，但在这一点上，我们有一个线性模型和一个损失函数，我们可以玩弄训练一个琐碎的模型。

以下代码创建一个输入[0.1, 0.2, 0.3]，然后训练模型以输出target[0] = 1.0。我们可以玩弄目标，看看我们是否可以得到一个模型来学习正确的输出。

%%writefile main.cu
#include "linear_layer.hh"
#include "bce_cost.hh"
#include 
#include "matrix.hh"

void printMatrix(Matrix& matrix, const std::string& name) {
    matrix.copyDeviceToHost();
    std::cout << name << ":" << std::endl;
    for (int i = 0; i < matrix.shape.x * matrix.shape.y; ++i) {
        std::cout << matrix[i] << " ";
    }
    std::cout << std::endl;
}

int main() {
    // Define input dimensions and initialize the layer
    Shape input_shape(1, 3); // (1 rows, 3 columns, transposed vector)
    Shape weight_shape(3, 1); // shape of weights, resulting in a 1x1 output

    LinearLayer layer("test_layer", weight_shape);

    // Allocate memory for input and output
    Matrix input(input_shape);
    input.allocateMemory();
    input[0] = 0.1f; input[1] = 0.2f; input[2] = 0.3f;
    input.copyHostToDevice();

    // Allocate memory for target
    Matrix target(Shape(1, 1)); // 1x1 target matrix
    target.allocateMemory();
    target[0] = 1.0f;
    target.copyHostToDevice();

    // Print initial weights and biases
    printMatrix(layer.getWeightsMatrix(), "Initial Weights");
    printMatrix(layer.getBiasVector(), "Initial Biases");

    // Training loop
    for (int i = 0; i < 10; ++i) {
        // Perform forward pass
        Matrix& output = layer.forward(input);
        output.copyDeviceToHost();

        // Print forward pass output
        std::cout << "Forward pass output:" << std::endl;
        for (int j = 0; j < output.shape.x * output.shape.y; ++j) {
            std::cout << output[j] << " ";
        }
        std::cout << std::endl;

        // Calculate BCE loss
        BCECost bce;
        float loss = bce.cost(output, target);
        std::cout << "Loss at iteration " << i << ": " << loss << std::endl;

        // Calculate gradient of BCE loss
        Matrix dZ(output.shape);
        dZ.allocateMemory();
        bce.dCost(output, target, dZ);

        // Perform backpropagation
        float learning_rate = 0.000001f;
        layer.backprop(dZ, learning_rate);
    }

    // Print updated weights and biases
    printMatrix(layer.getWeightsMatrix(), "Updated Weights");
    printMatrix(layer.getBiasVector(), "Updated Biases");

    return 0;
}

！nvcc main.cu matrix.cu shape.cu bce_cost.cu linear_layer.cu -o ll.out
!./ll.out

让我们训练模型对几个target值进行预测：

预测1.0的训练

预测0.0的训练

从我们之前的二进制交叉熵实现中可以看出，我们正在将预测剪切到零和1之间，这样就不会有奇怪的零误差除法。

//from the implementation of binary cross entropy
float pred = predictions[index];
pred = fmaxf(fminf(pred, 1.0f - 1e-7), 1e-7);

因此，就我们对二元交叉熵的定义而言，小于0的预测在功能上等价于0，高于1的预测在功能上等价于1。因此，这意味着我们的线性网络正在正确学习在本例中预测1和0。

乙状体激活功能

使密集网络学会预测复杂事物的事情之一是激活函数的存在。如果我们只是背靠背地串起一堆线性层，输入将始终与输出有线性关系，并且永远无法表达比线更复杂的关系。我们在混合中抛出非线性函数，如sigmoid，这样模型的线性输出就可以变成非线性值。

来自我关于YOLO的手工文章

如果您了解线性网络，激活函数的实现非常简单。没有可学习的参数，它只是在正向和向后传递中非线性映射值。在正向传递中，sigmoid函数使用以下表达式：

乙状体功能。来源。

在向后传递中，它使用以下表达式：

乙状函数的导数。来源。

因此，我们只需要在CUDA内核中实现这些功能，将它们粘在NNLayer版本中，我们非常好。

__device__ float sigmoid(float x) {
    return 1.0f / (1 + exp(-x));
}

__global__ void sigmoidActivationForward(float* Z, float* A,
                                         int Z_x_dim, int Z_y_dim) {

    int index = blockIdx.x * blockDim.x + threadIdx.x;

    if (index < Z_x_dim * Z_y_dim) {
        A[index] = sigmoid(Z[index]);
    }
}

__global__ void sigmoidActivationBackprop(float* Z, float* dA, float* dZ,
                                          int Z_x_dim, int Z_y_dim) {

    int index = blockIdx.x * blockDim.x + threadIdx.x;

    if (index < Z_x_dim * Z_y_dim) {
        dZ[index] = dA[index] * sigmoid(Z[index]) * (1 - sigmoid(Z[index]));
    }
}

这里在概念上没有什么新东西，但请随时查看要挖掘的代码。

ReLU激活功能

ReLU本质上是完全相同的故事，但具有整流的线性激活功能。整换线性只是一条有断裂的线。任何低于x=0的值都等于0，任何高于x=0的值都等于x。因此，这是一条正值的线，0的负值。

__global__ void reluActivationForward(float* Z, float* A,
                                      int Z_x_dim, int Z_y_dim) {

    int index = blockIdx.x * blockDim.x + threadIdx.x;

    if (index < Z_x_dim * Z_y_dim) {
        A[index] = fmaxf(Z[index], 0);
    }
}

__global__ void reluActivationBackprop(float* Z, float* dA, float* dZ,
                                       int Z_x_dim, int Z_y_dim) {

    int index = blockIdx.x * blockDim.x + threadIdx.x;

    if (index < Z_x_dim * Z_y_dim) {
        if (Z[index] > 0) {
            dZ[index] = dA[index];
        }
        else {
            dZ[index] = 0;
        }
    }
} 

和这个一样，我跳过了一点锅炉板，但在概念上没有什么新意。请随时查看代码。

训练我们的模型

最后，我们到了。我们只有几件事要实施。在我们结束之前，让我们评估一下我们到目前为止所做的工作。我们有这些公用事业：

• 处理矩阵的抽象
• 处理错误的抽象
• 一个用于定义形状的方便的小结构
• 二进制交叉熵的实现，包括正向和向后通道

我们已经定义了我们所有的层：

• 完全连接的层，包括向前和向后通过
• 乙状体激活函数，包括正向和向后传递
• ReLU激活功能，包括前进和后退通道

现在我们可以把这些都放在一起来训练一个模型。要做到这一点，我们将再定义两个抽象：

• 一个给模特
• 一个用于数据集

首先，让我们定义模型。实现相当微不足道，您可以参考完整代码了解所有细节。为了本文的目的，我们只需查看标题。

%%writefile neural_network.hh
#pragma once

#include 
#include "nn_layer.hh"
#include "bce_cost.hh"

class NeuralNetwork {
private:
 std::vector layers;
 BCECost bce_cost;

 Matrix Y;
 Matrix dY;
 float learning_rate;

public:
 NeuralNetwork(float learning_rate = 0.01);
 ~NeuralNetwork();

 Matrix forward(Matrix X);
 void backprop(Matrix predictions, Matrix target);

 void addLayer(NNLayer *layer);
 std::vector getLayers() const;

};

该模型包括一系列图层、一个损失函数、一个在通过模型时保存预测的地方（Y）、一个在模型中反向传播时存储梯度的地方（dY、学习率、一个向前传播的函数、一个反向传播的函数，以及一些用于添加图层和打印图层的实用程序。

数据集也很简单。它定义了一堆点，根据它们处于什么象限将它们标记为1或0（这是我们简单的神经网络将尝试学习的内容），然后公开一些功能来获取这些值。

%%writefile coordinates_dataset.hh
#pragma once

#include "matrix.hh"
#include 

class CoordinatesDataset {
private:
 size_t batch_size;
 size_t number_of_batches;

 std::vector batches;
 std::vector targets;

public:

 CoordinatesDataset(size_t batch_size, size_t number_of_batches);

 int getNumOfBatches();
 std::vector& getBatches();
 std::vector& getTargets();

};

有了这一点，我们终于可以实施我们的培训代码：

%%writefile main.cu
#include 
#include 

#include "neural_network.hh"
#include "linear_layer.hh"
#include "relu_activation.hh"
#include "sigmoid_activation.hh"
#include "nn_exception.hh"
#include "bce_cost.hh"

#include "coordinates_dataset.hh"

float computeAccuracy(const Matrix& predictions, const Matrix& targets);

int main() {

    srand( time(NULL) );

    CoordinatesDataset dataset(100, 21);
    BCECost bce_cost;

    NeuralNetwork nn;
    nn.addLayer(new LinearLayer("linear_1", Shape(2, 30)));
    nn.addLayer(new ReLUActivation("relu_1"));
    nn.addLayer(new LinearLayer("linear_2", Shape(30, 1)));
    nn.addLayer(new SigmoidActivation("sigmoid_output"));

    // network training
    Matrix Y;
    for (int epoch = 0; epoch < 1001; epoch++) {
        float cost = 0.0;

        for (int batch = 0; batch < dataset.getNumOfBatches() - 1; batch++) {
            Y = nn.forward(dataset.getBatches().at(batch));
            nn.backprop(Y, dataset.getTargets().at(batch));
            cost += bce_cost.cost(Y, dataset.getTargets().at(batch));
        }

        if (epoch % 100 == 0) {
            std::cout     << "Epoch: " << epoch
                        << ", Cost: " << cost / dataset.getNumOfBatches()
                        << std::endl;
        }
    }

    // compute accuracy
    Y = nn.forward(dataset.getBatches().at(dataset.getNumOfBatches() - 1));
    Y.copyDeviceToHost();

    float accuracy = computeAccuracy(
            Y, dataset.getTargets().at(dataset.getNumOfBatches() - 1));
    std::cout     << "Accuracy: " << accuracy << std::endl;

    return 0;
}

float computeAccuracy(const Matrix& predictions, const Matrix& targets) {
    int m = predictions.shape.x;
    int correct_predictions = 0;

    for (int i = 0; i < m; i++) {
        float prediction = predictions[i] > 0.5 ? 1 : 0;
        if (prediction == targets[i]) {
            correct_predictions++;
        }
    }

    return static_cast(correct_predictions) / m;
}

!nvcc main.cu matrix.cu shape.cu bce_cost.cu sigmoid_activation.cu relu_activation.cu linear_layer.cu coordinates_dataset.cu neural_network.cu -o main.out
!./main.out

简而言之，该代码具有以下作用：

1. 定义我们的数据集和神经网络

CoordinatesDataset dataset(100, 21);
BCECost bce_cost;

NeuralNetwork nn;
nn.addLayer(new LinearLayer("linear_1", Shape(2, 30)));
nn.addLayer(new ReLUActivation("relu_1"));
nn.addLayer(new LinearLayer("linear_2", Shape(30, 1)));
nn.addLayer(new SigmoidActivation("sigmoid_output"));

2.将模型暴露给我们的数据1000次

for (int epoch = 0; epoch < 1001; epoch++) {

3.对于每个纪元，它遍及所有批次的数据，通过模型传递批次以生成预测池，计算损失，并通过模型反向传播以更新参数。

for (int batch = 0; batch < dataset.getNumOfBatches() - 1; batch++) {
    Y = nn.forward(dataset.getBatches().at(batch));
    nn.backprop(Y, dataset.getTargets().at(batch));
    cost += bce_cost.cost(Y, dataset.getTargets().at(batch));
}

4.一旦完成，它会通过模型传递所有数据，并计算模型的准确性。

// compute accuracy
Y = nn.forward(dataset.getBatches().at(dataset.getNumOfBatches() - 1));
Y.copyDeviceToHost();

float accuracy = computeAccuracy(
        Y, dataset.getTargets().at(dataset.getNumOfBatches() - 1));
std::cout     << "Accuracy: " << accuracy << std::endl;

有了这一点，我们在CUDA中从头开始实现了一个简单的神经网络。

我的天啊

即使你没有读这篇文章，我希望你没有伤到手指，一直滚动到这一点。在这篇文章中，我们做了很多。

• 我们简要讨论了为什么CUDA值得学习
• 我们浏览了计算机的所有主要组件，发现主机（CPU和RAM）和设备（GPU和vRAM）本质上是独立的参与者，可以一起工作
• 我们讨论了为什么GPU在人工智能环境中很重要（它们比CPU更擅长运行可并行化任务）。
• 然后，我们引入了CUDA，在那里我们探索了定义CUDA内核和执行启动配置，以及设备和主机之间的通信。
• 我们利用新发现的CUDA知识在GPU上并行化一个简单的CPU程序。
• 我们简要地探索了NVIDIA分析器，它使我们能够了解某些操作在CUDA程序中需要多长时间。
• 我们定义了一些实用程序：形状、NNException、矩阵和二进制交叉熵。
• 我们定义了模型的层：线性、Sigmoid和ReLU。
• 然后，最后，我们训练了我们的模型，并观察到学习简单任务的能力。

在未来的文章中，我将讨论如何将CUDA注入PyTorch，以及一些使用CUDA提高人工智能效率的论文。敬请期待！