Try   HackMD

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 下面的補充,其中

Δ 是其步長,一般會定義為
Δ=2(b1)1
,在 q4_0 也就是 ((1 << 3) - 1) ,但是因為使用二補數,所以 4 位元表示範圍應該是 -8~7。我們現在令 n = 30000, max = 65504(FP16 表示的最大值) 作為舉例,分別進行 -8 與 7 的量化與解量化為
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.cppCosmopolitan 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 的資料型態。

.ggufllama.cpp 的作者所開發的資料格式,可以使用在其相關的專案中,如以下圖例:

image

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

xint=round(xΔ)

xQ=clamp(Nlevels2,Nlevels21,xint)if signed

xQ=clamp(0,Nlevels1,xint)if un-signed

其中

clamp(a,b,x)={aif xaxif axbbif xb

s 為scale factor,z 為 zero point

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

Affine Quantizer

量化:

xint=round(xΔ)+z

xQ=clamp(0,Nlevels1,xint)

解量化:

xfloat=(xQz)Δ

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

量化手法介紹

手法通常分為 Qn_0Qn_1Qn_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 的開發規則,在此專案之中的定義的矩陣相乘是

CT=ABTC=BAT ,這與傳統的矩陣相乘記憶體儲存手法不同。
image

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=BAT 的運算規則之下。對其中一個矩陣雖然因為 L1 cache 的優勢,但是以演算法的角度尚需要
O(n3)
的成本。因此,llama.cpp 原本的作法是將量化後的 uint 以相鄰的方式排列,但之後需要額外的轉置成本。改動後,量化時便將其以矩陣轉置
AT
的排序放進記憶體中,以減少後續的轉置成本。

Perplexity

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

PPL(X)=exp{1titlogpθ(xix<i)}

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

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

X<i 的已知 token 去推斷
Xi
的 token,
t
表示為 token 數,所以 perplexity 越低,表示模型的預測結果與實際 token 序列越接近。
image

文章中使用不同的 quantization 手法,有 Q4_0Q8_0FP16。以 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_0FP16 使用 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_0q8_0FP16 使用 llamafile 加速後,在 prompt (tok/sec) 上的性能分別提升了 440%、%、413% 和 364%,有顯著的進步。,在 eval (tok/sec) 的表現上則相近。

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

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

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

矩陣相乘

使用 perf 依照 cpu-cyclesfaultscache-referencescache-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=BAT 並且在量化的時候便已經將
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_INITGGML_TASK_TYPE_COMPUTEGGML_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,也就是多種矩陣相乘的計算模組,原有支援 CPUOpenCLVulkanKomputeMetalGPU BLASBLAS。在 commit 8cc91dc 加入 sgemm.cpp 這個乘法模組。文章中有提到其在 Q8FP16 有最好的效果,那可以從其計算熱點觀察。

  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(n3) 做計算:

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/max8+8.5)8)max8
=n(8max+8.58)max8

=n(8max+0.5)max8

=n(n+max8)

=max8

q4_0 * q4_0
令前項為

n1、block 中的最大值為
max1
,後項為
n2
、最大值為
max2

其相乘所帶來的誤差為:
n1n2((n1/max18+8.5)8)((n2/max28+8.5)8)(max18max28)

=n1n2(8n1max1+0.5)(8n2max2+0.5)(max18max28)

=n1n2(64n1n2max1max2+4n1max1+4n2max2+0.25)(max18max28)

=n1n2n1n2+max2n116+max1n216+max1max264

=max1n216+max2n116+max1max264

q4_0 * q8_0 時,誤差為:
令前項為

n1、block 中的最大值為
max1
,後項為
n2
、最大值為
max2

n1n2((n1/max18+8.5)8)(round(n2/max127))(max18max2127)

=n1n2(8n1max1+0.5)(round(127n2max2))(max18max2127)

=n1n2(8n1max1+0.5)(127n2max2+0.5)(max18max2127)

=n1n2(max18(8n1max1+0.5))(max2127127n2max2+0.5)

=n1n2(n1+max116)(max2127(127n2max2±1))
使用
±1
表示
n2>0,n2<0
的最大修正
=n1n2(n1+max116)((n2±max2127))

=n1n2n1n2±max2n1127max1n216±max1max22032

=±max2n1127max1n216±max1max22032

使用代數

k127n2max2+0.5<k+1 ?

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

參考

#909
#896