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

執行人: vestata
專題解說影片

Reviewed by yu-hsiennn

整個程式運作都的計算熱點會是矩陣乘法

指的是 整個程式運作的計算熱點都會是矩陣乘法? 還有 熱點 指的是什麼?

Wikipedia: hot spot 我翻譯為「熱點」，解釋為程式執行過程中大部分時間消耗的區域，〈LLaMA Now Goes Faster on CPUs〉中使用的描述是 Wikipedia: Bottleneck 但是我覺得在這邊使用熱點比較符合。

Reviewed by yenslife

既然知道 ((1 << 3) - 1) 表示 7，那為什麼一開始不要直接寫 7 就好？我剛剛在猜會不會與操作成本有關，但看到之後開發者修改成 -8 就更不懂了。

使用 ((1 << 3) - 1) 表示為 7 而不是直接使用 7 的具體原因…我還不知道，從第一個版本 commit 26c0846 便是使用 ((1 << 3) - 1) 延用到現在，我認為是為了增加程式可解讀性。若有發現我會再更新。
可以參考我之後的在 Quantization 下面的補充，其中

$\mathrm{\Delta }$ 是其步長，一般會定義為

$\mathrm{\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 納入考量，這是針對現代 x86 處理器的 AVX 和 FMA 指令集擴展開發的快速矩陣乘法實作，程式碼精簡又性能表現優異。參照該專案的介紹文章 Beating NumPy's matrix multiplication in 150 lines of C code。

任務簡介

探討 llama.cpp 的效能表現，分析其效能瓶頸並嘗試改進

TODO: 讀懂〈LLaMA Now Goes Faster on CPUs〉並重現實驗

不同 quantization 所產生的精度差異，如何實驗它對模型的影響？

改進你的漢語表達。

llamafile 是為了簡化使用 LLM 所生的專案。結合了 llama.cpp 與 Cosmopolitan Libc 使配置與運行壓縮到單一可執行檔中。

注意用語：

• file 是「檔案」，而非「文件」(document)
• thread 是「執行緒」，而非「線程」

好的，收到

Quantization

〈LLaMA Now Goes Faster on CPUs〉這篇文章是中使用不同的量化手法進行效能的比較，主要針對 prompt (tok/sec) 即讀取 token 的速度，與 eval (tok/sec) 輸出的 token 速度進行比較。一般訓練好並開放的大型語言模型權重 .pth 是使用 BF16 的資料型態，在 llama.cpp 需要將 .pth 權重轉成 .gguf 的資料型態。

.gguf 是 llama.cpp 的作者所開發的資料格式，可以使用在其相關的專案中，如以下圖例：

https://huggingface.co/docs/hub/gguf

TODO
gguf/tensor 的差異
補 gguf 的 mmap()

依據〈Quantizing deep convolutional networks for efficient inference: A whitepaper〉 這篇論文，可以將量化的分為兩種 Affine Quantizer 與 Symmetric Quantizer 兩種。

Symmetric quantize

$$x_{\text{int}}=\text{round}\left(\frac{x}{\mathrm{\Delta }}\right)$$

$$x_Q=\text{clamp}(-\frac{N_{\text{levels}}}{2},\frac{N_{\text{levels}}}{2}-1,x_{\text{int}})\text{if signed}$$

$$x_Q=\text{clamp}(0,N_{\text{levels}}-1,x_{\text{int}})\text{if un-signed}$$

其中

$$\text{clamp}(a,b,x)=\{\begin{array}{ll}a& \text{if }x\le a\\ x& \text{if }a\le x\le b\\ b& \text{if }x\ge b\end{array}$$

s 為scale factor，z 為 zero point

此方法會對應到等等的 Q8_0 手法。

Affine Quantizer

量化：

$$x_{\text{int}}=\text{round}\left(\frac{x}{\mathrm{\Delta }}\right)+z$$

$$x_Q=\text{clamp}(0,N_{\text{levels}}-1,x_{\text{int}})$$

解量化：

$x_{float}=(x_Q-z)\mathrm{\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

TODO

Q4_0 會儲存在一個 block_q4_0 結構體中，其大小會是 32，其中 d 是他的縮放比例（scalar），qs 是一個 uint8_t 陣列用來儲存兩個 4 位元的，在量化的過程中，原本 BF16 的模型權重以 block_q4_0 為單位去做儲存：

#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 的變更，將 ((1 << 3) - 1) 改為 -8 因為這邊是為了將原本 FP16 量化到 4 個位元，所以 4 位元的二的補數表示範圍為 -8 到 7。為了善用更大的分母所以使用 -8，其所導致的整個存放在 .gguf 的權重變成負數，在 perplexity 的表現上勝於 ((1 << 3) - 1) 的方法，見 #729。

而在後面要將其放進 4 位元之中，將原本權重參數乘上縮放比例 id，但是後面的順序的擺放是將 x[i*qk + 0 + j] 與 x[i*qk + qk/2 + j] 擺在一起。這是在 commit b9fd7ee 的變更。

         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=A{B}^T\iff C=B{A}^T$ ，這與傳統的矩陣相乘記憶體儲存手法不同。

README.md

在這邊補充一下，一般我們在做矩陣相乘的時候

$C=AB$ 會遇到的狀態會是

$A,B$ 在記憶體的排列是為 row-major。當處理

$A$ 時，可以利用其連續記憶體的特性，減少到主記憶體的成本，這樣能夠利用 L1 cache。然而，

$B$ 的矩陣乘法需要每一 column 的資料，這與 row-major 記憶體排列剛好相反，所以

$C=AB$ 這個乘法最大的效能瓶頸會是

$B$ 對主記憶體成本。有很多矩陣乘法手段可以更好的利用 L1 cache，一個簡單的方法是對

$B$ 做轉置，詳細可以參考我對矩陣相乘的筆記。

所以對於 quantize_row_q4_0_reference() 的變更就是為了符合在

$C=B{A}^T$ 的運算規則之下。對其中一個矩陣雖然因為 L1 cache 的優勢，但是以演算法的角度尚需要

$O(n^3)$ 的成本。因此，llama.cpp 原本的作法是將量化後的 uint 以相鄰的方式排列，但之後需要額外的轉置成本。改動後，量化時便將其以矩陣轉置

$A^T$ 的排序放進記憶體中，以減少後續的轉置成本。

Perplexity

經過量化數值的精確度一定會喪失，只是不同的量化手段會有不同精確度喪失，但是這精確度喪失會大型語言模型會有什麼影響呢？回到模型的根本，是使用模型權重去選擇生成下一個最有可能出現的 token，我們可以針對生出的 token 進行評估。在 llama.cpp 專案之中的方法是使用 perplexity。其定義如下：

$$\text{PPL}(X)=\exp \{-\frac{1}{t}\sum _i^t\log p_{\theta }(x_i\mid x_{<i})\}$$

闡述其特性、考量因素，以及對應的限制。

給定特定的文章，這邊使用的是 WikiText-2 在給定的文章中，模型使用

$X_{<i}$ 的已知 token 去推斷

$X_i$ 的 token，

$t$ 表示為 token 數，所以 perplexity 越低，表示模型的預測結果與實際 token 序列越接近。

文章中使用不同的 quantization 手法，有 Q4_0，Q8_0，FP16。以 Q4_0 為例，以下是主要實作的程式：

重現實驗

使用的測試工具是 llama-bench，這是在專案之中被廣泛的使用的測試工具，分別會計算輸入 token(prompt) 與輸出 token(eval) 的時間。主要的評估函式是 test_promt() 與 test_gen() 。其中 test_promt() 是以 rand() 去生成隨機數字，計算模型在消化 n_prompt 的隨機數字需要多少時間。

    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 次所花費的時間。

    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;
    }

測試設備：

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 使用舊版向量內積，與 commit 74f33ad，此版本是 2024-5-23 是我拿來做測試的版本，使用的是 sgemm 新的矩陣乘法模組（llamafile 的乘法模組）。我將以這兩個做實驗比較。在 x86 機器中使用 make 即可完成編譯，其中執行緒的數量皆是使用預設。

補上對應的超連結。

已補上

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) 的表現上則相近。

討論其精確度的犧牲狀況。

TODO: 影響 LLaMA 推理速度的因素

找出其中效能瓶頸的地方，可能是矩陣運算，或是記憶體 (特別是 data model)，予以量化，善用 perf 和 Ftrace 等工具

矩陣相乘

使用 perf 依照 cpu-cycles，faults，cache-references，cache-misses 進行比較

$ 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 的 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=B{A}^T$ 並且在量化的時候便已經將

$A$ 向量化排列在連續記憶體之中。有很多開源的 BLAS libraries 有很純熟的矩陣相乘手法，但是放在已經向量化的場景反而會降低速度，因為這些 BLAS libraries 有很多複雜的手段會增加許多延遲，是對應各層面的矩陣相乘場景，而 gg 作者的向量內積是為了這個專案特別設計的。

其中 ggml_vec_dot_q4_0_q8_0() 是實作整個向量內積的底層程式。但是細看程式碼，先以矩陣相乘的演算法的角度（SIMD 等技巧先忽略），我們可以發現這是一個 q4_0 乘上 q8_0 的混合精度乘法：

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() 這個函式。

typedef struct {
    ggml_half d;           // delta
    uint8_t qs[QK4_0 / 2]; // nibbles / quants
} block_q4_0;

	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 的量化。

    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()，降低了量化過程中所帶來的精度喪失。在之後底層運算的實作皆是使用這個向量乘法演算法。

建立數學模型，探討量化過程的精確度損失和補償方案。

數學分析！

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 核心的排程會降低計算效能。

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 然後移到排程佇列最末端，讓其他執行緒優先執行。

搭配 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.

if (atomic_fetch_sub(&state->shared->n_active, 1) == 1)

所以如果 atomic_fetch_sub() 回傳是 1 的話，一方面執行 n_active 執行緒數減一，一方面表示這是矩陣乘法多執行緒中的最後一個執行緒，其他工作的執行緒已完成，最後一個執行緒會有一個收尾的動作。以下這一段是處理收尾的主要程式：

            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 測試發現這邊預設的執行緒是 8。後來發現這邊也有緣由。在原本的版本使用 get_num_physical_cores() 去判斷現在的設備有多少邏輯核心，如果邏輯核心數超過 4 個便會使用 邏輯核心數 / 2 個執行緒。

以下是我的設備為例有 24 核，所以在原本向量矩陣相乘的方法將 12 配置給執行中的程式。

而 sgemm.cpp 改成 8 是因為 13th Gen Intel(R) Core(TM) i7-13700 其中 24 個核心包含 8 個雙執行緒的效能核心與 8 個單執行緒的效率核心，希望避免將大量的矩陣計算移到效率核心去做計算。
執行緒選取的更動 這邊有對核心的選取，我對這邊核心相關的程式細節不了解(cpuid(), is_hybrid_cpu(), is_running_on_efficiency_core())，希望知道人可以補充。這邊使用 pthread_setaffinity_np()
將設定在效能核心上。

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.
man page

矩陣相乘

llama.cpp 支援許多種 backend，也就是多種矩陣相乘的計算模組，原有支援 CPU，OpenCL，Vulkan，Kompute，Metal，GPU BLAS，BLAS。在 commit 8cc91dc 加入 sgemm.cpp 這個乘法模組。文章中有提到其在 Q8 與 FP16 有最好的效果，那可以從其計算熱點觀察。

  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 去運算。

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 中的版本，後續不斷有新的 llmmnxn 不斷更新）。

那這樣的分配的確增加原本一個一個 block_qn_0 去做計算更快。它做的更動還有加入更多的 SIMD 指令，若是以演算法的角度去看 vec_dot() 也是基本的

$O(n^3)$ 做計算：

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 指令展開。

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 展開：

    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 這一個專案中
在編譯 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)\ast \frac{max}{-8}$

$=n-(\frac{-8}{max}+8.5-8)\ast \frac{max}{-8}$

$=n-(\frac{-8}{max}+0.5)\ast \frac{max}{-8}$

$=n-(n+\frac{max}{-8})$

$=\frac{max}{-8}$

在 q4_0 * q4_0 時
令前項為

$n_1$、block 中的最大值為

$ma{x}_1$，後項為

$n_2$、最大值為

$ma{x}_2$
其相乘所帶來的誤差為：

$n_1\ast n_2-((n_1/\frac{ma{x}_1}{-8}+8.5)-8)\ast ((n_2/\frac{ma{x}_2}{-8}+8.5)-8)\ast (\frac{ma{x}_1}{-8}\ast \frac{ma{x}_2}{-8})$

$=n_1\ast n_2-(\frac{-8n_1}{ma{x}_1}+0.5)\ast (\frac{-8n_2}{ma{x}_2}+0.5)\ast (\frac{ma{x}_1}{-8}\ast \frac{ma{x}_2}{-8})$

$=n_1\ast n_2-(\frac{64n_1{n}_2}{max1max2}+\frac{-4n_1}{ma{x}_1}+\frac{-4n_2}{ma{x}_2}+0.25)\ast (\frac{ma{x}_1}{-8}\ast \frac{ma{x}_2}{-8})$

$=n_1\ast n_2-n_1\ast n_2+\frac{ma{x}_2\ast n_1}{16}+\frac{ma{x}_1\ast n_2}{16}+\frac{ma{x}_1\ast ma{x}_2}{64}$

$=\frac{ma{x}_1\ast n_2}{16}+\frac{ma{x}_2\ast n_1}{16}+\frac{ma{x}_1\ast ma{x}_2}{64}$

在 q4_0 * q8_0 時，誤差為：
令前項為

$n_1$、block 中的最大值為

$ma{x}_1$，後項為

$n_2$、最大值為

$ma{x}_2$

$n_1\ast n_2-((n_1/\frac{ma{x}_1}{-8}+8.5)-8)(round(n_2/\frac{max}{127}))\ast (\frac{ma{x}_1}{-8}\ast \frac{ma{x}_2}{127})$

$=n_1\ast n_2-(\frac{-8n_1}{ma{x}_1}+0.5)\ast (round(\frac{127n_2}{ma{x}_2}))\ast (\frac{ma{x}_1}{-8}\ast \frac{ma{x}_2}{127})$

$=n_1\ast n_2-(\frac{-8n_1}{ma{x}_1}+0.5)\ast (\lfloor \frac{127n_2}{ma{x}_2}+0.5\rfloor )\ast (\frac{ma{x}_1}{-8}\ast \frac{ma{x}_2}{127})$

$=n_1\ast n_2-(\frac{ma{x}_1}{-8}\ast (\frac{-8n_1}{ma{x}_1}+0.5))\ast (\frac{ma{x}_2}{127}\ast \lfloor \frac{127n_2}{ma{x}_2}+0.5\rfloor )$

$=n_1\ast n_2-(n_1+\frac{ma{x}_1}{-16})\ast (\frac{ma{x}_2}{127}\ast (\frac{127n_2}{ma{x}_2}\pm 1))$ 使用

$\pm 1$ 表示

$n_2>0,n_2<0$ 的最大修正

$=n_1\ast n_2-(n_1+\frac{ma{x}_1}{-16})\ast ((n_2\pm \frac{ma{x}_2}{127}))$

$=n_1\ast n_2-n_1\ast n_2\pm \frac{ma{x}_2\ast n_1}{127}-\frac{ma{x}_1\ast n_2}{16}\pm \frac{ma{x}_1\ast ma{x}_2}{-2032}$

$=\pm \frac{ma{x}_2\ast n_1}{127}-\frac{ma{x}_1\ast n_2}{16}\pm \frac{ma{x}_1\ast ma{x}_2}{-2032}$

使用代數

$k\le \frac{127n_2}{ma{x}_2}+0.5<k+1$ ?

為何 q8_0 做量化使用 Symmetric quantize，而不是跟 q4_0 一樣使用 Affine Quantizer？

參考

#909
#896