K-Means agrupado con Python Numba y CUDA C

Agrupamiento K-Means con Python Numba y CUDA C

Cómo acelerar tu análisis de datos 1600x vs Scikit-Learn (¡con código!)

Paralelizar cargas de trabajo de análisis de datos puede ser una tarea desafiante, especialmente cuando no hay una implementación lista para usar eficiente disponible para tu caso de uso específico. En este tutorial, te guiaré a través de los principios de escritura de kernels CUDA tanto en C como en Python Numba, y cómo esos principios pueden aplicarse al algoritmo clásico de clustering k-means. Al final de este artículo, podrás escribir una implementación personalizada y paralelizada del algoritmo batched k-means tanto en C como en Python, logrando una aceleración de hasta 1600x en comparación con la implementación estándar de scikit-learn. Si quieres ir directamente al código, está disponible en Colab.

Introducción

Estudiaremos dos frameworks para paralelizar en GPUs NVIDIA: Python Numba y CUDA C. Luego implementaremos el mismo algoritmo batched k-means en ambos y los compararemos con scikit-learn.

Numba es una curva de aprendizaje más suave para aquellos que prefieren Python sobre C. Numba compila porciones de tu código en funciones CUDA especializadas llamadas kernels. CUDA C está un nivel más bajo en el árbol de abstracción, lo que ofrece un control más refinado. Los ejemplos de Colab están diseñados para reflejar estrechamente las implementaciones de los algoritmos, por lo que una vez que comprendas uno de ellos, podrás entender fácilmente el otro.

Este tutorial surgió de una implementación personalizada del algoritmo k-means batched que tuve que escribir, y lo utilizaré como ejemplo ilustrativo. Normalmente, las bibliotecas de k-means están optimizadas para ejecutar una sola instancia del algoritmo en un conjunto de datos extremadamente grande. Sin embargo, en nuestro proyecto necesitábamos clusterizar millones de conjuntos de datos pequeños e individuales en paralelo, algo para lo que no pudimos encontrar una biblioteca lista para usar.

• Coraje para aprender ML Desmitificando la regularización L1 y L2 (parte 1)
• Eleva tus análisis de negocios Guía paso a paso para el ajuste estacional
• Valentía para aprender ML Descifrando la Regularización L1 y L2 (parte 2)

Este artículo se divide en cuatro secciones principales: Conceptos básicos de CUDA, Algoritmo k-means paralelizado, Python Numba y CUDA C. Intentaré mantener el material relativamente contenido en sí mismo y repasar brevemente los conceptos según sea necesario.

Conceptos básicos de CUDA

Las CPUs están diseñadas para un procesamiento secuencial rápido y las GPUs están diseñadas para un paralelismo masivo. Para nuestro algoritmo, necesitábamos ejecutar k-means en millones de pequeños conjuntos de datos independientes, lo que se presta perfectamente para la implementación en GPU.

Nuestra implementación utilizará CUDA, que es el acrónimo de Compute Unified Device Architecture, una biblioteca y plataforma de cómputo desarrollada por NVIDIA para aprovechar las GPUs en la computación paralela.

Para comprender los kernels de GPU y las funciones del dispositivo, aquí tienes algunas definiciones:

Un hilo es una secuencia de instrucciones que se puede ejecutar de forma independiente.

Un núcleo es una unidad de procesamiento que puede ejecutar un solo hilo.

Un warp es la unidad más pequeña de programación de hilos. Cada warp consta de 32 hilos y los programa en los núcleos.

Un multiprocesador, también conocido como Streaming Multiprocessor (SM), está compuesto por un número fijo de núcleos.

Un kernel de GPU es una función diseñada para ejecutarse en paralelo en muchos hilos de GPU. Las funciones de kernel son invocadas por la CPU y ejecutadas en la GPU.

Una función del dispositivo es una función que puede ser llamada por una función kernel u otra función del dispositivo. Esta función se ejecuta en la GPU.

Estos kernels se organizan en una jerarquía de grid, blocks y threads. En el contexto de una GPU, cada hilo está conectado a un solo núcleo y ejecuta una copia del kernel.

Un bloque es una colección de hilos que se ejecutan en un multiprocesador.

Una grilla es una matriz abstracta para organizar bloques de hilos. Las grillas se utilizan para asignar las instancias del kernel a los hilos por índice.

En la GPU, existen dos tipos de memoria: memoria global y memoria compartida.

La memoria global se almacena en la DRAM (Memoria de Acceso Aleatorio Dinámica). Todos los hilos pueden acceder a su contenido a costa de una alta latencia de memoria.

La memoria compartida se almacena en la caché y es privada para los hilos dentro del mismo bloque. Se puede acceder mucho más rápido que la memoria global.

Algoritmo K-Means Paralelizado

Ahora vamos a introducir el algoritmo k-means. K-means es un algoritmo de agrupamiento no supervisado que divide un conjunto de datos en k clústeres distintos y no superpuestos. Dado un conjunto de puntos de datos, primero inicializamos k centroides, o los centros iniciales de los clústeres:

Luego, después de elegir los centroides iniciales, realizamos iterativamente dos pasos:

1. Paso de Asignación: Cada punto de datos se asigna a su centroide más cercano basado en la distancia euclidiana.
2. Paso de Actualización: Las ubicaciones de los centroides se reasignan al promedio de todos los puntos asignados a ese centroide en el paso anterior.

Estos pasos se repiten hasta converger, es decir, cuando las ubicaciones de los centroides ya no cambian significativamente. El resultado es el conjunto de coordenadas de los centroides de los k clústeres junto con un arreglo que marca los índices de los clústeres (llamados etiquetas) para cada uno de los puntos de datos originales.

Para un conjunto de datos grande, la inicialización de los centroides puede afectar significativamente el resultado del algoritmo. Por lo tanto, el programa intenta múltiples centroides iniciales, llamados semillas iniciales, y devuelve el resultado de la mejor semilla. Cada semilla se selecciona del conjunto de datos inicial sin repetición, lo que significa que no se repetirá ningún centroide inicial. El número óptimo de semillas para nuestro algoritmo es un tercio del número de puntos de datos. En nuestro programa, debido a que ejecutamos k-means en filas individuales de 100 puntos de datos, el número óptimo de semillas sería 33.

En nuestra función de k-means, las un millón de filas se representan mediante bloques, mientras que los hilos representan las semillas. Los hilos en un bloque están organizados en warps, que son la unidad más pequeña de programación de hilos en la arquitectura del hardware. Cada warp consta de 32 hilos, por lo que es óptimo establecer los tamaños de bloque como múltiplos de 32. Luego, cada bloque devuelve los datos de la semilla que tiene la inercia más pequeña, medida por la suma de la distancia euclidiana entre los centros de los clústeres y sus puntos asignados.

Aquí se encuentra el Colab en el que puedes seguir en Python o C.

variables globales: numRows, lineSize, numClustersdef hostKMeans:    inputData = initializeInputData(numRows, lineSize)    outputCentroids = createEmptyArray(numRows, numClusters)    outputLabels = createEmptyArray(numRows, lineSize)        sendToDevice(inputData, outputCentroids, outputLabels)    cuda_kmeans[blockDimensions, threadDimensions](inputData, outputCentroids, outputLabels)    waitForKernelCompletion()    copyFromDevice(outputCentroids, outputLabels)

En nuestro algoritmo k-means, comenzamos estableciendo variables globales. Necesitaremos hacer referencia a ellas tanto desde la CPU como desde la GPU.

Mientras que las variables globales son accesibles desde el kernel, el kernel no puede devolver valores directamente al host. Para sortear esta limitación, la porción de CPU de nuestro código pasa dos arrays vacíos al kernel junto con los datos de entrada. Estos arrays se utilizarán para copiar los arrays finales de centroides y etiquetas de vuelta a la CPU.

Con nuestras estructuras de datos instanciadas, invocamos el kernel, definiendo las dimensiones de la cuadrícula y los bloques como una tupla. La llamada al kernel incluye pasar nuestros datos, asignación de memoria para centroides y etiquetas, y un estado para la inicialización aleatoria de centroides.

def KMeansKernel(datos, centroidesSalida, etiquetasSalida)    fila = bloqueActual()    semilla = hiloActual()    filaCompartidaEntrada = arrayCompartido(shape=(tamañoLinea))    inerciaCompartida = arrayCompartido(shape=(numSemillas))    centroidesCompartidos = arrayCompartido(shape=(numSemillas, numClusters))    etiquetasCompartidas = arrayCompartido(shape=(numSemillas, tamañoLinea))    filaCompartidaEntrada = datos[fila]    sincronizarHilos()    si semilla == 0        centroides = inicializarCentroides(datos)    sincronizarHilos()        algoritmoKMeans(filaCompartidaEntrada, centroidesCompartidos, etiquetasCompartidas)        inerciaCompartida[semilla] = calcularInercia(filaCompartidaEntrada, centroidesCompartidos, etiquetasCompartidas)        sincronizarHilos()    si semilla == 0        indiceMinimaInercia = encontrarMinimaInercia(Inercia)    centroidesSalida = centroides[indiceMinimaInercia]    etiquetasSalida = etiquetas[indiceMinimaInercia]

Una vez en el dispositivo (también conocido como GPU), nuestro código ahora existe simultáneamente en toda la cuadrícula. Para saber dónde estamos en la cuadrícula, accedemos a los índices de bloque y hilo.

En la GPU, la memoria se predetermina a la memoria global, que se almacena en DRAM. La memoria compartida es privada para los hilos dentro del mismo bloque.

Al transferir la sola fila de datos (que todos los hilos en un bloque deben leer) a la memoria compartida, reducimos el tiempo de acceso a la memoria. En la función cuda_kmeans, creamos memoria compartida para almacenar la fila de datos, centroides, etiquetas y la medida de precisión de cada semilla, llamada inercia.

En nuestro programa, cada hilo se corresponde con una semilla, y todos los hilos en un bloque trabajan en la misma fila de datos. Para cada bloque, un hilo crea secuencialmente 32 semillas y agrega sus resultados en una única estructura de datos para los otros hilos en el bloque.

Cuando un paso subsiguiente en el algoritmo depende de que esta agregación se haya completado, los hilos deben ser sincronizados utilizando la función interna de CUDA syncthreads(). NB: hay que tener mucho cuidado al ubicar las llamadas a syncthreads(), ya que intentar sincronizar hilos cuando no todos han tenido la oportunidad de completar puede provocar bloqueos y que colapse todo el programa.

Nuestra función de kernel, que se describe en seudocódigo a continuación, se llama cuda_kmeans. Esta función se encarga de organizar el proceso descrito anteriormente, reservando espacio para los resultados de las 32 semillas y seleccionando la mejor semilla para la salida final en centroides y etiquetas.

def KMeansDispositivo(filaDatos, centroides, etiquetas)    semilla = hiloActual()    filaCentroides = centroides[semilla]    filaEtiquetas = etiquetas[semilla]     filaCentroides = ordenar(filaCentroides)    yardStick = calcularYardstick(centroidesOrdenados)     centroidesAntiguos = arrayLocal(shape=(numSemillas, numClusters))    para iteracion en rango(100):        si convergencia(centroidesAntiguos, filaCentroides)            romper        centroidesAntiguos = copiar(filaCentroides)        asignarEtiquetas(filaDatos, filaCentroides, filaEtiquetas)        actualizarCentroides(filaDatos, filaCentroides, filaEtiquetas)

Desde cuda_kmeans luego llamamos al algoritmo k-means real, pasando la nueva memoria compartida instanciada. En nuestro algoritmo k-means, comenzamos seleccionando los centroides iniciales y luego los ordenamos de menor a mayor. Asignamos iterativamente los puntos de datos al centroide más cercano y actualizamos las ubicaciones de los centroides hasta la convergencia.

Para poder determinar si se ha alcanzado la convergencia, utilizamos una función auxiliar llamada encontrar_yardstick. Esta función calcula y devuelve la distancia más pequeña entre dos centroides iniciales (yardstick). Nuestra condición de convergencia se cumple cuando ninguno de los centroides se ha movido más de yardstick veces epsilon durante una sola iteración.

Después de la convergencia, volvemos a cuda_kmeans. Aquí, determinamos la mejor semilla calculando la distancia euclidiana al cuadrado entre cada centroide y sus puntos de datos. La semilla con el agrupamiento más efectivo, indicado por la menor inercia, se considera la mejor. Luego tomamos los centroides y etiquetas de esta semilla y los copiamos en una sola fila en nuestros arrays de salida. Una vez que todos los bloques han terminado, estos resultados se copian de vuelta al host (CPU).

Una introducción a Numba

La forma más sencilla de diseñar núcleos personalizados es utilizando Numba. Numba es una biblioteca de Python que se puede utilizar para compilar código Python en núcleos de CUDA.

En el fondo, Numba se conecta con CUDA. Para paralelizar tu código, Numba compila tu código designado para la GPU en un núcleo y lo pasa a la GPU, dividiendo la lógica del programa en dos partes principales:

1. Código a nivel de CPU
2. Código a nivel de GPU

Utilizando Numba, separas y transfieres las partes secuenciales y paralelizables del código a la CPU y la GPU. Para compilar una función para la GPU, el programador utiliza el decorador @cuda.jit encima de la definición de la función, transformando así esta función en un núcleo que se llama desde la CPU (anfitrión) pero se ejecuta en paralelo en la GPU (dispositivo).

Python Numba

Enlace al Colab.

Numba sirve como puente entre el código Python y la plataforma CUDA. Debido a que el código Python es casi idéntico al pseudocódigo del algoritmo anterior, solo voy a proporcionar un par de ejemplos de sintaxis relevante clave.

cuda_kmeans[(NUM_ROWS,), (NUM_SEEDS,)](input_rows, output_labels, output_centroids, random_states)

Después de instanciar las variables globales y las estructuras de datos necesarias, podemos llamar al núcleo, cuda_kmeans, desde el anfitrión. Numba requiere dos tuplas para la dimensionalidad de los bloques y los hilos. Como utilizaremos bloques y hilos unidimensionales, el segundo índice de cada tupla está vacío. También pasamos nuestras estructuras de datos y una matriz de estados aleatorios para la instanciación de la semilla aleatoria.

@cuda.jit()def cuda_kmeans(input, output_labels, output_centroids, random_states):    row = cuda.blockIdx.x    seed = cuda.threadIdx.x    shared_input_row = cuda.shared.array(shape=(LINE_SIZE), dtype=np.float32)    shared_inertia = cuda.shared.array(shape=(NUM_SEEDS), dtype=np.float32)    shared_centroids = cuda.shared.array(shape=(NUM_SEEDS, NUM_CLUSTERS), dtype=np.float32)    shared_labels = cuda.shared.array(shape=(NUM_SEEDS, LINE_SIZE), dtype=np.int32)    if seed == 0:        get_initial_centroids(shared_input_row, shared_centroids, random_states)    cuda.syncthreads()    ...    kmeans(shared_input_row, shared_labels, shared_centroids)

Utilizamos el decorador Numba @cuda.jit() para marcar la función para la compilación de la GPU. Utilizando la notación cuda.blockIdx.x y cuda.threadIdx.x, obtenemos el índice actual del núcleo. Las matrices compartidas se instancian utilizando cuda.shared.array con dos argumentos, shape y el tipo, ambos deben ser conocidos en tiempo de compilación. Después de obtener los centroides y llenar la fila con datos, llamamos a la función kmeans, poblamos las matrices compartidas y hacemos una llamada a cuda.syncthreads().

@cuda.jit(device=True)def kmeans(data_row, output_labels, output_centroids):     seed = cuda.threadIdx.x    labels_row = output_labels[seed]    centroids_row = output_centroids[seed]        ...    old_centroids = cuda.local.array(shape=(NUM_CLUSTERS), dtype=np.float32)    for iteration in range(NUM_ITERATIONS):            if iteration > 0:                if converged(centroids_row, old_centroids, yard_stick * EPSILON_PERCENT, iteration):                    break      # Assign labels and update centroids

K-means es una función del dispositivo porque se llama desde la función del núcleo. Por lo tanto, debemos especificar device=True en el decorador: @cuda.jit(device=True). La función k-means luego obtiene la fila actual para las etiquetas y los centroides y se ejecuta hasta la convergencia.

Con solo una docena de líneas de código adicionales y un poco de esfuerzo, tu código de Python puede convertirse en un núcleo optimizado listo para su uso en paralelo.

Nuestro k-means paralelizado reduce drásticamente el tiempo de cálculo, sin embargo, envolver y compilar un lenguaje de alto nivel como Python no es necesariamente óptimo. Curioso por ver si escribir el código en C aceleraría nuestro proyecto, me sumergí en el mundo de CUDA C.

Una Introducción a C

En Python, la asignación de memoria y tipos es automática, mientras que en C, la memoria puede asignarse en la pila o en el montón, ambas requieren una declaración explícita de tipo y se les asigna una cantidad fija de memoria. La memoria de la pila se asigna y libera automáticamente por el compilador, mientras que la memoria del montón se asigna manualmente en tiempo de ejecución utilizando funciones como malloc(), siendo el programador responsable de su liberación.

Los punteros son herramientas que almacenan las direcciones de memoria de las variables. El tipo de datos al que se refiere el puntero se define durante su declaración. Los punteros se especifican utilizando un asterisco (*). Para obtener la dirección de una variable, llamada referencia, se utiliza un ampersand (&). Para acceder al valor de un puntero, se desreferencia nuevamente mediante un asterisco.

Un puntero doble, o un puntero a un puntero, almacena la dirección de un puntero. Esto es útil para modificar las direcciones de otros punteros al pasar matrices a funciones. Al pasar matrices a funciones, se pasan como punteros a su primer elemento, sin información de tamaño y confiando en la aritmética de punteros para poder indexar a través de la matriz. Para devolver una matriz desde una función, se pasaría la dirección de un puntero utilizando & y se recibiría con un puntero doble **, lo que te permite modificar la dirección del puntero original, pasando así la matriz.

int var = 100; // declaración de tipoint *ptr = &var; // uso de un puntero y referenciaint **double_ptr = &ptr; // ejemplo de puntero dobleprintf(“Desreferenciar double_ptr y ptr: %d %d \n:”, **double_ptr, *ptr)int *ptr = 100; // inicializar puntero de tipo int

CUDA C

Enlace al Colab.

CUDA es una plataforma informática que utiliza las GPUs de NVIDIA para paralelizar problemas computacionales complejos. Basándonos en tu recién adquirido dominio de C (bromeando), la estructura del código CUDA C es literalmente idéntica a la estructura del pseudocódigo que revisamos.

En el lado de la CPU, configuramos algunas constantes que le indican al algoritmo qué esperar, importamos librerías, inicializamos nuestras variables y definimos algunas macros.

#define NUM_ROWS 00000        // dimensión y de nuestro conjunto de datos, número de bloques#define LINE_SIZE 100         // dimensión x de nuestro conjunto de datos#define NUM_ITERATIONS 100    // número máximo de iteraciones#define NUM_CLUSTERS 3        // estamos ejecutando k = 3#define MAX_INPUT_VALUE 100   // límite superior en los datos#define NUM_SEEDS 32          // El número de semillas/hilos es 1/3 de LINE_SIZE#define EPSILON_PERCENT 0.02  // Condición de convergenciavoid initInputData(float** input) {    srand(1);     // asignar memoria para los datos        ... // inicializar los datos usando malloc y rand    // asignar memoria en la GPU    cudaMalloc(input, NUM_ROWS * LINE_SIZE * sizeof(float));     // copiar los datos desde la CPU (sample_data) a la memoria de la GPU    cudaMemcpy(*input, sample_data, NUM_ROWS * LINE_SIZE * sizeof(float), cudaMemcpyHostToDevice);    free(sample_data);}int main() {    float* inputData; // inicializar los datos de entrada, las dimensiones serán NUM_ROWS x LINE_SIZE    initInputData(&inputData); // desreferenciar y pasar a la función    // inicializar etiquetas de salida y centroides    cudaExtent labelsExtent = make_cudaExtent(sizeof(int), LINE_SIZE, NUM_ROWS);    cudaPitchedPtr outputLabels; // crear el puntero necesario para la siguiente llamada    cudaMalloc3D(&outputLabels, labelsExtent); // asignar memoria en la GPU        cudaExtent centroidsExtent = make_cudaExtent(sizeof(float), NUM_CLUSTERS, NUM_ROWS);    cudaPitchedPtr outputCentroids; // crear el puntero necesario para la siguiente llamada    cudaMalloc3D(&outputCentroids, centroidsExtent); // asignar memoria en la GPU    cuda_kmeans <<<NUM_ROWS, NUM_SEEDS>>> (inputData, outputLabels, outputCentroids);    cudaDeviceSynchronize();        ... // copiar la salida del dispositivo de nuevo a la CPU}

Descompongamos las diferencias.

La función principal comienza inicializando los datos creando un puntero y pasándolo a initInputData como una dirección a un puntero. La función recibe inputData como un puntero a un puntero (float** input), lo que permite que la función modifique la dirección mantenida por el puntero original. Luego, se apunta a una dirección de memoria GPU que se inicializa usando cudaMalloc y se llena utilizando cudaMemcpy, copiando datos del arreglo host temporal sample_data, ya poblado con números aleatorios.

Luego, el código asigna memoria en el dispositivo para almacenar los resultados de nuestra función k-means. La función utiliza make_cudaExtent para crear un objeto cudaExtent cuyo propósito es encapsular las dimensiones del arreglo multidimensional.

El tipo cudaPitchedPointer se utiliza para definir un puntero que puede direccionar este espacio de memoria pitch. Este tipo de puntero está específicamente diseñado para funcionar con memoria asignada por cudaMalloc3D, que toma tanto un cudaPitchedPtr como un objeto cudaExtent para asignar la memoria lineal en la GPU.

cuda_kmeans <<<NUM_ROWS, NUM_SEEDS>>> (inputData, outputLabels, outputCentroids);

Entrando en el código de la GPU, definimos la cuadrícula de manera que cada bloque corresponda a una fila de datos y cada hilo corresponda a una semilla.

Las semillas se inicializan por un solo hilo en cada bloque, asegurando semillas completamente diferentes.

__global__ void cuda_kmeans(float* input, cudaPitchedPtr outputLabels, cudaPitchedPtr outputCentroids) {    int row = blockIdx.x;    int seed = threadIdx.x;    // shared memory is shared between threads in blocks    __shared__ float input_shm[LINE_SIZE];    __shared__ float error_shm[NUM_SEEDS];    __shared__ float seed_centroids[NUM_SEEDS][NUM_CLUSTERS];    __shared__ int seed_labels[NUM_SEEDS][LINE_SIZE];        ... // get a single row of data    ... // populate input_shm    ... // populating the struct core_params    // the actual k-means function    kmeans(core_params);     // find seed with smallest error    calcError(core_params);    __syncthreads();    if (seed == 0) {        int* labels_line = LABELS_LINE(outputLabels, row);        float* centroids_line = CENTROIDS_LINE(outputCentroids, row);        labels_line[threadIdx.x] = seed_labels[seed][threadIdx.x];        centroids_line[threadIdx.x] = seed_centroids[seed][threadIdx.x];    }}

El código de CUDA C utiliza memoria compartida para los datos, centroides, etiquetas y errores. Sin embargo, a diferencia de Python, el código toma los punteros a memoria compartida y los almacena en una estructura, que es simplemente un método para pasar variables en masa. Finalmente, cuda_kmeans llama al algoritmo real k-means y pasa core_params.

__device__ void kmeans(core_params_t& core_params) {    DECLARE_CORE_PARAMS(core_params);    getInitialCentroids(core_params);    sort_centroids(centroids, num_clusters);    float yard_stick = findYardStick(core_params);    float* oldCentroids = (float*)malloc(NUM_CLUSTERS * sizeof(float));    struct work_params_t work_params;    work_params.min = find_min(line, LINE_SIZE);    work_params.max = find_max(line, LINE_SIZE);    work_params.aux_buf1 = (int*)malloc(NUM_CLUSTERS * sizeof(int));    work_params.aux_buf2 = (int*)malloc(NUM_CLUSTERS * sizeof(int));    work_params.aux_buf3 = (float*)malloc(NUM_CLUSTERS * sizeof(float));    for (int iterations = 0; true; iterations++) {        bool stop = (iterations > 100) || (iterations > 0 && (converged(core_params, oldCentroids, yard_stick * EPSILON_PERCENT)));        if (stop)            break;        memcpy(oldCentroids, core_params.centroids, NUM_CLUSTERS * sizeof(float));        getLabels(core_params);        getCentroids(core_params, work_params);    }    free(work_params.aux_buf1);    free(work_params.aux_buf2);    free(work_params.aux_buf3);    free(oldCentroids);}

En la función del dispositivo, antes que nada, extraemos los valores de la estructura core_params en variables usando la macro DECLARE_CORE_PARAMS.

Luego, ejecutamos el mismo algoritmo k-means que en Python, con la única diferencia de que pasamos la estructura en lugar de variables y tenemos que gestionar memoria, punteros y tipos.

Benchmarks

Para comparar nuestros algoritmos con k-means no paralelizado, importamos el módulo k-means de scikit-learn.

Para nuestro benchmark, ejecutamos 100,000 filas por 100 columnas con tres clusters. Debido a que scikit-learn no tiene un k-means paralelizado para diferentes filas, ejecutamos las filas secuencialmente en un bucle for.

Para nuestros benchmarks en colab, utilizamos la instancia gratuita T4 GPU Colab.

Los resultados son buenos: el código Python Numba es dos órdenes de magnitud más rápido que el código de CPU no paralelizado y el código CUDA C es tres órdenes de magnitud más rápido. Las funciones del kernel se escalan fácilmente y el algoritmo se puede modificar para admitir clústering en dimensiones más altas.

Tenga en cuenta que la generación de centroides iniciales aleatorios en ambos algoritmos, C y Python, no está completamente optimizada para utilizar todos los núcleos. Cuando se mejore, es posible que el algoritmo Python se acerque al tiempo de ejecución del código C.

Después de ejecutar la función k-means cien veces en diferentes conjuntos de datos y registrar los tiempos resultantes, observamos que la primera iteración es significativamente más lenta debido al tiempo que lleva compilar tanto C como Python en Colab.

Conclusión

Ahora deberías estar listo para escribir tus propios kernels de GPU personalizados. Una pregunta que queda es: ¿deberías usar CUDA C o Numba para procesamiento de datos paralelizado? Depende. Ambos son mucho más rápidos que el scikit-learn estándar. Aunque en mi caso, la implementación de k-means en lote mediante CUDA C resultó ser aproximadamente 3.5 veces más rápida que una implementación equivalente usando Numba, Python ofrece algunas ventajas importantes, como legibilidad y menos dependencia de habilidades especializadas en programación C en equipos que trabajan principalmente en Python. Además, el tiempo de ejecución de tu implementación específica dependerá de cuán bien hayas optimizado tu código para, por ejemplo, no desencadenar operaciones serializadas en la GPU. En conclusión, si eres nuevo tanto en C como en programación paralelizada, te sugiero que comiences con Numba para prototipar tu algoritmo y luego lo traduzcas a CUDA C si necesitas una mayor aceleración.

Referencias

1. Scikit-learn: Aprendizaje automático en Python, Pedregosa et al., JMLR 12, pp. 2825–2830, 2011.
2. NVIDIA, Vingelmann, P. & Fitzek, F.H.P., 2020. CUDA, versión: 10.2. 89, Disponible en: https://developer.nvidia.com/cuda-toolkit.
3. Lam, Siu Kwan, Antoine Pitrou y Stanley Seibert. “Numba: Un compilador de JIT basado en LLVM para Python”. Actas del Segundo Taller sobre Infraestructura del Compilador LLVM en HPC. 2015.
4. Harris, C.R., Millman, K.J., van der Walt, S.J. et al. “Programación de matrices con NumPy.” Nature 585, 357–362 (2020). DOI: 10.1038/s41586–020–2649–2. (Enlace del editor).

We will continue to update Zepes; if you have any questions or suggestions, please contact us!