平行化分析

[5000user5000]

層次 1：SIMD 向量化（資料層級平行化）

位置：SimdUtils.h - L2 距離計算

程式碼位置：include/zenann/SimdUtils.h:5-14

當前實現（Naive）：

  float d = 0.f;
  for (size_t i = 0; i < dim; ++i) {
      float diff = a[i] - b[i];
      d += diff * diff;
  }

平行化潛力：

• 高度適合 SIMD：完全獨立的資料運算
• 預期加速：4-8x（AVX2 每次處理 8 個 float）
• 技術：AVX2 intrinsics (mm256*)
• 影響範圍：所有距離計算（整個系統的瓶頸）

被調用處：

• IVFFlatIndex::search() - 第 52, 69 行
• IVFFlatIndex::train() - 第 32-34 行（內嵌版本）
• IVFFlatIndex::kmeans() - 第 144-147 行（內嵌版本）

優化方向：

  // 使用 AVX2 一次處理 8 個 float
  __m256 va = _mm256_loadu_ps(a + i);
  __m256 vb = _mm256_loadu_ps(b + i);
  __m256 diff = _mm256_sub_ps(va, vb);
  acc = _mm256_fmadd_ps(diff, diff, acc);  // FMA：融合乘加

層次 2：OpenMP 多線程（任務層級平行化）

2.1 IVFFlatIndex::search() - 質心距離計算

位置：src/IVFFlatIndex.cpp:51-54

當前實現：

  for (size_t c = 0; c < nlist_; ++c) {
      float d = l2_naive(query.data(), centroids_[c].data(), dimension_);
      cdist[c] = {d, c};
  }

平行化潛力：

• 完美平行：無數據依賴（每個質心獨立計算）
• 預期加速：接近線性（N 核心 ≈ N 倍）
• 技術：OpenMP parallel for（static schedule）
  -條件：nlist 足夠大（通常 100+）

優化方向：

  #pragma omp parallel for schedule(static)
  for (size_t c = 0; c < nlist_; ++c) {
      float d = l2_simd(...);  // 結合 SIMD
      cdist[c] = {d, c};
  }

2.2 IVFFlatIndex::search() - 列表探測

位置：src/IVFFlatIndex.cpp:65-81

當前實現：

  for (size_t pi = 0; pi < nprobe_; ++pi) {
      size_t c = cdist[pi].second;
      for (size_t id : lists_[c]) {
          float dist = l2_naive(...);
          // 更新 heap
      }
  }

平行化潛力：

• 中等複雜度：有共享數據（heap）
• 預期加速：3-8x（需要 thread-local heap + 合併）
• 技術：OpenMP parallel for + critical section
• 策略：每個線程維護 local heap，最後合併

優化方向：

  #pragma omp parallel for schedule(dynamic)
  for (size_t pi = 0; pi < nprobe_; ++pi) {
      std::vector<Pair> local_heap;  // thread-local
      // ... 搜尋 ...
      #pragma omp critical
      {
          // 合併到全局 heap
      }
  }

2.3 IVFFlatIndex::search_batch() - 批量查詢

位置：src/IVFFlatIndex.cpp:108-110

當前實現：

  for (size_t i = 0; i < nq; ++i) {
      results[i] = search(queries[i], k);
  }

平行化潛力：

• 最容易平行化：完全獨立的查詢
• 預期加速：線性（N 核心 ≈ N 倍）
• 技術：OpenMP parallel for（dynamic schedule）
• 高影響：批量查詢是典型使用場景

優化方向：

  #pragma omp parallel for schedule(dynamic)
  for (size_t i = 0; i < nq; ++i) {
      results[i] = search(queries[i], k);
  }

2.4 IVFFlatIndex::kmeans() - K-means 訓練

位置：src/IVFFlatIndex.cpp:139-155 和 160-165

2.4.1 Assignment 階段（E-step）

  // 第 139-155 行
  for (size_t i = 0; i < n; ++i) {
      // 為每個點找最近的質心
  }

平行化潛力：

• 完美平行：每個點獨立計算
• 預期加速：線性
• 技術：OpenMP parallel for

2.4.2 Update 階段（M-step）

  // 第 160-165 行
  for (size_t i = 0; i < n; ++i) {
      size_t c = assignment[i];
      for (size_t k = 0; k < dimension_; ++k) {
          sums[c][k] += data[i][k];  // ⚠️ Race condition
      }
  }

平行化潛力：

• 需要同步：多線程可能寫入同一個 sums[c]
• 預期加速：2-4x（需要 reduction 或 atomic）
• 技術：OpenMP reduction 或 thread-local 累加

優化方向：

// E-step
 #pragma omp parallel for
 for (size_t i = 0; i < n; ++i) { ... }

 // M-step - 使用 reduction
 #pragma omp parallel for reduction(+:sums[:nlist_*dimension_])
 for (size_t i = 0; i < n; ++i) { ... }

2.5 IVFFlatIndex::train() - 構建倒排表

位置：src/IVFFlatIndex.cpp:26-42

當前實現：

  for (size_t id = 0; id < data.size(); ++id) {
      // 為每個向量找最近的質心
      for (size_t c = 0; c < nlist_; ++c) {
          // 距離計算（可 SIMD）
      }
      lists_[best_c].push_back(id);  //  Race condition
  }

平行化潛力：

• 中等難度：寫入共享的 lists_
• 預期加速：3-6x
• 技術：OpenMP + atomic 或 thread-local 再合併

層次 3：其他平行化機會

3.1 Pipeline 平行化

概念：將索引構建和搜尋流水線化

場景：

• 邊構建邊搜尋（適合動態數據）
• I/O 與計算重疊（讀取數據 + 計算距離）

3.2 GPU 加速（未來方向）

適合 GPU 的操作：

• ✅ 大規模距離矩陣計算
• ✅ K-means 聚類（大數據集）
• ❌ 圖遍歷（HNSW）- 不規則訪問模式

平行化優先級

高優先級（容易 + 高回報）

位置	技術	難度	預期加速	理由
1. L2 距離（SIMD）	AVX2	⭐⭐	4-8x	系統瓶頸，影響所有操作
2. 批量查詢	OpenMP	⭐	Nx	最簡單，線性加速
3. 質心距離計算	OpenMP	⭐	Nx	簡單，無同步問題

中優先級（中等難度 + 中等回報）

位置	技術	難度	預期加速	挑戰
4. 列表探測	OpenMP + critical	⭐⭐⭐	3-8x	需要 heap 合併邏輯
5. K-means E-step	OpenMP	⭐⭐	Nx	簡單平行化

低優先級（複雜 + 回報不確定）

位置	技術	難度	預期加速	挑戰
6. K-means M-step	OpenMP reduction	⭐⭐⭐⭐	2-4x	需要處理 race condition
7. 構建倒排表	OpenMP + atomic	⭐⭐⭐⭐	3-6x	不規則寫入模式

組合優化效果預估

假設 8 核心 CPU + AVX2：

優化組合	預期總加速	備註
僅 SIMD	4-8x	距離計算占主要時間
僅 OpenMP（批量）	8x	僅限批量查詢
SIMD + OpenMP（質心）	10-15x	相乘效果
SIMD + OpenMP（全部）	20-50x	最佳情況

實作順序

階段 1：基礎 SIMD（1-2 天）
SimdUtils.h: l2_naive → l2_simd (AVX2)

階段 2：簡單 OpenMP（1 天）

1. search_batch() 批量查詢
2. search() 質心距離計算

階段 3：複雜 OpenMP（2-3 天）
3. search() 列表探測（thread-local heap）
4. kmeans() E-step

階段 4：進階優化（選做）
5. kmeans() M-step (reduction)
6. train() 倒排表構建

[5000user5000]

CUDA 加速面向筆記

1️⃣ 批量距離計算

為什麼適合：

• 大量獨立的距離計算（天生並行）
• 計算密集、記憶體訪問規則
• 批量操作可分攤 GPU 記憶體傳輸成本

應該加速的地方：
IVFFlatIndex.cpp:54-60
計算 query 到所有質心的距離（nlist 個距離計算）

CUDA 優勢：

• 假設 nlist=1000，dimension=128，可以啟動 1000+ 個 thread 並行計算
• 比 OpenMP 的 8-16 核心並行度高得多

2️⃣ 批量查詢優化（GPU 最佳使用場景）

為什麼重要：
IVFFlatIndex.cpp:155-168
批量查詢：每個 query 獨立處理

CUDA 優勢：

• 100 個查詢 × 1000 個質心 = 100,000 個距離計算
• GPU 可以同時處理所有計算
• 數據傳輸成本被分攤

3️⃣ K-means 訓練加速（中等優先）

IVFFlatIndex.cpp:194-210
Assignment step: 每個向量找最近的質心

CUDA 優勢：

• n 個向量完全獨立，可完全並行
• 訓練通常只做一次，但對大數據集仍有價值

實作建議

📋 CUDA 實作優先順序建議

● 🏗️ 實作架構設計

1. 檔案結構（類似現有 SIMD 設計）

  include/zenann/
    ├── SimdUtils.h      # 現有的 SIMD 優化
    └── CudaUtils.h      # 新增：CUDA kernel 聲明

  src/
    ├── CudaUtils.cu     # 新增：CUDA kernel 實作
    └── IVFFlatIndex.cpp # 修改：加入 CUDA 分支

2. 條件編譯策略

// CudaUtils.h
#pragma once
#include <cstddef>

namespace zenann {

#ifdef ENABLE_CUDA
// GPU 版本
void l2_distance_batch_gpu(
    const float* queries,     // [nq x dim]
    const float* centroids,   // [nc x dim]
    float* distances,         // [nq x nc] 輸出
    size_t nq,               // query 數量
    size_t nc,               // centroid 數量
    size_t dim               // 維度
);
#endif

// CPU fallback（使用現有的 l2_distance）
inline void l2_distance_batch_cpu(/* ... */) {
    // 使用現有的 OpenMP/SIMD 實作
}

} // namespace zenann

3. IVFFlatIndex 整合

// src/IVFFlatIndex.cpp
SearchResult IVFFlatIndex::search(const Vector& query, size_t k) const {
    // ...

#ifdef ENABLE_CUDA
    // 使用 GPU 計算所有質心距離
    if (nlist_ > CUDA_THRESHOLD) {  // 例如 nlist_ > 100
        l2_distance_batch_gpu(query.data(), centroids_data,
                             distances, 1, nlist_, dimension_);
    } else {
        // Fallback to CPU
#endif

#ifdef ENABLE_OPENMP
    #pragma omp parallel for schedule(static)
#endif
    for (size_t c = 0; c < nlist_; ++c) {
        float d = l2_distance(query.data(), centroids_[c].data(), dimension_);
        cdist[c] = {d, c};
    }

#ifdef ENABLE_CUDA
    }
#endif
}

📊 性能預期

場景	CPU (OpenMP+SIMD)	GPU (CUDA)	預期加速比
單查詢 (nlist=100)	基準	可能更慢*	0.5x
單查詢 (nlist=1000)	基準	較快	2-3x
批量查詢 (100個, nlist=1000)	基準	顯著加速	10-50x
K-means 訓練 (10000向量)	基準	顯著加速	20-100x

註： GPU 有記憶體傳輸開銷，小規模計算反而可能更慢

🎯 關鍵設計決策

1. 何時使用 GPU？

  // 自動選擇策略
  bool should_use_gpu(size_t nq, size_t nc, size_t dim) {
      // 計算量夠大才使用 GPU
      size_t compute_ops = nq * nc * dim;
      return compute_ops > 1000000;  // 100萬次 FLOPs
  }

2. 記憶體管理

  class CudaMemoryManager {
      // 預分配 GPU memory pool
      // 減少頻繁的 cudaMalloc/cudaFree
  };

3. 混合策略

• 小數據集: CPU (OpenMP + SIMD)
• 大數據集: GPU (CUDA)
• 批量查詢: 優先 GPU

🔧 Makefile 修改

需要加入 CUDA 編譯器

  NVCC := nvcc
  NVCC_FLAGS := -std=c++17 -O3 --compiler-options -fPIC

  # CUDA sources
  CUDA_SOURCES := src/CudaUtils.cu
  CUDA_OBJECTS := $(CUDA_SOURCES:.cu=.o)

  # CUDA version
  cuda: prepare $(CUDA_OBJECTS)
        $(CXX) $(FULL_CXXFLAGS) -DENABLE_CUDA $(ALL_INCLUDES) -shared -o $(TARGET) \
            $(SOURCES) $(CUDA_OBJECTS) \
            -L$(FAISS_ROOT)/lib -lfaiss \
            -L/usr/local/cuda/lib64 -lcudart \
            $(ALL_LIBS) \
            $(FULL_LDFLAGS)

  %.o: %.cu
        $(NVCC) $(NVCC_FLAGS) -DENABLE_CUDA -c $< -o $@ $(PROJECT_INCLUDE)

總結建議

✅ 優先實作

1. GPU 批量距離計算 - 最大收益
2. 批量查詢優化 - 實用性最高
3. 自動 CPU/GPU 選擇 - 提升使用體驗

⚠️ 注意事項

• GPU 記憶體傳輸是瓶頸，需要批量操作
• 小規模計算不要用 GPU（反而變慢）
• 考慮使用 CUDA Streams 重疊計算和傳輸

📈 預期效果

• 批量查詢場景: 10-50x 加速
• 訓練場景: 20-100x 加速
• 單查詢小數據: 可能更慢（fallback to CPU）

[5000user5000]

注意 openMP 和 gpu 可能會互斥,
這裡的做法可以是編譯 cuda 時不使用 openmp simd,
如果要混合(make full),可以 GPU 做批次距離計算, openMP 做結果處理和 list probing, 或是如上文所敘,根據數據大小來判斷.

[5000user5000]

或是開一個分支,把 openMP SIMD 拔掉,只用 cuda

[Seco1024]

Profiling OpenMP + SIMD：

nlist = 1024, nprobe = 10

Stage	Mean	Min	Max	% Total
Centroid Distance	0.6273	0.0215	11.7319	42.8%
Centroid Selection	0.0027	0.0020	0.0041	0.2%
List Scanning	0.8368	0.1725	8.6686	57.0%
Final Sorting	0.0003	0.0001	0.0007	0.0%
TOTAL	1.4671	0.1987	16.3619	100%

nlist = 1024, nprobe = 50

Stage	Mean	Min	Max	% Total
Centroid Distance	0.3880	0.0172	8.0563	31.7%
Centroid Selection	0.0088	0.0068	0.0119	0.7%
List Scanning	0.8255	0.2709	12.6484	67.5%
Final Sorting	0.0004	0.0002	0.0012	0.0%
TOTAL	1.2227	0.3026	19.9063	100%

nlist = 1024, nprobe = 100

Stage	Mean	Min	Max	% Total
Centroid Distance	0.3355	0.0166	8.4438	23.5%
Centroid Selection	0.0152	0.0117	0.0244	1.1%
List Scanning	1.0791	0.4304	19.6904	75.4%
Final Sorting	0.0004	0.0002	0.0013	0.0%
TOTAL	1.4303	0.4650	22.5965	100%

CUDA 加速可以以 Centroid Distance （coarse-grained search）+ List Scanning （fine-grained search）為主。

[Seco1024]

11/25 update

已完成 CUDA 版本的大致雛形，
四個 kernel function 對照到上面 profiling 數據的四個階段。

現在的問題是 k=100 會產生 shared memory 不夠的 error 而無法執行，
以及效能會變很差，可以再看看怎麼 workaround。