Documents

Beyond CUDA: How DeepSeek Optimized at the PTX Level

41 min readFeb 5, 2025Feb 5, 2025

https://dev.to/datamonk_/how-deepseek-is-making-high-performance-ai-accessible-to-all-26fp

I don't know much about low-level internals, but

I came across an interesting article worth digging into.

DeepSeek has taken the LLM and AI world by storm, and its architecture is now something even non-specialists have heard of.

That said, it's fair to say DeepSeek's most impactful contributions don't actually lie in those areas.

The things I reviewed before aren't really the most important parts:

1. Attention architectures like MLA
2. MoE structure with shared experts
3. Self-evolution (a bold approach, but not the most "wow" part, in my view)

So what actually matters most?
The real question is: how did they achieve that level of efficiency on H800 GPUs?


Beyond CUDA

Working with AI ultimately means working with GPUs,
and the standard way to do that is through Python and PyTorch — not low-level assembly.

Even dropping down to the CUDA level is no small feat.

But the real question is whether these languages for GPU computation actually guarantee maximum efficiency.


How did DeepSeek construct that mysterious mixed-precision training architecture?

DeepSeek was always that kind of company, but —

whether driven by GPU export restrictions or not — they appear to have gone deep into the low-level stack in pursuit of maximum efficiency.

They moved beyond PyTorch to CUDA, and when that wasn't enough, they pushed further by customizing at the PTX level.


The GH100 die (144 SMs). The H100 SXM5 has 132 SMs; the PCIe variant has 114.

This is the H100 processor — the same architecture as the H800 DeepSeek used.

(What I know is that an SM, or Streaming Multiprocessor, is a kind of functional unit and the primary performance metric of a GPU.)

A single Streaming Multiprocessor looks like the above,
and the individual Tensor Cores inside it deliver the TeraFLOPS figures shown.

Looking at even a simple performance table,
you naturally start thinking about where the bottlenecks are
and how to distribute resources efficiently.

But when you're also setting up a distributed training environment and connecting multiple nodes together, keeping all of that in mind is…

incredibly hard — which is why most people just rely on PyTorch.

That's exactly where DeepSeek sets itself apart.

They open-sourced everything — model architecture, parameters, all of it — except for one thing.

Specifically: how they actually trained the model using something other than standard CUDA.


What does that customized PTX actually look like?

https://docs.nvidia.com/cuda/inline-ptx-assembly/contents.html

PTX is essentially a near-assembly-level language.

Following the stack from top to bottom:

CUDA → PTX (virtual assembly level) → GPU driver-level compiler → GPU

The GPU's own driver-level ISA is not publicly disclosed and is device-specific —
NVIDIA never publishes it, and there is no standardized ISA across devices.

PTX instructions, on the other hand, are device-independent.

DeepSeek targeted exactly this layer to push GPU optimization further.


The paper describes a range of techniques in this space, but…

Even CUDA-level programming is already low-level C-style code, so working at this layer is genuinely difficult.
(And then there's assembly-level code on top of that…)

  1. Used 2,048 NVIDIA H800 GPUs (8 GPUs per node).
    Applied various techniques to optimize NVLink and NVSwitch communication bandwidth.
    (Communication itself in distributed training was identified as a bottleneck.)
  2. Overlapped computation and communication.
  3. Implemented a dual pipeline.
  4. Dedicated 20 Streaming Multiprocessors exclusively to communication.
  5. Designed a custom mixed-precision training loop:
    • Multiply FP32 tensors by a scaling factor and convert to FP8
    • Apply GEMM (FP8 General Matrix Multiplication), Bias, and RMSNorm
    • Restore back to FP32
  6. And more…
// fp8_mixed_precision_training.cu
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cuda_runtime.h>
#include <cublas_v2.h>

// 오류 체크 매크로
#define CUDA_CHECK(err) do { \
    cudaError_t err_ = (err); \
    if (err_ != cudaSuccess) { \
        std::cerr << "CUDA Error: " << cudaGetErrorString(err_) \
                  << " at " << __FILE__ << ":" << __LINE__ << std::endl; \
        exit(EXIT_FAILURE); \
    } \
} while(0)

#define CUBLAS_CHECK(err) do { \
    cublasStatus_t err_ = (err); \
    if (err_ != CUBLAS_STATUS_SUCCESS) { \
        std::cerr << "CUBLAS Error at " << __FILE__ << ":" << __LINE__ << std::endl; \
        exit(EXIT_FAILURE); \
    } \
} while(0)

// FP8 타입 정의 (실제 FP8 연산은 최신 GPU에서 지원하지만, 여기서는 char로 시뮬레이션)
typedef char fp8_t;
const float FP8_MAX = 127.0f;  // FP8에서 표현 가능한 최대 절대값

// --------------------------------------------------------------------------
// [1] FP32 -> FP8 quantization / FP8 -> FP32 dequantization 커널
// --------------------------------------------------------------------------

// FP32 배열을 주어진 scaling factor(scale)를 곱해 FP8로 quantize (각 값은 [-127,127]로 clamp)
__global__ void quantizeKernel(const float* input, fp8_t* output, float scale, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if(idx < n) {
        float scaled = input[idx] * scale;  // scale = FP8_MAX / (max_abs)
        int q = (int)roundf(scaled);
        if(q > 127) q = 127;
        if(q < -127) q = -127;
        output[idx] = (fp8_t)q;
    }
}

// FP8 배열을 dequantize하여 FP32 배열로 변환 (invScale = 1/scale)
__global__ void dequantizeKernel(const fp8_t* input, float* output, float invScale, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if(idx < n) {
        int q = (int)input[idx];
        output[idx] = ((float)q) * invScale;
    }
}

// --------------------------------------------------------------------------
// [2] RMSNorm Forward 커널 (y = x / sqrt(mean(x^2)+epsilon))
//    – 실제로는 저장하지 않고, 역전파 시 재계산(recomputation)을 수행하여 메모리 절약 효과를 얻습니다.
//    (여기서는 데모용으로 순전파에서도 RMSNorm 결과를 출력합니다.)
// --------------------------------------------------------------------------
__global__ void rmsNormForwardKernel(const float* x, float* y, int n, float epsilon) {
    // 각 블록이 한 샘플을 처리한다고 가정하고, 블록 내 쓰레드가 협력해 reduction 수행 (간단히 구현)
    extern __shared__ float sdata[];
    int tid = threadIdx.x;
    int idx = blockIdx.x * n + tid;
    float val = (tid < n) ? x[idx] : 0.0f;
    sdata[tid] = val * val;
    __syncthreads();
    // reduction: (단순 구현 – n이 쓰레드 수와 같다고 가정)
    for (int s = n/2; s > 0; s >>= 1) {
        if(tid < s && tid + s < n) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
    float rms = sqrtf(sdata[0] / n + epsilon);
    if(tid < n) {
        y[idx] = x[idx] / rms;
    }
}

// --------------------------------------------------------------------------
// [3] FP8 GEMM 함수
//    – 입력 A와 가중치 B는 FP8 형식으로 저장되어 있으며, 각각 scaling factor scaleA, scaleB를 갖습니다.
//    – 내부에서는 임시 FP32 버퍼로 dequantize한 후, cuBLAS의 Sgemm을 통해 FP32 누적으로 GEMM을 수행합니다.
//    – (참고: 실제 Nvidia FTX 구현에서는 Tensor Core와 FP8 전용 커널을 사용하며, 일정 간격마다 CUDA Core로 promotion을 수행합니다.)
// --------------------------------------------------------------------------
void fp8GEMM(cublasHandle_t handle, int M, int N, int K,
             const fp8_t* d_A, float scaleA,
             const fp8_t* d_B, float scaleB,
             float* d_C) {  // 결과는 FP32
    float *d_A_fp32, *d_B_fp32;
    CUDA_CHECK(cudaMalloc(&d_A_fp32, M * K * sizeof(float)));
    CUDA_CHECK(cudaMalloc(&d_B_fp32, K * N * sizeof(float)));

    int total_A = M * K;
    int total_B = K * N;
    int blockSize = 256;
    int numBlocksA = (total_A + blockSize - 1) / blockSize;
    int numBlocksB = (total_B + blockSize - 1) / blockSize;
    float invScaleA = 1.0f / scaleA;
    float invScaleB = 1.0f / scaleB;
    dequantizeKernel<<<numBlocksA, blockSize>>>(d_A, d_A_fp32, invScaleA, total_A);
    CUDA_CHECK(cudaGetLastError());
    dequantizeKernel<<<numBlocksB, blockSize>>>(d_B, d_B_fp32, invScaleB, total_B);
    CUDA_CHECK(cudaGetLastError());

    float alpha = 1.0f, beta = 0.0f;
    // cuBLAS는 기본적으로 열우선 저장을 가정하므로, 여기서는 단순화를 위해 행렬 크기를 그대로 사용합니다.
    CUBLAS_CHECK(cublasSgemm(handle,
                             CUBLAS_OP_N, CUBLAS_OP_N,
                             N, M, K,
                             &alpha,
                             d_B_fp32, N,
                             d_A_fp32, K,
                             &beta,
                             d_C, N));
    cudaFree(d_A_fp32);
    cudaFree(d_B_fp32);
}

// --------------------------------------------------------------------------
// [4] Simple FP32 matrix bias-add kernel
// --------------------------------------------------------------------------
__global__ void addBiasKernel_FP32(float* mat, const float* bias, int cols, int total) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < total) {
        int col = idx % cols;
        mat[idx] += bias[col];
    }
}

// --------------------------------------------------------------------------
// [5] MSE loss kernel (simple demo)
// --------------------------------------------------------------------------
__global__ void mseLossKernel_FP32(const float* pred, const float* target, float* loss, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if(idx < n) {
        float diff = pred[idx] - target[idx];
        atomicAdd(loss, diff * diff);
    }
}

// --------------------------------------------------------------------------
// [6] CPU-side EMA (Exponential Moving Average) update (synchronous)
// --------------------------------------------------------------------------
void updateEMA(const float* param, float* ema_param, int size, float decay) {
    for (int i = 0; i < size; i++) {
        ema_param[i] = decay * ema_param[i] + (1.0f - decay) * param[i];
    }
}

// --------------------------------------------------------------------------
// [7] Main training code
// --------------------------------------------------------------------------
int main() {
    // Hyperparameters and size configuration
    const int batch_size   = 64;
    const int input_dim    = 512;  // e.g., embedding dimension
    const int output_dim   = 256;  // e.g., output dimension (output head)
    const float learning_rate = 0.001f;
    const int epochs       = 5;
    const int num_batches  = 10;

    // Create cuBLAS handle
    cublasHandle_t handle;
    CUBLAS_CHECK(cublasCreate(&handle));

    // ----------------------------------------------------------------------
    // [A] Host memory: inputs and targets (same data reused each batch for demo)
    // ----------------------------------------------------------------------
    float *h_input  = new float[batch_size * input_dim];
    float *h_target = new float[batch_size * output_dim];
    for (int i = 0; i < batch_size * input_dim; i++) {
        h_input[i] = static_cast<float>(rand()) / RAND_MAX;
    }
    for (int i = 0; i < batch_size * output_dim; i++) {
        h_target[i] = static_cast<float>(rand()) / RAND_MAX;
    }

    // ----------------------------------------------------------------------
    // [B] Device memory allocation: input (FP32)
    // ----------------------------------------------------------------------
    float *d_input;
    CUDA_CHECK(cudaMalloc(&d_input, batch_size * input_dim * sizeof(float)));
    CUDA_CHECK(cudaMemcpy(d_input, h_input, batch_size * input_dim * sizeof(float), cudaMemcpyHostToDevice));

    // ----------------------------------------------------------------------
    // [C] Store linear layer weights in FP8 for intermediate-layer GEMMs.
    //      (Embedding and output head can be kept in FP32; here only bias is FP32 for simplicity.)
    // ----------------------------------------------------------------------
    fp8_t *d_W_fp8;
    CUDA_CHECK(cudaMalloc(&d_W_fp8, input_dim * output_dim * sizeof(fp8_t)));
    // Initialize FP32 master weights on the host, then quantize to FP8
    float *h_W = new float[input_dim * output_dim];
    for (int i = 0; i < input_dim * output_dim; i++) {
        h_W[i] = ((float)rand() / RAND_MAX - 0.5f) * 0.1f;
    }
    // Compute a global scaling factor (production code would use per-tile scaling)
    float max_val = 0.0f;
    for (int i = 0; i < input_dim * output_dim; i++) {
        float abs_val = fabsf(h_W[i]);
        if (abs_val > max_val) max_val = abs_val;
    }
    float scale_W = FP8_MAX / max_val;  // weight quantization scale
    // FP8 quantization (on host)
    fp8_t *h_W_fp8 = new fp8_t[input_dim * output_dim];
    for (int i = 0; i < input_dim * output_dim; i++) {
        int q = (int)roundf(h_W[i] * scale_W);
        if(q > 127) q = 127;
        if(q < -127) q = -127;
        h_W_fp8[i] = (fp8_t)q;
    }
    CUDA_CHECK(cudaMemcpy(d_W_fp8, h_W_fp8, input_dim * output_dim * sizeof(fp8_t), cudaMemcpyHostToDevice));

    // ----------------------------------------------------------------------
    // [D] Bias allocation (FP32)
    // ----------------------------------------------------------------------
    float *d_bias;
    CUDA_CHECK(cudaMalloc(&d_bias, output_dim * sizeof(float)));
    float *h_bias = new float[output_dim];
    for (int i = 0; i < output_dim; i++) {
        h_bias[i] = 0.0f;
    }
    CUDA_CHECK(cudaMemcpy(d_bias, h_bias, output_dim * sizeof(float), cudaMemcpyHostToDevice));

    // ----------------------------------------------------------------------
    // [E] EMA parameters (CPU): stores EMA values for master weights and bias
    // ----------------------------------------------------------------------
    float *ema_W = new float[input_dim * output_dim];
    float *ema_bias = new float[output_dim];
    // Initialize EMA from the current FP32 master weights and bias (weights are dequantized)
    for (int i = 0; i < input_dim * output_dim; i++) {
        ema_W[i] = ((float)h_W_fp8[i]) / scale_W;
    }
    for (int i = 0; i < output_dim; i++) {
        ema_bias[i] = h_bias[i];
    }
    float ema_decay = 0.999f;

    // ----------------------------------------------------------------------
    // [F] Training loop
    // ----------------------------------------------------------------------
    for (int epoch = 0; epoch < epochs; epoch++) {
        float epoch_loss = 0.0f;
        for (int batch = 0; batch < num_batches; batch++) {
            // (Same input reused each batch for demo purposes)
            // Input is already on device in FP32; proceed with FP8 quantization.

            // [F-1] Input quantization: FP32 -> FP8
            fp8_t *d_input_fp8;
            CUDA_CHECK(cudaMalloc(&d_input_fp8, batch_size * input_dim * sizeof(fp8_t)));
            // Compute global scaling factor for the input on the host (production code uses tile-wise scaling)
            float max_input = 0.0f;
            for (int i = 0; i < batch_size * input_dim; i++) {
                float v = fabsf(h_input[i]);
                if(v > max_input) max_input = v;
            }
            float scale_input = FP8_MAX / max_input;
            int total_input = batch_size * input_dim;
            int blockSize = 256;
            int numBlocks = (total_input + blockSize - 1) / blockSize;
            quantizeKernel<<<numBlocks, blockSize>>>(d_input, d_input_fp8, scale_input, total_input);
            CUDA_CHECK(cudaGetLastError());

            // [F-2] FP8 GEMM: d_hidden = d_input_fp8 (batch_size x input_dim) * d_W_fp8 (input_dim x output_dim)
            // Output d_hidden is FP32 (batch_size x output_dim)
            float *d_hidden;
            CUDA_CHECK(cudaMalloc(&d_hidden, batch_size * output_dim * sizeof(float)));
            fp8GEMM(handle, batch_size, output_dim, input_dim, d_input_fp8, scale_input, d_W_fp8, scale_W, d_hidden);

            // [F-3] Add bias: d_hidden += bias
            int total_hidden = batch_size * output_dim;
            numBlocks = (total_hidden + blockSize - 1) / blockSize;
            addBiasKernel_FP32<<<numBlocks, blockSize>>>(d_hidden, d_bias, output_dim, total_hidden);
            CUDA_CHECK(cudaGetLastError());

            // [F-4] RMSNorm (to save memory, the forward pass does not store intermediate results;
            // they are recomputed during backprop — for this demo we store the output in d_rmsnorm_out)
            float *d_rmsnorm_out;
            CUDA_CHECK(cudaMalloc(&d_rmsnorm_out, batch_size * output_dim * sizeof(float)));
            // One block per sample, blockDim = output_dim, with dynamic shared memory
            for (int i = 0; i < batch_size; i++) {
                rmsNormForwardKernel<<<1, output_dim, output_dim * sizeof(float)>>>(d_hidden + i * output_dim,
                                                                                    d_rmsnorm_out + i * output_dim,
                                                                                    output_dim, 1e-5f);
            }
            CUDA_CHECK(cudaDeviceSynchronize());

The example above was written at my request, since I don't know C/C# well enough to write CUDA code myself.

Even so, it's clear that even without PTX ISA — purely at the CUDA level — virtually everything has to be defined and allocated manually.

(If this had been written in PyTorch, you'd define a few layers, specify the data types, and be done with it...)

DeepSeek appears to have pushed exactly this kind of low-level optimization to its absolute limit (though the specifics haven't been fully disclosed).

![](/uploads/images/20260221_bc9ef5252d72bd51.png)
*Mirae Asset's Digital Research AI Weekly #45 — assessment of DeepSeek*

I would have liked to go through the papers in full and examine what the team was thinking at each individual stage, but...

That's something to attempt when the opportunity arises.

---

![](/uploads/images/20260221_de8c85f644876455.png)
*This is probably not the end of the CUDA moat.*

This topic is already getting a lot of coverage.

Some people have even argued that Nvidia's stock dropped because someone suggested moving beyond CUDA to PTX — even though both are ultimately just languages for programming Nvidia GPUs...

DeepSeek almost certainly didn't write everything at the PTX level either.

My guess is they started from a CUDA codebase and applied PTX ISA selectively to achieve a handful of targeted optimizations.
Tags
deepseek