C++ for AI / Machine Learning

Why C++ in AI/ML

Python is the training language; C++ is the deployment language. Most ML frameworks have C++ backends (PyTorch/libtorch, TensorFlow, ONNX Runtime) and expose native C++ APIs. Use C++ when you need:

Low-latency inference — no Python GIL, no interpreter overhead
Embedded / edge deployment — no Python runtime on device
Custom operators / CUDA kernels — extend frameworks in C++
Production serving — integrate ML directly into your C++ service

ONNX Runtime — universal inference engine

ONNX Runtime runs models exported from PyTorch, TensorFlow, scikit-learn, and more.

cpp

Godbolt

// CMakeLists.txt:
// find_package(onnxruntime REQUIRED)
// target_link_libraries(myapp onnxruntime::onnxruntime)

#include <onnxruntime_cxx_api.h>
#include <vector>

int main() {
    Ort::Env env(ORT_LOGGING_LEVEL_WARNING, "inference");
    Ort::SessionOptions opts;
    opts.SetIntraOpNumThreads(4);
    // opts.AppendExecutionProvider_CUDA(0);  // enable GPU

    Ort::Session session(env, L"model.onnx", opts);

    // Input tensor setup
    std::vector<float> input_data(1 * 3 * 224 * 224, 0.5f);  // batch=1, CHW
    std::vector<int64_t> input_shape = {1, 3, 224, 224};

    Ort::MemoryInfo mem_info = Ort::MemoryInfo::CreateCpu(
        OrtArenaAllocator, OrtMemTypeDefault);

    Ort::Value input_tensor = Ort::Value::CreateTensor<float>(
        mem_info, input_data.data(), input_data.size(),
        input_shape.data(), input_shape.size());

    // Run inference
    const char* input_names[] = {"input"};
    const char* output_names[] = {"output"};

    auto outputs = session.Run(
        Ort::RunOptions{nullptr},
        input_names, &input_tensor, 1,
        output_names, 1);

    // Read output
    float* out_data = outputs[0].GetTensorMutableData<float>();
    auto out_shape = outputs[0].GetTensorTypeAndShapeInfo().GetShape();
    // out_data[0..N-1] = logits
}

libtorch — PyTorch C++ API

cpp

Godbolt

#include <torch/torch.h>
#include <torch/script.h>  // TorchScript

// Load a TorchScript model (export with torch.jit.script or torch.jit.trace)
torch::jit::script::Module model = torch::jit::load("model.pt");
model.eval();

// Create input tensor
auto input = torch::randn({1, 3, 224, 224});
if (torch::cuda::is_available()) input = input.to(torch::kCUDA);

// Run inference
std::vector<torch::jit::IValue> inputs = {input};
auto output = model.forward(inputs).toTensor();

// Postprocess
auto probs = torch::softmax(output, 1);
auto [top_prob, top_idx] = probs.topk(5, 1);

Defining a model in C++ (nn::Module)

cpp

Godbolt

struct Net : torch::nn::Module {
    torch::nn::Linear fc1{nullptr}, fc2{nullptr};

    Net(int in_features, int hidden, int num_classes) {
        fc1 = register_module("fc1", torch::nn::Linear(in_features, hidden));
        fc2 = register_module("fc2", torch::nn::Linear(hidden, num_classes));
    }

    torch::Tensor forward(torch::Tensor x) {
        x = torch::relu(fc1->forward(x));
        x = fc2->forward(x);
        return x;
    }
};

Net net(784, 256, 10);
auto optimizer = torch::optim::Adam(net.parameters(), 1e-3);

// Training step
optimizer.zero_grad();
auto loss = torch::nn::functional::cross_entropy(net.forward(batch_x), batch_y);
loss.backward();
optimizer.step();

BLAS/LAPACK for linear algebra

For inference without a full ML framework — matrix multiplications power transformers, CNNs, MLPs.

cpp

Godbolt

// OpenBLAS / MKL — dgemm for double, sgemm for float

// C = alpha*A*B + beta*C
// A: m×k, B: k×n, C: m×n (column-major convention)
#include <cblas.h>

void matmul(const float* A, const float* B, float* C, int m, int k, int n) {
    cblas_sgemm(
        CblasRowMajor,   // row-major storage
        CblasNoTrans, CblasNoTrans,
        m, n, k,
        1.0f,            // alpha
        A, k,
        B, n,
        0.0f,            // beta
        C, n
    );
}
// ~10-100× faster than naive loops on large matrices (SIMD + threading)

Custom CUDA kernel

cpp

Godbolt

// softmax kernel — simple illustrative version
__global__ void softmax_kernel(const float* input, float* output, int n) {
    int tid = threadIdx.x;
    extern __shared__ float smem[];

    // Load and find max (numerically stable softmax)
    smem[tid] = (tid < n) ? input[tid] : -FLT_MAX;
    __syncthreads();

    for (int stride = blockDim.x / 2; stride > 0; stride /= 2) {
        if (tid < stride) smem[tid] = fmaxf(smem[tid], smem[tid + stride]);
        __syncthreads();
    }
    float max_val = smem[0];

    smem[tid] = (tid < n) ? expf(input[tid] - max_val) : 0.0f;
    __syncthreads();

    for (int stride = blockDim.x / 2; stride > 0; stride /= 2) {
        if (tid < stride) smem[tid] += smem[tid + stride];
        __syncthreads();
    }
    float sum = smem[0];

    if (tid < n) output[tid] = expf(input[tid] - max_val) / sum;
}

// Launch
int n = 1024;
softmax_kernel<<<1, n, n * sizeof(float)>>>(d_input, d_output, n);

Embedding a tokenizer with HuggingFace tokenizers C++ bindings

cpp

Godbolt

// tokenizers-cpp (HuggingFace tokenizers Rust → C++ via FFI)
#include "tokenizers_cpp.h"

auto tokenizer = tokenizers::Tokenizer::FromFile("tokenizer.json");

std::string text = "The quick brown fox";
auto encoding = tokenizer->Encode(text, /*add_special_tokens=*/true);

std::vector<int> input_ids(encoding.GetIds().begin(), encoding.GetIds().end());
// input_ids ready for transformer model

On-device inference: NCNN (Tencent, mobile/edge)

cpp

Godbolt

#include "ncnn/net.h"

ncnn::Net net;
net.opt.use_vulkan_compute = true;  // GPU on Android/iOS
net.load_param("mobilenet.param");
net.load_model("mobilenet.bin");

ncnn::Mat in = ncnn::Mat::from_pixels_resize(rgb_data, ncnn::Mat::PIXEL_RGB,
    orig_w, orig_h, 224, 224);

// Normalize
const float mean_vals[3] = {103.94f, 116.78f, 123.68f};
const float norm_vals[3] = {0.017f, 0.017f, 0.017f};
in.substract_mean_normalize(mean_vals, norm_vals);

ncnn::Extractor ex = net.create_extractor();
ex.input("data", in);

ncnn::Mat out;
ex.extract("prob", out);
// out[i] = probability of class i

Key libraries

Library	Purpose	Notes
ONNX Runtime	Universal inference (CPU/GPU/TensorRT)	Microsoft, production-grade
libtorch	PyTorch C++ API	Training + inference
TensorFlow C API	TF inference	Official C API, C++ wrapper available
NCNN	Mobile/edge inference	No dependencies, Vulkan support
MNN	Mobile neural network	Alibaba, fast mobile inference
OpenBLAS / MKL	BLAS operations	Core of all matrix computation
Eigen	Matrix/vector math	Header-only, no BLAS needed
FAISS	Vector similarity search	Facebook AI, used for embeddings
Triton Inference Server	gRPC/HTTP serving	NVIDIA, ONNX/TensorRT/TorchScript

SIMD for inference hot paths

cpp

Godbolt

#include <immintrin.h>

// AVX2 dot product — 8 floats per instruction
float dot_avx2(const float* a, const float* b, int n) {
    __m256 sum = _mm256_setzero_ps();
    int i = 0;
    for (; i + 8 <= n; i += 8) {
        __m256 va = _mm256_loadu_ps(a + i);
        __m256 vb = _mm256_loadu_ps(b + i);
        sum = _mm256_fmadd_ps(va, vb, sum);  // FMA: sum += va * vb
    }
    // Horizontal sum
    float tmp[8];
    _mm256_storeu_ps(tmp, sum);
    float result = tmp[0]+tmp[1]+tmp[2]+tmp[3]+tmp[4]+tmp[5]+tmp[6]+tmp[7];
    for (; i < n; ++i) result += a[i] * b[i];  // remainder
    return result;
}

Quantization (INT8/FP16)

cpp

Godbolt

// ONNX Runtime INT8 quantization at load time
Ort::SessionOptions opts;
// ORT_ENABLE_ALL applies all quantization passes
opts.SetGraphOptimizationLevel(ORT_ENABLE_ALL);

// Or load a pre-quantized INT8 model (quantized offline with onnxruntime tools)
// python -m onnxruntime.quantization.quantize_static model.onnx model_int8.onnx calibration_data

// libtorch FP16
auto model_half = model.to(torch::kHalf);
auto input_half = input.to(torch::kHalf);
auto output = model_half.forward({input_half}).toTensor();

ONNX Runtime execution provider selection

Choosing the right execution provider unlocks 10–100× latency improvements over the CPU default. ONNX Runtime dispatches ops to the first provider that supports them, falling back to the next in the list.

cpp

Godbolt

#include <onnxruntime_cxx_api.h>

Ort::Env env(ORT_LOGGING_LEVEL_WARNING, "myapp");
Ort::SessionOptions opts;
opts.SetIntraOpNumThreads(1);                 // latency-optimised: 1 thread
opts.SetInterOpNumThreads(1);
opts.SetGraphOptimizationLevel(ORT_ENABLE_ALL);

// --- CPU only (default, always available) ---
// No extra call needed; CPU provider is implicit.

// --- CUDA (requires onnxruntime-gpu package) ---
OrtCUDAProviderOptions cuda_opts{};
cuda_opts.device_id = 0;
cuda_opts.cudnn_conv_algo_search = OrtCudnnConvAlgoSearchDefault;
cuda_opts.gpu_mem_limit = 2ULL * 1024 * 1024 * 1024; // 2 GB
opts.AppendExecutionProvider_CUDA(cuda_opts);

// --- TensorRT (highest throughput for NVIDIA GPUs) ---
OrtTensorRTProviderOptions trt_opts{};
trt_opts.device_id = 0;
trt_opts.trt_max_workspace_size = 1ULL << 30;   // 1 GB
trt_opts.trt_fp16_enable = 1;                   // FP16 kernels
trt_opts.trt_engine_cache_enable = 1;           // cache compiled engines
trt_opts.trt_engine_cache_path = "./trt_cache";
opts.AppendExecutionProvider_TensorRT(trt_opts);

Ort::Session session(env, L"model.onnx", opts);

Provider	Typical speedup vs CPU	Requires
CPU	1× (baseline)	nothing
CUDA	5–30×	CUDA toolkit, cuDNN
TensorRT	10–100×	TensorRT SDK, NVIDIA GPU
DirectML	5–20× (Windows)	DirectX 12
CoreML	5–15× (macOS/iOS)	Apple Silicon or ANE
OpenVINO	3–10×	Intel CPU/iGPU

Quantization and precision (FP32 → FP16 → INT8)

Reducing numerical precision shrinks model size and raises throughput at the cost of a small accuracy drop. Always validate on a held-out dataset after quantization.

bash

# Offline static INT8 quantization with a calibration dataset
python -m onnxruntime.quantization.quantize_static \
    --model_input  model_fp32.onnx \
    --model_output model_int8.onnx \
    --calibration_data_reader calib_reader.py \
    --quant_format QDQ \
    --activation_type QInt8 \
    --weight_type QInt8

cpp

Godbolt

// Runtime: load the pre-quantized model — no code change needed
Ort::Session session(env, L"model_int8.onnx", opts);

// libtorch: dynamic quantization (weights only, no calibration data)
auto quantized = torch::quantization::quantize_dynamic(
    model,
    {at::kLinear},     // quantize Linear layers
    torch::kQInt8
);

Precision	Size vs FP32	Typical accuracy drop	Speed (GPU)
FP32	1×	0 %	1×
FP16	0.5×	< 0.1 %	1.5–2×
INT8	0.25×	0.1–1 %	2–4×
INT4	0.125×	1–3 %	3–6×

Registering a custom C++ operator in LibTorch

When a built-in op is unavailable or too slow, you can register a C++ kernel that TorchScript calls by name. The kernel participates in autograd if you supply a backward function.

cpp

Godbolt

#include <torch/extension.h>

// Forward kernel
torch::Tensor my_relu_forward(const torch::Tensor& input) {
    return input.clamp_min(0);
}

// Backward kernel
torch::Tensor my_relu_backward(
    const torch::Tensor& grad_output,
    const torch::Tensor& input)
{
    return grad_output * (input > 0).to(grad_output.dtype());
}

// Register with the dispatcher
TORCH_LIBRARY(my_ops, m) {
    m.def("my_relu(Tensor input) -> Tensor");
}

TORCH_LIBRARY_IMPL(my_ops, CPU, m) {
    m.impl("my_relu", &my_relu_forward);
}

// Python-side (for hybrid workflows):
// import torch; torch.ops.load_library("libmy_ops.so")
// out = torch.ops.my_ops.my_relu(tensor)

For a CUDA kernel, add a TORCH_LIBRARY_IMPL(my_ops, CUDA, m) block that dispatches to your __global__ function.

Matrix library comparison (Eigen vs xtensor vs Armadillo)

Eigen is the default choice; xtensor suits NumPy-style broadcasting; Armadillo targets MATLAB users. All three are header-friendly and integrate cleanly into CMake projects.

cpp

Godbolt

// --- Eigen ---
#include <Eigen/Dense>
Eigen::MatrixXf A = Eigen::MatrixXf::Random(512, 512);
Eigen::MatrixXf B = Eigen::MatrixXf::Random(512, 512);
Eigen::MatrixXf C = A * B;                  // BLAS-backed on large matrices

// --- xtensor ---
#include <xtensor/xtensor.hpp>
#include <xtensor-blas/xlinalg.hpp>
xt::xtensor<float,2> X = xt::random::randn<float>({512, 512});
auto result = xt::linalg::dot(X, X);        // NumPy-like API

// --- Armadillo ---
#include <armadillo>
arma::fmat M = arma::randu<arma::fmat>(512, 512);
arma::fmat N = M * M;                       // LAPACK-backed

Library	API style	Header-only	Auto-BLAS	Sparse	Notes
Eigen	C++ expressions	Yes	Optional (MKL)	Yes	Fastest, widest adoption
xtensor	NumPy-like	Yes	Via xblas	No	Best for N-d arrays
Armadillo	MATLAB-like	No (link arma)	Yes (OpenBLAS)	Yes	Easiest for linear algebra
Blaze	Expression tmpl	Yes	Yes	Yes	High raw throughput

GPU memory management with CUDA

Picking the right memory type is the single biggest factor in CUDA kernel throughput. Pinned memory doubles PCIe transfer bandwidth; unified memory simplifies code at some performance cost.

cpp

Godbolt

#include <cuda_runtime.h>

// 1. Device memory — fastest for GPU-only access
float* d_buf;
cudaMalloc(&d_buf, N * sizeof(float));
cudaFree(d_buf);

// 2. Pinned (page-locked) host memory — enables async H→D transfers
float* h_pinned;
cudaMallocHost(&h_pinned, N * sizeof(float));   // or cudaHostAlloc
// Fill h_pinned, then async copy
cudaMemcpyAsync(d_buf, h_pinned, N * sizeof(float),
                cudaMemcpyHostToDevice, stream);
cudaFreeHost(h_pinned);

// 3. Zero-copy (mapped) memory — GPU reads directly from host RAM
//    Good for small buffers read once; avoids explicit copy
float* h_mapped;
cudaHostAlloc(&h_mapped, N * sizeof(float), cudaHostAllocMapped);
float* d_mapped;
cudaHostGetDevicePointer(&d_mapped, h_mapped, 0);
// Pass d_mapped to kernel — it reads from host DRAM over PCIe

// 4. Unified / managed memory — single pointer, driver migrates pages
float* um_buf;
cudaMallocManaged(&um_buf, N * sizeof(float));
// Accessible on both host and device; prefetch for performance:
cudaMemPrefetchAsync(um_buf, N * sizeof(float), device_id, stream);
cudaFree(um_buf);

Type	Allocation	Bandwidth	Use case
Device	cudaMalloc	Peak GPU	All persistent GPU data
Pinned	cudaMallocHost	2× H↔D	Frequent H↔D transfers
Zero-copy	cudaHostAllocMapped	PCIe limited	Occasional reads, small
Unified	cudaMallocManaged	Variable	Prototype / correctness

Batching strategies for inference

Batching amortises kernel launch overhead and saturates GPU compute units. The optimal batch size is the largest that fits in GPU memory without increasing per-sample latency beyond your SLA.

cpp

Godbolt

// Dynamic batching with a request queue (thread-safe accumulator)
#include <deque>
#include <mutex>
#include <condition_variable>

struct InferRequest { std::vector<float> input; std::promise<float> result; };

class BatchingEngine {
    std::deque<InferRequest*> queue_;
    std::mutex mu_;
    std::condition_variable cv_;
    int max_batch_ = 32;
    std::chrono::microseconds max_wait_{500};

public:
    void push(InferRequest* req) {
        std::lock_guard lk{mu_};
        queue_.push_back(req);
        cv_.notify_one();
    }

    void worker_loop(Ort::Session& session) {
        while (true) {
            std::unique_lock lk{mu_};
            cv_.wait_for(lk, max_wait_,
                [&]{ return queue_.size() >= (size_t)max_batch_; });

            int bs = std::min((int)queue_.size(), max_batch_);
            if (bs == 0) continue;

            std::vector<InferRequest*> batch(
                queue_.begin(), queue_.begin() + bs);
            queue_.erase(queue_.begin(), queue_.begin() + bs);
            lk.unlock();

            // Build batched tensor [bs, feature_dim], run session, scatter results
            run_batch(session, batch);
        }
    }
};

Rules of thumb for optimal batch size selection:

Start at batch 1 and double until latency exceeds SLA or GPU memory runs out.
Profile with nsys profile to see if the GPU is compute-bound (use larger batch) or memory-bound (already optimal).
For transformer models: sequence length × batch size determines activation memory; prefer shorter sequences over larger batches.

Deploying C++ inference servers

Triton Inference Server (NVIDIA) is the fastest path to production; a custom gRPC server gives full control. Both can serve ONNX, TensorRT, and TorchScript models.

bash

# Triton — model repository layout
models/
  my_model/
    config.pbtxt
    1/
      model.onnx       # or model.plan (TensorRT)

yaml

# config.pbtxt (Protocol Buffer text format)
name: "my_model"
backend: "onnxruntime"
max_batch_size: 64

input  [{ name: "input",  data_type: TYPE_FP32, dims: [3, 224, 224] }]
output [{ name: "output", data_type: TYPE_FP32, dims: [1000]        }]

dynamic_batching {
  preferred_batch_size: [8, 16, 32]
  max_queue_delay_microseconds: 500
}
instance_group [{ count: 2, kind: KIND_GPU, gpus: [0] }]

bash

# Launch
docker run --gpus all -p 8000:8000 -p 8001:8001 \
  -v $(pwd)/models:/models \
  nvcr.io/nvidia/tritonserver:24.01-py3 \
  tritonserver --model-repository=/models

# gRPC client
tritonclient.grpc.InferenceServerClient("localhost:8001")

For a custom gRPC server: use grpc++ with a .proto defining your Infer RPC, then call session.Run() inside the handler. Keep one Ort::Session per server process and protect it with a thread pool.

Profiling AI workloads

Profile before optimising — the bottleneck is rarely where you expect it. For CUDA workloads, Nsight Systems gives a timeline of CPU/GPU overlap; for CPU, perf and FLOP counting reveal compute efficiency.

bash

# Nsight Systems — full CPU+GPU timeline
nsys profile --trace=cuda,nvtx,osrt \
     --output=report \
     ./my_inference_app

nsys stats report.nsys-rep   # summary tables

# Nsight Compute — per-kernel roofline analysis
ncu --set full --target-processes all \
    --output kernel_report \
    ./my_inference_app

# Linux perf — CPU-side bottlenecks
perf stat -e cycles,instructions,cache-misses,fp_arith_inst_retired.256b_packed_single \
    ./my_inference_app

cpp

Godbolt

// NVTX ranges — annotate your own ranges in the Nsight timeline
#include <nvtx3/nvToolsExt.h>

nvtxRangePushA("preprocessing");
preprocess(input);
nvtxRangePop();

nvtxRangePushA("inference");
session.Run(...);
nvtxRangePop();

FLOP counting rule of thumb: a matrix multiply of shape (M, K) × (K, N) performs 2·M·K·N floating-point operations. Divide by wall-clock time to get achieved TFLOP/s, then compare against GPU peak (e.g., A100 = 77.6 TFLOP/s FP32) to compute arithmetic intensity.