Linux 核心專題: llama.cpp 效能分析

# Linux 核心專題: llama.cpp 效能分析 > 執行人: vestata > [專題解說影片](https://youtu.be/WM6qmB4tQJs) ### Reviewed by `yu-hsiennn` > 整個程式運作都的計算熱點會是矩陣乘法指的是整個程式運作的計算熱點都會是矩陣乘法? 還有 `熱點` 指的是什麼? > [Wikipedia: hot spot](https://en.wikipedia.org/wiki/Hot_spot_(computer_programming)) 我翻譯為「熱點」，解釋為程式執行過程中大部分時間消耗的區域，〈[LLaMA Now Goes Faster on CPUs](https://justine.lol/matmul/)〉中使用的描述是 [Wikipedia: Bottleneck](https://en.wikipedia.org/wiki/Bottleneck_(software)) 但是我覺得在這邊使用熱點比較符合。 ### Reviewed by `yenslife` 既然知道 `((1 << 3) - 1)` 表示 7，那為什麼一開始不要直接寫 7 就好？我剛剛在猜會不會與操作成本有關，但看到之後開發者修改成 -8 就更不懂了。 > 使用 `((1 << 3) - 1)` 表示為 7 而不是直接使用 7 的具體原因...我還不知道，從第一個版本 [commit 26c0846](https://github.com/ggerganov/llama.cpp/commit/26c084662903ddaca19bef982831bfb0856e8257#diff-6d9ce99fcb6f51ff76f59e479f6e6fc0bb62edef7442805d7a5bb15b23996b5d) 便是使用 `((1 << 3) - 1)` 延用到現在，我認為是為了增加程式可解讀性。若有發現我會再更新。 > 可以參考我之後的在 [Quantization](https://hackmd.io/8VsaB9JESziJyp5lfmSvRQ#Quantization) 下面的補充，其中 $\Delta$ 是其步長，一般會定義為 $\Delta = 2^{(b-1)}-1$ ，在 `q4_0` 也就是 `((1 << 3) - 1)` ，但是因為使用二補數，所以 4 位元表示範圍應該是 -8~7。我們現在令 `n = 30000, max = 65504(FP16 表示的最大值)` 作為舉例，分別進行 -8 與７的量化與解量化為 > 7:((30000 / (65504 / 7) + 8.5) - 8) * (65504 / 7) = 34678.85714285714 > -8:((30000 / (65504 / -8) + 8.5) - 8) * (65504 / -8) = 25905.999999999996 誤差分別為 -4678.857 與 4094，使用 -8 有較小的誤差。 ### Reviewed by `jserv` 將 [matmul.c](https://github.com/salykova/matmul.c) 納入考量，這是針對現代 x86 處理器的 AVX 和 FMA 指令集擴展開發的快速矩陣乘法實作，程式碼精簡又性能表現優異。參照該專案的介紹文章 [Beating NumPy's matrix multiplication in 150 lines of C code](https://salykova.github.io/matmul-cpu)。 ## 任務簡介探討 llama.cpp 的效能表現，分析其效能瓶頸並嘗試改進 ## TODO: 讀懂〈[LLaMA Now Goes Faster on CPUs](https://justine.lol/matmul/)〉並重現實驗 > 不同 quantization 所產生的精度差異，如何實驗它對模型的影響？ :::danger 改進你的漢語表達。 ::: [llamafile](https://github.com/Mozilla-Ocho/llamafile) 是為了簡化使用 LLM 所生的專案。結合了 [llama.cpp](https://github.com/ggerganov/llama.cpp) 與 [Cosmopolitan Libc](https://github.com/jart/cosmopolitan) 使配置與運行壓縮到單一可執行檔中。 :::danger 注意用語： * file 是「檔案」，而非「文件」(document) * thread 是「執行緒」，而非「線程」 > 好的，收到 ::: ### Quantization 〈[LLaMA Now Goes Faster on CPUs](https://justine.lol/matmul/)〉這篇文章是中使用不同的量化手法進行效能的比較，主要針對 prompt (tok/sec) 即讀取 token 的速度，與 eval (tok/sec) 輸出的 token 速度進行比較。一般訓練好並開放的大型語言模型權重 `.pth` 是使用 BF16 的資料型態，在 [llama.cpp](https://github.com/ggerganov/llama.cpp) 需要將 `.pth` 權重轉成 `.gguf` 的資料型態。 `.gguf` 是 [llama.cpp](https://github.com/ggerganov/llama.cpp) 的作者所開發的資料格式，可以使用在其相關的專案中，如以下圖例： ![image](https://hackmd.io/_uploads/r1bLTzDSR.png) > https://huggingface.co/docs/hub/gguf :::info TODO gguf/tensor 的差異補 gguf 的 `mmap()` ::: 依據〈[Quantizing deep convolutional networks for efficient inference: A whitepaper](https://arxiv.org/pdf/1806.08342)〉這篇論文，可以將量化的分為兩種 Affine Quantizer 與 Symmetric Quantizer 兩種。 ##### Symmetric quantize \begin{equation} x_{\text{int}} = \text{round}\left(\frac{x}{\Delta}\right) \end{equation} \begin{equation} x_Q = \text{clamp}\left(-\frac{N_{\text{levels}}}{2}, \frac{N_{\text{levels}}}{2} - 1, x_{\text{int}}\right) \quad \text{if signed} \end{equation} \begin{equation} x_Q = \text{clamp}(0, N_{\text{levels}} - 1, x_{\text{int}}) \quad \text{if un-signed} \end{equation} 其中 \begin{equation*} \text{clamp}(a, b, x) = \begin{cases} a & \text{if } x \leq a \\ x & \text{if } a \leq x \leq b \\ b & \text{if } x \geq b \end{cases} \end{equation*} s 為scale factor，z 為 zero point 此方法會對應到等等的 `Q8_0` 手法。 ##### Affine Quantizer 量化： \begin{equation} x_{\text{int}} = \text{round}\left(\frac{x}{\Delta}\right) + z \end{equation} \begin{equation} x_Q = \text{clamp}(0, N_{\text{levels}} - 1, x_{\text{int}}) \end{equation} 解量化： $x_{float}=(x_Q-z)\Delta$ 此方法會對應到等等的 `Q4_0`、`Q4_1` 手法。 #### 量化手法介紹手法通常分為 `Qn_0`，`Qn_1`，`Qn_k`。 - [ ] `Qn_0` 將權重縮小到 n 位元表示，一次進行一個 block，大小為 32，選擇每一 block 局部最大值進行縮放，區間為 0 到最大值之間。 - [ ] `Qn_1` 將權重縮小到 n 位元表示，一次進行一個 block，大小為 32，選擇每一 block 局部最大值與最小值進行縮放，區間為最小值到最大值之間。 - [ ] `Qn_k` :::info TODO ::: `Q4_0` 會儲存在一個 `block_q4_0` 結構體中，其大小會是 32，其中 `d` 是他的縮放比例（scalar），`qs` 是一個 `uint8_t` 陣列用來儲存兩個 4 位元的，在量化的過程中，原本 BF16 的模型權重以 `block_q4_0` 為單位去做儲存： ```c #define QK4_0 32 typedef struct { ggml_half d; // delta is uint16 uint8_t qs[QK4_0 / 2]; // nibbles / quants } block_q4_0; ``` 量化函式 `quantize_row_q4_0_reference()` 是負責量化的程式，其中發現很有趣的地方。以下是量化函式中的細節 `const float d = max / -8;` 是一個在 [commit dd0eabc](https://github.com/ggerganov/llama.cpp/commit/dd0eabc049fb1efc631cab8eb0a646808d704e18) 的變更，將 `((1 << 3) - 1)` 改為 -8 因為這邊是為了將原本 FP16 量化到 4 個位元，所以 4 位元的二的補數表示範圍為 -8 到 7。為了善用更大的分母所以使用 -8，其所導致的整個存放在 `.gguf` 的權重變成負數，在 perplexity 的表現上勝於 `((1 << 3) - 1)` 的方法，見 [#729](https://github.com/ggerganov/llama.cpp/pull/729)。而在後面要將其放進 4 位元之中，將原本權重參數乘上縮放比例 `id`，但是後面的順序的擺放是將 `x[i*qk + 0 + j]` 與 `x[i*qk + qk/2 + j]` 擺在一起。這是在 [commit b9fd7ee](https://github.com/ggerganov/llama.cpp/commit/b9fd7eee57df101d4a3e3eabc9fd6c2cb13c9ca1) 的變更。 ```diff y[i].d = d; - for (int l = 0; l < QK4_0; l += 2) { - const float v0 = x[i*QK4_0 + l + 0]*id; - const float v1 = x[i*QK4_0 + l + 1]*id; + for (int j = 0; j < qk/2; ++j) { + const float x0 = x[i*qk + 0 + j]*id; + const float x1 = x[i*qk + qk/2 + j]*id; - const uint8_t vi0 = MIN(15, (int8_t)roundf(v0) + 8); - const uint8_t vi1 = MIN(15, (int8_t)roundf(v1) + 8); + const uint8_t xi0 = MIN(15, (int8_t)(x0 + 8.5f)); + const uint8_t xi1 = MIN(15, (int8_t)(x1 + 8.5f)); - assert(vi0 < 16); - assert(vi1 < 16); - - pp[l/2] = vi0 | (vi1 << 4); + y[i].qs[j] = xi0; + y[i].qs[j] |= xi1 << 4; } - - memcpy(y[i].qs, pp, sizeof(pp)); } ``` 這個變更讓我困惑很久，但是這必須提到 `llama.cpp` 的開發規則，在此專案之中的定義的矩陣相乘是 $C^T=AB^T \Leftrightarrow C = BA^T$ ，這與傳統的矩陣相乘記憶體儲存手法不同。 ![image](https://hackmd.io/_uploads/ry1NAY_UC.png) > [README.md](https://github.com/ggerganov/llama.cpp/blob/master/README.md#coding-guidelines) 在這邊補充一下，一般我們在做矩陣相乘的時候 $C = AB$ 會遇到的狀態會是 $A,B$ 在記憶體的排列是為 row-major。當處理 $A$ 時，可以利用其連續記憶體的特性，減少到主記憶體的成本，這樣能夠利用 L1 cache。然而，$B$ 的矩陣乘法需要每一 column 的資料，這與 row-major 記憶體排列剛好相反，所以 $C = AB$ 這個乘法最大的效能瓶頸會是 $B$ 對主記憶體成本。有很多矩陣乘法手段可以更好的利用 L1 cache，一個簡單的方法是對 $B$ 做轉置，詳細可以參考我對[矩陣相乘的筆記](https://hackmd.io/@tony268596/Hy4EWvXM0)。所以對於 `quantize_row_q4_0_reference()` 的變更就是為了符合在 $C = BA^T$ 的運算規則之下。對其中一個矩陣雖然因為 L1 cache 的優勢，但是以演算法的角度尚需要 $O(n^3)$ 的成本。因此，llama.cpp 原本的作法是將量化後的 `uint` 以相鄰的方式排列，但之後需要額外的轉置成本。改動後，量化時便將其以矩陣轉置 $A^T$ 的排序放進記憶體中，以減少後續的轉置成本。 ### Perplexity 經過量化數值的精確度一定會喪失，只是不同的量化手段會有不同精確度喪失，但是這精確度喪失會大型語言模型會有什麼影響呢？回到模型的根本，是使用模型權重去選擇生成下一個最有可能出現的 token，我們可以針對生出的 token 進行評估。在 llama.cpp 專案之中的方法是使用 perplexity。其定義如下： $$\text{PPL}(X) = \exp \left\{ -\frac{1}{t} \sum_{i}^t \log p_\theta (x_i \mid x_{<i}) \right\} $$ :::danger 闡述其特性、考量因素，以及對應的限制。 ::: 給定特定的文章，這邊使用的是 [WikiText-2](https://paperswithcode.com/dataset/wikitext-2) 在給定的文章中，模型使用 $X_{<i}$ 的已知 token 去推斷 $X_i$ 的 token，$t$ 表示為 token 數，所以 perplexity 越低，表示模型的預測結果與實際 token 序列越接近。 ![image](https://hackmd.io/_uploads/H1xvZ8nL0.png =70%x) 文章中使用不同的 quantization 手法，有 `Q4_0`，`Q8_0`，`FP16`。以 `Q4_0` 為例，以下是主要實作的程式： ### 重現實驗使用的測試工具是 llama-bench，這是在專案之中被廣泛的使用的測試工具，分別會計算輸入 token(prompt) 與輸出 token(eval) 的時間。主要的評估函式是 `test_promt()` 與 `test_gen()` 。其中 `test_promt()` 是以 `rand()` 去生成隨機數字，計算模型在消化 n_prompt 的隨機數字需要多少時間。 ```c while (n_processed < n_prompt) { int n_tokens = std::min(n_prompt - n_processed, n_batch); tokens[0] = n_processed == 0 && llama_add_bos_token(model) ? llama_token_bos(model) : std::rand() % n_vocab; for (int i = 1; i < n_tokens; i++) { tokens[i] = std::rand() % n_vocab; } llama_decode(ctx, llama_batch_get_one(tokens.data(), n_tokens, n_past + n_processed, 0)); n_processed += n_tokens; } ``` 在 `test_gen()` 則是計算模型每次處理一個 `rand()` 隨機數，進行 n_gen 次所花費的時間。 ```c for (int i = 0; i < n_gen; i++) { llama_decode(ctx, llama_batch_get_one(&token, 1, n_past + i, 0)); token = std::rand() % n_vocab; } ``` 測試設備： ```bash Architecture: x86_64 CPU op-mode(s): 32-bit, 64-bit Address sizes: 46 bits physical, 48 bits virtual Byte Order: Little Endian CPU(s): 24 On-line CPU(s) list: 0-23 Vendor ID: GenuineIntel Model name: 13th Gen Intel(R) Core(TM) i7-13700 CPU family: 6 Model: 183 Thread(s) per core: 2 Core(s) per socket: 16 Socket(s): 1 Stepping: 1 CPU max MHz: 5200.0000 CPU min MHz: 800.0000 BogoMIPS: 4224.00 Virtualization features: Virtualization: VT-x Caches (sum of all): L1d: 640 KiB (16 instances) L1i: 768 KiB (16 instances) L2: 24 MiB (10 instances) L3: 30 MiB (1 instance) ``` 我這邊選擇兩個版本一個是 [commit dbceec8](https://github.com/ggerganov/llama.cpp/commit/dbceec87c0221ec952e69448df6a71f1372a7487) 使用舊版向量內積，與 [commit 74f33ad](https://github.com/ggerganov/llama.cpp/commit/74f33adf5f8b20b08fc5a6aa17ce081abe86ef2f)，此版本是 2024-5-23 是我拿來做測試的版本，使用的是 sgemm 新的矩陣乘法模組（llamafile 的乘法模組）。我將以這兩個做實驗比較。在 x86 機器中使用 `make` 即可完成編譯，其中執行緒的數量皆是使用預設。 :::danger 補上對應的超連結。 > 已補上 ::: | prompt (tok/sec) | eval (tok/sec) | model | weights data type | hardware | software | | ---------------- | -------------- | -------------- | ----------------- | ----------- | ------------------- | | 255.90 ± 12.60 | 26.29 ± 0.56 | TinyLlama 1.1B | q8_0 | Raptor Lake | llama.cpp 2024-5-23 | | 129.67 ± 2.42 | 22.70 ± 0.61 | TinyLlama 1.1B | q8_0 | Raptor Lake | llama.cpp old | | 199.98 ± 33.37 | 14.96 ± 0.36 | TinyLlama 1.1B | FP16 | Raptor Lake | llama.cpp 2024-5-23 | | 48.99 ± 0.19 | 12.97 ± 0.27 | TinyLlama 1.1B | FP16 | Raptor Lake | llama.cpp old | 在 `TinyLlama 1.1B` 權重之下 `q8_0`，`FP16` 使用 llamafile 加速後，在 `prompt (tok/sec)` 上的性能分別提升了 197% 和 414%，有顯著的進步。，在 `eval (tok/sec)` 的表現上則相近。 | prompt (tok/sec) | eval (tok/sec) | model | weights data type | hardware | software | | ---------------- | -------------- | ---------- | ----------------- | ----------- | ------------------- | | 25.08 ± 1.10 | 2.05 ± 0.08 | mistral-7b | q8_0 | Raptor Lake | llama.cpp 2024-5-23 | | 19.40 ± 0.75 | 3.55 ± 0.02 | mistral-7b | q8_0 | Raptor Lake | llama.cpp old | | 41.51 ± 0.68 | 4.21 ± 0.05 | mistral-7b | FP16 | Raptor Lake | llama.cpp 2024-5-23 | | 6.86 ± 0.07 | 1.84 ± 0.05 | mistral-7b | FP16 | Raptor Lake | llama.cpp old | 在 `mistral-7b` 權重之下 `q4_0`，`q8_0`，`FP16` 使用 llamafile 加速後，在 `prompt (tok/sec)` 上的性能分別提升了 440%、%、413% 和 364%，有顯著的進步。，在 `eval (tok/sec)` 的表現上則相近。 :::danger 討論其精確度的犧牲狀況。 :::   ## TODO: 影響 LLaMA 推理速度的因素 > 找出其中效能瓶頸的地方，可能是矩陣運算，或是記憶體 (特別是 data model)，予以量化，善用 perf 和 Ftrace 等工具 ### 矩陣相乘使用 perf 依照 `cpu-cycles`，`faults`，`cache-references`，`cache-misses` 進行比較 ```bash $ sudo perf record -e cpu-cycles,faults,cache-references,cache-misses ./llama-bench -m models/tinyllama-1.1b-chat-v1.0.Q4_0.gguf $ sudo perf report ``` 在 llama.cpp 併入 llamafile 的 sgemm 之前，也就是 [commit dbceec8](https://github.com/ggerganov/llama.cpp/commit/dbceec87c0221ec952e69448df6a71f1372a7487) 的 backend cpu 實作手法，整個程式運作都的計算熱點會是矩陣乘法： ``` # Samples: 1M of event 'cpu_core/cpu-cycles/' 53.96% ,llama-bench,llama-bench ,[.] ggml_graph_compute_thread 37.68% ,llama-bench,llama-bench ,[.] ggml_vec_dot_q4_0_q8_0 1.86% ,llama-bench,llama-bench ,[.] ggml_vec_dot_q6_K_q8_K 1.60% ,llama-bench,llama-bench ,[.] ggml_vec_dot_f16 0.66% ,llama-bench,llama-bench ,[.] ggml_compute_forward_mul_mat # Samples: 699K of event 'cpu_atom/cpu-cycles/' 78.01% ,llama-bench,llama-bench ,[.] ggml_vec_dot_q4_0_q8_0 6.88% ,llama-bench,llama-bench ,[.] ggml_vec_dot_f16 6.21% ,llama-bench,llama-bench ,[.] ggml_graph_compute_thread 2.67% ,llama-bench,llama-bench ,[.] ggml_vec_dot_q6_K_q8_K 0.95% ,llama-bench,llama-bench ,[.] ggml_compute_forward_mul_mat 0.74% ,llama-bench,llama-bench ,[.] ggml_compute_forward_soft_max ``` ### 矩陣計算底層程式在 llama.cpp 中有 gg 作者所撰寫的手法是簡單的將矩陣轉為內積去做計算，由上述提到的 $C=BA^T$ 並且在量化的時候便已經將 $A$ 向量化排列在連續記憶體之中。有很多開源的 BLAS libraries 有很純熟的矩陣相乘手法，但是放在已經向量化的場景反而會降低速度，因為這些 BLAS libraries 有很多複雜的手段會增加許多延遲，是對應各層面的矩陣相乘場景，而 gg 作者的向量內積是為了這個專案特別設計的。其中 `ggml_vec_dot_q4_0_q8_0()` 是實作整個向量內積的底層程式。但是細看程式碼，先以矩陣相乘的演算法的角度（SIMD 等技巧先忽略），我們可以發現這是一個 q4_0 乘上 q8_0 的混合精度乘法： ```c void ggml_vec_dot_q4_0_q8_0(int n, float * restrict s, size_t bs, const void * restrict vx, size_t bx, const void * restrict vy, size_t by, int nrc) { const int qk = QK8_0; const int nb = n / qk; const block_q4_0 * restrict x = vx; const block_q8_0 * restrict y = vy; // scalar float sumf = 0.0; for (int i = 0; i < nb; i++) { int sumi = 0; for (int j = 0; j < qk/2; ++j) { const int v0 = (x[i].qs[j] & 0x0F) - 8; const int v1 = (x[i].qs[j] >> 4) - 8; sumi += (v0 * y[i].qs[j]) + (v1 * y[i].qs[j + qk/2]); } sumf += sumi*GGML_FP16_TO_FP32(x[i].d)*GGML_FP16_TO_FP32(y[i].d); // multipyy the scalars } *s = sumf; } ``` 此類函式是整個專案中負責最重要的矩陣計算，所以 gg 作者希望可以將這段程式最佳化。它負責計算 `z = x * y`，其中 `x` 為量化完的數值，`y` 是 FP32，`z` 為 FP32。這邊就令人好奇為什麼是 q4_0 * q8_0 的混合精度乘法呢？為什麼不是直覺的 q4_0 * q4_0 呢？在一開始的確是使用 `ggml_vec_dot_q4_0()` 這個函式。 ```c typedef struct { ggml_half d; // delta uint8_t qs[QK4_0 / 2]; // nibbles / quants } block_q4_0; ``` ```c const block_q4_0 * restrict x = vx; const block_q4_0 * restrict y = vy; float sumf = 0.0; ... for (int i = 0; i < nb; i++) { const float d0 = x[i].d; const float d1 = y[i].d; const uint8_t * restrict p0 = x[i].qs; const uint8_t * restrict p1 = y[i].qs; int sumi = 0; for (int j = 0; j < QK/2; j++) { const uint8_t v0 = p0[j]; const uint8_t v1 = p1[j]; const int8_t i0 = (int8_t) (v0 & 0xf) - 8; const int8_t i1 = (int8_t) (v0 >> 4) - 8; const int8_t i2 = (int8_t) (v1 & 0xf) - 8; const int8_t i3 = (int8_t) (v1 >> 4) - 8; sumi += i0*i2 + i1*i3; } sumf += d0 * d1 * sumi; ``` 在 caller `ggml_compute_forward_mul_mat` 發現原先在計算 z = x * y 時，會先將 y (FP32) 做一次 q4_0 的量化。 ```c if (params->type == GGML_TASK_INIT) { char * wdata = params->wdata; const size_t row_size = ne10*GGML_TYPE_SIZE[type]/GGML_BLCK_SIZE[type]; for (int64_t i13 = 0; i13 < ne13; ++i13) { for (int64_t i12 = 0; i12 < ne12; ++i12) { for (int64_t i11 = 0; i11 < ne11; ++i11) { quantize_row_q((float *)((char *) src1->data + i13*nb13 + i12*nb12 + i11*nb11), (void *) wdata, ne10); wdata += row_size; } } } ``` 在這邊可以看到原本在計算 `y` 的步驟是：  ``` z(FP32) = x(Q4_0) * y(FP32) ^ | | | | |------這個 Q4_0 量化過程會導致精度的損失 | | | | x(int8_t) y(Q4_0) | | | | | |------內積時還是要將 4 bit 擴充成 int8_t 做計算 | | | --------------- * ---------- y(int8_t) ``` 改為 ``` z(FP32) = x(Q4_0) * y(FP32) ^ | | | | |------ 改為 Q8_0 的量化 | | | | x(int8_t) y(Q8_0) | | | | | | | | | --------------- * ---------- y(int8_t) ``` 這個步驟在每一個向量相乘都會發生，且因為他是一個底層運算函式，不斷的做量化會不斷放大誤差，導致模型的 perplexity 居高不下。 >But at the end, avoiding quantization of intermediate results will always be more important than any improvements to model weight quantization. 這是在此專案場景矩陣相乘遇到的困難點，所以才衍生出之後的 `ggml_vec_dot_q4_0_q8_0()`，降低了量化過程中所帶來的精度喪失。在之後底層運算的實作皆是使用這個向量乘法演算法。 :::danger 建立數學模型，探討量化過程的精確度損失和補償方案。數學分析！ ::: ```c void ggml_vec_dot_q4_0_q8_0(int n, float * restrict s, size_t bs, const void * restrict vx, size_t bx, const void * restrict vy, size_t by, int nrc) { const int qk = QK8_0; const int nb = n / qk; assert(n % qk == 0); #if defined(__ARM_FEATURE_MATMUL_INT8) assert((nrc == 2) || (nrc == 1)); #else assert(nrc == 1); #endif UNUSED(nrc); UNUSED(bx); UNUSED(by); UNUSED(bs); const block_q4_0 * restrict x = vx; const block_q8_0 * restrict y = vy; // scalar float sumf = 0.0; for (int i = 0; i < nb; i++) { int sumi = 0; for (int j = 0; j < qk/2; ++j) { const int v0 = (x[i].qs[j] & 0x0F) - 8; const int v1 = (x[i].qs[j] >> 4) - 8; sumi += (v0 * y[i].qs[j]) + (v1 * y[i].qs[j + qk/2]); } sumf += sumi*GGML_FP16_TO_FP32(x[i].d)*GGML_FP16_TO_FP32(y[i].d); } *s = sumf; } ``` ### 矩陣計算 thread 執行緒的工作分為三種，`GGML_TASK_TYPE_INIT`，`GGML_TASK_TYPE_COMPUTE`，`GGML_TASK_TYPE_FINALIZE`，分別處理工作的初始化、矩陣的計算以及不同執行緒結果的收尾。在這邊使用的同步機制是 spinlock，執行緒會上鎖並分頭計算不同區塊的矩陣。這邊並行程式是為了處理底層計算的部份，因此沒有使用 futexes 或是 semaphores，因為 Linux 核心的排程會降低計算效能。 ```c static void ggml_graph_compute_thread_sync_node(int * node_n, struct ggml_compute_state * state, const bool do_yield) { // wait for other threads to finish const int last_node_n = * node_n; while (true) { if (do_yield) { sched_yield(); } * node_n = atomic_load(&state->shared->node_n); if (* node_n != last_node_n) break; } } ``` 首先紀錄 node_n 編號。這段程式會一直查看佇列中是否有新的工作，使用 atomic_load 查看若是 node_n 值有所不同，表示其他執行緒已更新了節點，可以開始處理其他節點了。如果設定 `do_yeild` 為真它會放棄目前 CPU 然後移到排程佇列最末端，讓其他執行緒優先執行。 :::warning 搭配 taskset 觀察多執行緒程式的執行。 ::: #### `ggml_graph_compute_thread()` 此函式負責的處理並行矩陣計算所要的執行緒所使用的所有動作，逐一解釋：翻閱 `atomic_fetch_sub()` man page 說明如下： >The atomic_fetch_sub() macro subtracts the value operand from atomic variable object and returns the original contents of the atomic variable. ```c if (atomic_fetch_sub(&state->shared->n_active, 1) == 1) ``` 所以如果 `atomic_fetch_sub()` 回傳是 1 的話，一方面執行 n_active 執行緒數減一，一方面表示這是矩陣乘法多執行緒中的最後一個執行緒，其他工作的執行緒已完成，最後一個執行緒會有一個收尾的動作。以下這一段是處理收尾的主要程式： ```c if (node_n != -1) { /* FINALIZE */ struct ggml_tensor * node = cgraph->nodes[node_n]; if (GGML_OP_HAS_FINALIZE[node->op]) { params.nth = ggml_get_n_tasks(node, n_threads, state->shared->n_threads); ggml_compute_forward(&params, node); } ggml_graph_compute_perf_stats_node(node, state->shared); } ``` `node_n` 初始化為 -1。調用 `ggml_get_n_tasks()` （存放所有動作會需要多少執行緒）來獲取現在的動作會需要多少執行緒。並使用 `ggml_compute_forward()` 完成最後的計算。 ## TODO: 探討 llamafile 相較於 llama.cpp 進行哪些調整，得以加速？ > 重現實驗並分析 sgemm.cpp 是 llamafile 的矩陣計算模組，大幅提昇 prompt evaluation 的速度，這個矩陣計算模組也併入 llam.cpp 之中，是目前 llam.cpp 的預設計算模組。在我使用原本 `vec_dot()` 版本進行效能測試時它預設開啟的執行緒數量是 12，但是當我到 [commit 8cc91dc](https://github.com/ggerganov/llama.cpp/commit/8cc91dc63c0df397d644a581b2cbeea74eb51ae0) 測試發現這邊預設的執行緒是 8。後來發現這邊也有緣由。在原本的版本使用 `get_num_physical_cores()` 去判斷現在的設備有多少邏輯核心，如果邏輯核心數超過 4 個便會使用 `邏輯核心數 / 2` 個執行緒。以下是我的設備為例有 24 核，所以在原本向量矩陣相乘的方法將 12 配置給執行中的程式。而 sgemm.cpp 改成 8 是因為 `13th Gen Intel(R) Core(TM) i7-13700` 其中 24 個核心包含 8 個雙執行緒的效能核心與 8 個單執行緒的效率核心，希望避免將大量的矩陣計算移到效率核心去做計算。 [執行緒選取的更動](https://github.com/ggerganov/llama.cpp/commit/8cc91dc63c0df397d644a581b2cbeea74eb51ae0#diff-201cbc8fd17750764ed4a0862232e81503550c201995e16dc2e2766754eaa57a) 這邊有對核心的選取，我對這邊核心相關的程式細節不了解(`cpuid()`, `is_hybrid_cpu()`, `is_running_on_efficiency_core()`)，希望知道人可以補充。這邊使用 `pthread_setaffinity_np()` 將設定在效能核心上。 ```c int get_math_cpu_count() { #if defined(__x86_64__) && defined(__linux__) int cpu_count = sysconf(_SC_NPROCESSORS_ONLN); if (cpu_count < 1) { return get_num_physical_cores(); } if (is_hybrid_cpu()) { cpu_set_t affinity; if (!pthread_getaffinity_np(pthread_self(), sizeof(affinity), &affinity)) { int result = count_math_cpus(cpu_count); pthread_setaffinity_np(pthread_self(), sizeof(affinity), &affinity); if (result > 0) { return result; } } } #endif return get_num_physical_cores(); } ``` > On success, these functions return 0; on error, they return a nonzero error number. > [name=man page] ### 矩陣相乘 llama.cpp 支援許多種 backend，也就是多種矩陣相乘的計算模組，原有支援 `CPU`，`OpenCL`，`Vulkan`，`Kompute`，`Metal`，`GPU BLAS`，`BLAS`。在 [commit 8cc91dc](https://github.com/ggerganov/llama.cpp/commit/8cc91dc63c0df397d644a581b2cbeea74eb51ae0) 加入 sgemm.cpp 這個乘法模組。文章中有提到其在 `Q8` 與 `FP16` 有最好的效果，那可以從其計算熱點觀察。 ```c 66.22% llama-bench llama-bench [.] void (anonymous namespace)::tinyBLAS_Q0_AVX<block_q8_0, block_q8_0, float>::gemm<4, 1>(long, long, long, long) ◆ 20.80% llama-bench llama-bench [.] void (anonymous namespace)::tinyBLAS_Q0_AVX<block_q8_0, block_q8_0, float>::gemm<4, 2>(long, long, long, long) 6.46% llama-bench llama-bench [.] ggml_graph_compute_thread ▒ 1.43% llama-bench llama-bench [.] void (anonymous namespace)::tinyBLAS<8, float __vector(8), float __vector(8), unsigned short, float, float>::gemm<4, 3>(long, long, long, long) ▒ 1.29% llama-bench llama-bench [.] ggml_vec_dot_f16 ▒ 0.56% llama-bench llama-bench [.] ggml_compute_forward_soft_max ▒ 0.56% llama-bench llama-bench [.] ggml_compute_forward_mul_mat ``` gg 作者的 `vec_dot()` 向量化方法目前最佳化 token 生成的速度，但是在 llam.cpp 中還是遇到瓶頸，那便是消化輸入的 token 的速度。畢竟在要總結大量文字或是在 LLaVa 消化圖片時還是需要依靠一個完整的矩陣相乘。這也是為什麼現在依舊有許多外部 BLAS Library。可以觀察到這邊運用到很多像是之前的 `vec_dot_q8_0_q8_0()` 的函式，其實兩者的實作方式有些許相似，以下解釋 sgemm 是怎麼運作的。 ### sgemm.cpp 這個矩陣模組的核心是 `mnpack()` 這一個函式，其作用是將矩陣分解成更小的塊。依據這些分解的小區塊再進入 `gemm` 去運算。 ```cpp class tinyBLAS { public: tinyBLAS(int k, const TA *A, int lda, const TB *B, int ldb, TC *C, int ldc, int ith, int nth) : A(A), B(B), C(C), k(k), lda(lda), ldb(ldb), ldc(ldc), ith(ith), nth(nth) { } void matmul(int m, int n, int task) { if (task == GGML_TASK_TYPE_COMPUTE) mnpack(0, m, 0, n); } private: NOINLINE void mnpack(int m0, int m, int n0, int n) { int mc, nc, mp, np; if (m - m0 <= 0 || n - n0 <= 0) return; if (VECTOR_REGISTERS >= 32 && n - n0 >= 5 && m - m0 >= 5) { mc = 5; nc = 5; gemm5x5(m0, m, n0, n); } else if (n - n0 >= 4 && m - m0 >= 3) { mc = 3; nc = 4; gemm3x4(m0, m, n0, n); } else if (n - n0 >= 4) { mc = 1; nc = 4; gemm1x4(m0, m, n0, n); } else if (m - m0 >= 4) { mc = 4; nc = 1; gemm4x1(m0, m, n0, n); } else { mc = 1; nc = 1; gemm1x1(m0, m, n0, n); } mp = m0 + (m - m0) / mc * mc; np = n0 + (n - n0) / nc * nc; mnpack(mp, m, n0, np); mnpack(m0, mp, np, n); mnpack(mp, m, np, n); } ``` 這邊 `m` 代表行數 `n` 代表列數。在 `mnpack()` 之中會使用遞迴方式分配每一個矩陣計算到它相對應 `llmmnxn` 去做計算。（以上的程式是第一次合併進 llama.cpp 的版本 [commit 8cc91dc](https://github.com/ggerganov/llama.cpp/commit/8cc91dc63c0df397d644a581b2cbeea74eb51ae0) 中的版本，後續不斷有新的 `llmmnxn` 不斷更新）。那這樣的分配的確增加原本一個一個 `block_qn_0` 去做計算更快。它做的更動還有加入更多的 SIMD 指令，若是以演算法的角度去看 `vec_dot()` 也是基本的 $O(n^3)$ 做計算： ```c template <typename T> void LLMM(int m, int n, int k, const T *A, int lda, const T *B, int ldb, T *C, int ldc) { for (int i = 0; i < m; ++i) for (int j = 0; j < n; ++j) { T d = 0; for (int l = 0; l < k; ++l) d += A[lda * i + l] * B[ldb * j + l]; C[ldc * j + i] = d; } } ``` 在 llam.cpp 中使用的技巧是將最裡面的迴圈使用 SIMD 指令展開。 ```c for (int i = 0; i < m; ++i) for (int j = 0; j < n; ++j) { __m256 c0 = _mm256_setzero_ps(); __m256 c1 = _mm256_setzero_ps(); for (int l = 0; l < k; l += 16) { c0 = _mm256_fmadd_ps(_mm256_loadu_ps(A + lda * i + l + 0), _mm256_loadu_ps(B + ldb * j + l + 0), c0); c1 = _mm256_fmadd_ps(_mm256_loadu_ps(A + lda * i + l + 8), _mm256_loadu_ps(B + ldb * j + l + 8), c1); } C[ldc * j + i] = hsum(c0) + hsum(c1); } ``` 但實際上效能提升有限，因為編譯器在最佳化時，即使原本沒有 SIMD 的程式，也會做類似的處理。而在 sgemm.cpp 中，中間層的迴圈使用 SIMD 展開。以 `gemm1x4` 為範例，直接將三層迴圈中的中下兩層直接使用 SIMD 展開： ```c NOINLINE void gemm1x4(int m0, int m, int n0, int n) { BEGIN_KERNEL(1, 4) __m256 c0 = _mm256_setzero_ps(); __m256 c1 = _mm256_setzero_ps(); __m256 c2 = _mm256_setzero_ps(); __m256 c3 = _mm256_setzero_ps(); const TB *Bp0 = B + ldb * (j + 0); const TB *Bp1 = B + ldb * (j + 1); const TB *Bp2 = B + ldb * (j + 2); const TB *Bp3 = B + ldb * (j + 3); const TA *Ap = A + lda * i; for (int l = 0; l < k; ++l) { float da0 = unhalf(Ap[l].d); __m256i f = load(Ap + l); __m256i u = _mm256_sign_epi8(f, f); // make all the number >=0 __m256 d0 = _mm256_set1_ps(unhalf(Bp0[l].d) * da0); // b.scalar * a.scalar __m256 d1 = _mm256_set1_ps(unhalf(Bp1[l].d) * da0); // make it a vector __m256 d2 = _mm256_set1_ps(unhalf(Bp2[l].d) * da0); __m256 d3 = _mm256_set1_ps(unhalf(Bp3[l].d) * da0); __m256 g0 = updot(u, _mm256_sign_epi8(load(Bp0 + l), f)); __m256 g1 = updot(u, _mm256_sign_epi8(load(Bp1 + l), f)); __m256 g2 = updot(u, _mm256_sign_epi8(load(Bp2 + l), f)); __m256 g3 = updot(u, _mm256_sign_epi8(load(Bp3 + l), f)); c0 = madd(d0, g0, c0); c1 = madd(d1, g1, c1); c2 = madd(d2, g2, c2); c3 = madd(d3, g3, c3); } C[ldc * (j + 0) + i] = hsum(c0); C[ldc * (j + 1) + i] = hsum(c1); C[ldc * (j + 2) + i] = hsum(c2); C[ldc * (j + 3) + i] = hsum(c3); END_KERNEL() } ``` 假設在 `Q8_0` 中，以 usigned int8 儲存，一次計算 256 個位元，會處理 256 / 8 = 32 筆資料，而 gemm1x4 就會處理 32 * 128 大小的矩陣。以此類推，尚有 `gemm4x3()`，`gemm4x1()`，`gemm1x4()`，`gemm1x1()`。可以說 sgemm 便是將大矩陣不斷分割成小矩陣進行 SIMD 並行指令計算。這也是得益於大型語言模型穩定的矩陣大小，讓針對性的此專案的矩陣相乘得以發展。 ## TODO: 選定 data type 並聚焦在 matmul 實作 > w/ SIMD, OpenMP, oneMKL, OpenBLAS, …) 的效益 => 我們需要有個 wrapper/tester，得以涵蓋現有實作並整合進 llama.cpp，在 Linux 平台驗證，perf (搭配第六次作業提到的效能分析方式，排除系統干擾) ### 實驗暫存 `make LLAMA_NO_ACCELERATE=1` tiny-llama q4_0 | model | size | params | backend | threads | test | t/s | | ------------- | ----------:| ------:| ------- | -------:| -----:| --------------:| | llama 1B Q4_0 | 606.53 MiB | 1.10 B | CPU | 8 | pp512 | 235.26 ± 13.62 | | llama 1B Q4_0 | 606.53 MiB | 1.10 B | CPU | 8 | tg16 | 43.79 ± 4.99 | ### SIMD ### OpenMP OpenMP 在 llama.cpp 的實作在於矩陣相乘 backend 使用 [BLIS](https://github.com/flame/blis) 這一個專案中在編譯 BLIS 的時候使用 `make LLAMA_BLIS=1` 可以順利編譯，但是使用 `./llama-bench` 出現 `./llama-bench: error while loading shared libraries: libblis.so.4: cannot open shared object file: No such file or directory` ### oneMKL ### OpenBLAS `make LLAMA_OPENBLAS=1` tiny-llama q4_0 | model | size | params | backend | threads | test | t/s | | ------------- | ----------:| ------:| ------- | -------:| -----:| ------------:| | tinyllama 1B Q4_0 | 606.53 MiB | 1.10 B | BLAS | 8 | pp512 | 72.94 ± 0.08 | | tinyllama 1B Q4_0 | 606.53 MiB | 1.10 B | BLAS | 8 | tg16 | 50.04 ± 0.89 | ### CUDA 這邊順便附上在 NVIDIA GeForce GTX 1070 的實驗，使用 `make LLAMA_CUDA=1` 去編譯，會出現 error，跟著他的步驟去 Nvidia 官往查詢顯示卡的計算能力 https://developer.nvidia.com/cuda-gpus ，這邊看到 GTX 1070 的計算能力為 6.1，便使用命令 `export CUDA_DOCKER_ARCH=compute_61`，即可完成編譯。 `ggml_cuda_init: GGML_CUDA_FORCE_MMQ: no ggml_cuda_init: CUDA_USE_TENSOR_CORES: yes ggml_cuda_init: found 1 CUDA devices: Device 0: NVIDIA GeForce GTX 1070, compute capability 6.1, VMM: yes` | model | size | params | backend | test | t/s | | ------------- | ----------:| ------:| ------- | -----:| --------------:| | tinyllama 1B Q4_0 | 606.53 MiB | 1.10 B | CUDA | pp512 | 505.89 ± 77.20 | | tinyllama 1B Q4_0 | 606.53 MiB | 1.10 B | CUDA | tg16 | 24.76 ± 0.33 | ## TODO: 提出更快的 matmul > 搭配作業系統核心層級的調整，讓運算能佔用更多 CPU 資源 ## 數學分析探討量化過程的精確度損失令 n 為原始數值，max 為該區塊的最大值以`q4_0` 為例：整個量化與解量化所產生的誤差為： $n-((n/\frac{max}{-8}+8.5)-8)*\frac{max}{-8}$ $=n-(\frac{-8}{max}+8.5-8)*\frac{max}{-8}$ $=n-(\frac{-8}{max}+0.5)*\frac{max}{-8}$ $=n-(n+\frac{max}{-8})$ $=\frac{max}{-8}$ 在 `q4_0 * q4_0` 時令前項為 $n_1$、block 中的最大值為 $max_1$，後項為 $n_2$、最大值為 $max_2$ 其相乘所帶來的誤差為： $n_1*n_2-((n_1/\frac{max_1}{-8}+8.5)-8)*((n_2/\frac{max_2}{-8}+8.5)-8)*(\frac{max_1}{-8}*\frac{max_2}{-8})$ $=n_1*n_2-(\frac{-8n_1}{max_1}+0.5)*(\frac{-8n_2}{max_2}+0.5)*(\frac{max_1}{-8}*\frac{max_2}{-8})$ $=n_1*n_2-(\frac{64n_1n_2}{max1max2}+\frac{-4n_1}{max_1}+\frac{-4n_2}{max_2}+0.25)*(\frac{max_1}{-8}*\frac{max_2}{-8})$ $=n_1*n_2-n_1*n_2+\frac{max_2*n_1}{16}+\frac{max_1*n_2}{16}+\frac{max_1*max_2}{64}$ $=\frac{max_1*n_2}{16}+\frac{max_2*n_1}{16}+\frac{max_1*max_2}{64}$ 在 `q4_0 * q8_0` 時，誤差為：令前項為 $n_1$、block 中的最大值為 $max_1$，後項為 $n_2$、最大值為 $max_2$ $n_1*n_2-((n_1/\frac{max_1}{-8}+8.5)-8)(round(n_2/\frac{max}{127}))*(\frac{max_1}{-8}*\frac{max_2}{127})$ $=n_1*n_2-(\frac{-8n_1}{max_1}+0.5)*(round(\frac{127n_2}{max_2}))*(\frac{max_1}{-8}*\frac{max_2}{127})$ $=n_1*n_2-(\frac{-8n_1}{max_1}+0.5)*(\lfloor\frac{127n_2}{max_2}+0.5\rfloor)*(\frac{max_1}{-8}*\frac{max_2}{127})$ $=n_1*n_2-(\frac{max_1}{-8}*(\frac{-8n_1}{max_1}+0.5))*(\frac{max_2}{127}*\lfloor\frac{127n_2}{max_2}+0.5\rfloor)$ $=n_1*n_2-(n_1+\frac{max_1}{-16})*(\frac{max_2}{127}*(\frac{127n_2}{max_2}\pm1))$ 使用 $\pm1$ 表示 $n_2>0,n_2<0$ 的最大修正 $=n_1*n_2-(n_1+\frac{max_1}{-16})*((n_2\pm\frac{max_2}{127}))$ $=n_1*n_2-n_1*n_2\pm\frac{max_2*n_1}{127}-\frac{max_1*n_2}{16}\pm\frac{max_1*max_2}{-2032}$ $=\pm\frac{max_2*n_1}{127}-\frac{max_1*n_2}{16}\pm\frac{max_1*max_2}{-2032}$ 使用代數 $k\leq\frac{127n_2}{max_2}+0.5<k+1$ ? 為何 `q8_0` 做量化使用 Symmetric quantize，而不是跟 `q4_0` 一樣使用 Affine Quantizer？ ## 參考 [#909](https://github.com/ggerganov/llama.cpp/issues/909) [#896](https://github.com/ggerganov/llama.cpp/pull/896)