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

執行人: LIAO-JIAN-PENG
專題解說錄影

Reviewed by Petakuo

能否解釋為何逐字處理輸入序列會容易遇到梯度消失和爆炸，以及控制門 (Gate) 具體解決此問題的運作方式?

LIAO-JIAN-PENG
已附上說明

Reviewed by marvin0102

transformer 在字串序列增加時，矩陣運算也會增加，運算效率與CNN RNN 相比，優勢在哪裡，可否舉例說明？
LIAO-JIAN-PENG
已附上說明

Reviewed by vestata

為什麼在 matmul 的 loop tiling 選擇 x=16？ loop tiling 大小是如何選擇的呢？

主流的 CPU Cache 的 Cache Line 大小多為 64 Bytes
一個 float 是 4 bytes
因此一個 cache line 可以存放 16 個 float 數字 (64/4 = 16)
設定 x = 16 可以最大化 cache line 的利用率，減少 cache miss
參考：現代處理器設計: Cache 原理和實際影響

Reviewed by jserv

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

任務簡介

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

TODO: 探討提升矩陣乘法的手法

矩陣乘法演算法

修正圖片的存取權限。

注意用語：

• algorithm 是「演算法」，而非「算法」(the way to calculate)

Matrix multiplication algorithm wiki

Strassen 矩陣乘法

Strassen 矩陣乘法是一種高效的矩陣乘法演算法，它比傳統的矩陣乘法具有更低的時間複雜度

傳統的矩陣乘法對於兩個

$n\times n$ 矩陣 A 和 B，需要

$O(n^3)$ 的時間。Strassen 算法則透過將矩陣劃分成較小的子矩陣，並用一種特殊的計算方式，將時間複雜度降至約

$O(n^{2.81})$。

在

$2\times 2$ 矩陣相乘中，傳統矩陣相乘使用 8 個乘法，4 個加法，在 Strassen 算法中，使用到 7 個乘法 18 加法，因為乘法跟加法的運算複雜度有差，當運算矩陣大型矩陣乘法時，乘法的運算成本也就會體現出來。

雖然說一般矩陣會是剛好

$2\times 2$ 的矩陣，但是我們可以使用 Block matrix 的方式將矩陣從

$n\times n$ 變成

$2\times 2$ 的矩陣。

資料平行處理

進行效能比較和實作分析，從而歸納出提升矩陣乘法的手法

matmul

Implementation	Long description
Naive	Most obvious implementation
Transposed	Transposing the second matrix for cache efficiency
sdot w/o hints	Replacing the inner loop with BLAS sdot()
sdot with hints	sdot() with a bit unrolled loop
SSE sdot	vectorized sdot() with explicit SSE instructions
SSE+tiling	sdot SSE sdot() with loop tiling
OpenBLAS sdot	sdot() provided by OpenBLAS
OpenBLAS sgemm	sgemm() provided by OpenBLAS

sdot w/o hints: 處理兩個向量相乘

float sdot_1(int n, const float *x, const float *y)
{
    int i;
    float s = 0.0f;
    for (i = 0; i < n; ++i) s += x[i] * y[i];
    return s;
}

在一些數值運算庫（如 BLAS，Basic Linear Algebra Subprograms），sdot 函數用於計算兩個向量的內積。例如，給定兩個向量

$x$ 和

$y$，sdot 會返回以下結果：

$$sdot(x,y)=\sum _{i=1}^n{x}_i\times y_i$$

sdot with hints: 用 loop unrolling 展開

float sdot_8(int n, const float *x, const float *y)
{
    int i, n8 = n>>3<<3;
    float s, t[8];
    t[0] = t[1] = t[2] = t[3] = t[4] = t[5] = t[6] = t[7] = 0.0f;
    for (i = 0; i < n8; i += 8) {
        t[0] += x[i+0] * y[i+0];
        t[1] += x[i+1] * y[i+1];
        t[2] += x[i+2] * y[i+2];
        t[3] += x[i+3] * y[i+3];
        t[4] += x[i+4] * y[i+4];
        t[5] += x[i+5] * y[i+5];
        t[6] += x[i+6] * y[i+6];
        t[7] += x[i+7] * y[i+7];
    }
    for (s = 0.0f; i < n; ++i) s += x[i] * y[i];
    s += t[0] + t[1] + t[2] + t[3] + t[4] + t[5] + t[6] + t[7];
    return s;
}

透過 Loop unrolling 的技巧，編譯器或組譯器可以預先計算出迴圈中每次反覆執行所需要的記憶體位置，並直接參照到矩陣的變數位置。這樣可以減少執行時的位址計算開銷，產生較高效能的機器碼。不過，Loop unrolling 也可能會導致程式碼大小增加，因此編譯器會根據情況決定是否應用這個技術。

避免濫用「優化」ㄧ詞，務必使用課程教材規範的術語。參見：

• 資訊科技詞彙翻譯
• 你所不知道的 C 語言：編譯器和最佳化原理篇

編譯器在最佳化階段會自動嘗試各種最佳化手法，包括 Loop unrolling，來提升程式的執行效能。然而，編譯器是否真正採用 Loop unrolling，端視迴圈的大小、複雜度、硬體特性等多種因素而定。對於一些簡單的迴圈，編譯器較有可能進行展開，但對於大型迴圈，展開可能會導致程式碼膨脹，因此編譯器可能不會採用此最佳化手法。

Loop unrolling wiki

SSE sdot: 使用 SSE

注意用語：

• data 是「資料」，而非「數據」

SSE(Streaming SIMD Extensions)是 Intel 處理器的一組指令集擴展，該指令集對於處理大量資料在相同操作(例如多媒體處理、科學計算)下特別有效。

SSE（Streaming SIMD Extensions）最初新增 8 個新的 128 位元暫存器，稱為 XMM0 到 XMM7。

XMM0 through 7 from Streaming SIMD Extensions wiki

float sdot_sse(int n, const float *x, const float *y)
{
    int i, n8 = n>>3<<3;             
    __m128 vs1, vs2;
    float s, t[4];
    vs1 = _mm_setzero_ps();          // {0.f, 0.f, 0.f, 0.f}
    vs2 = _mm_setzero_ps();          // {0.f, 0.f, 0.f, 0.f}
    for (i = 0; i < n8; i += 8) {
        __m128 vx1, vx2, vy1, vy2;
        vx1 = _mm_loadu_ps(&x[i]);   // load {x[0], x[1], x[2], x[3]}
        vx2 = _mm_loadu_ps(&x[i+4]); // load {x[4], x[5], x[6], x[7]}
        vy1 = _mm_loadu_ps(&y[i]);   // load {y[0], y[1], y[2], y[3]}
        vy2 = _mm_loadu_ps(&y[i+4]); // load {y[4], y[5], y[6], y[7]}
        vs1 = _mm_add_ps(vs1, _mm_mul_ps(vx1, vy1)); // vs1 += vx1*vy1
        vs2 = _mm_add_ps(vs2, _mm_mul_ps(vx2, vy2));
    }
    for (s = 0.0f; i < n; ++i) s += x[i] * y[i]; // Complete the remaining (insufficient 8) operations.
    _mm_storeu_ps(t, vs1);
    s += t[0] + t[1] + t[2] + t[3]; // {t[0], t[1], t[2], t[3]}
    _mm_storeu_ps(t, vs2);
    s += t[0] + t[1] + t[2] + t[3]; // {t[4], t[5], t[6], t[7]}
    return s;
}

loop tiling

int i, j, ii, jj, x = 16, n_b_rows = n_a_cols;
float **m, **bT;
m = mat_init(n_a_rows, n_b_cols);
bT = mat_transpose(n_b_rows, n_b_cols, b);
for (i = 0; i < n_a_rows; i += x) {
    for (j = 0; j < n_b_cols; j += x) {
        // loop tiling
        int je = n_b_cols < j + x? n_b_cols : j + x;
        int ie = n_a_rows < i + x? n_a_rows : i + x;
        for (ii = i; ii < ie; ++ii)
            for (jj = j; jj < je; ++jj)
                m[ii][jj] += sdot_sse(n_a_cols, a[ii], bT[jj]);
    }
}

Loop tiling 將大迴圈分成多個小塊 (tile) 迴圈，有效的利用記憶體快取，提升資料在快取中的命中率。

Loop nest optimization

測量不同實作方式的時間

$ lscpu
Architecture:                       x86_64
CPU op-mode(s):                     32-bit, 64-bit
Address sizes:                      39 bits physical, 48 bits virtual
Byte Order:                         Little Endian
CPU(s):                             12
On-line CPU(s) list:                0-11
Vendor ID:                          GenuineIntel
Model name:                         Intel(R) Core(TM) i5-10400 CPU @ 2.90GHz
CPU family:                         6
Model:                              165
Thread(s) per core:                 2
Core(s) per socket:                 6
Socket(s):                          1
Stepping:                           3
CPU max MHz:                        4300.0000
CPU min MHz:                        800.0000
BogoMIPS:                           5799.77
Virtualization:                     VT-x
L1d cache:                          192 KiB (6 instances)
L1i cache:                          192 KiB (6 instances)
L2 cache:                           1.5 MiB (6 instances)
L3 cache:                           12 MiB (1 instance)
NUMA node(s):                       1
NUMA node0 CPU(s):                  0-11

• 1024+i*64(i=0~50)
  
  
  忽略 Naive
  
  

matmul_bench

相較於 matmul 專案的只有單執行緒比較，在 matmul_bench 實作多執行緒的處理

• 使用更多的執行緒可以讓運算時間大幅縮短，以下為 6 個執行緒的結果
  

  

• thread function: 矩陣 inL 分成數等份，從 i_start 到 i_end 的區間資料做平行處理。

for (unsigned long i=i_start; i<i_end; i++) {
    for (unsigned long j=0; j<n; j++) {
        float v = 0;
        for (unsigned long k=0; k<n; k++) {
            v += inL[i*n+k] * inR[k*n+j];
        }
        out[i*n+j] = v;
    }
}

outer function

將向量總和操作放到迴圈外面，讓 out,inR,inL 皆能存取連續空間

for (unsigned long i=0; i<n; i++) {
    for (int k=0; k<n; k++) {
        float lik = inL[i*n+k];
        for (int j=0; j<n; j++) {
            out[i*n+j] += lik * inR[k*n + j];
        }
    }
}

• outer function: 不用轉置也能讓資料連續存取
  
  

使用 perf 分析 cache-misses, caceh-references, instructions, cycles 的指標，也觀察其中執行時間 sudo perf stat --repeat 5 -e cache-misses,cache-references,instructions,cycles ./matmul -n 2048 -a 1 ，在矩陣相乘的部份使用的 outer function 的方式，測量後的 outer function 比 transposed 更快

• matmul transposed

 Performance counter stats for './matmul -n 2048 -a 1' (5 runs):

       562,795,487      cache-misses                     #   47.72% of all cache refs           ( +-  3.09% )
     1,179,384,120      cache-references                                                        ( +-  2.47% )
    51,870,785,897      instructions                     #    1.40  insn per cycle              ( +-  0.00% )
    37,013,974,362      cycles                                                                  ( +-  0.08% )

            9.0445 +- 0.0109 seconds time elapsed  ( +-  0.12% )

• 在 matmul 中實作 outer function

 Performance counter stats for './matmul -n 2048 -a 8' (5 runs):

       464,438,845      cache-misses                     #   40.17% of all cache refs           ( +-  2.43% )
     1,156,173,078      cache-references                                                        ( +-  0.38% )
    60,405,629,722      instructions                     #    2.91  insn per cycle              ( +-  0.00% )
    20,736,333,706      cycles                                                                  ( +-  0.13% )

            5.0596 +- 0.0111 seconds time elapsed  ( +-  0.22% )

SIMD GCC vector type

Using Vector Instructions through Built-in Functions

GCC 提供了向量擴展（Vector Extensions），允許以相對較高階的方式使用向量指令，可以使用類似 C 語言的語法來操作向量，而編譯器會將這些操作轉換為目標架構的特定 SIMD 指令。

透過 __attribute__((vector_size(n))) 可以建立自訂大小的向量型態，例如 128-bit 和 256-bit 的向量，然後直接對這些向量進行運算，GCC 會自動將其轉換成 SIMD 指令。在以下範例中，展示 v4sf（128-bit 單精度浮點數向量）和 v8sf（256-bit 單精度浮點數向量）進行 SIMD 計算。

使用 v4sf (128-bit 向量)

這段程式碼使用 128-bit 的向量型態 v4sf 來加速矩陣運算。每次迴圈會同時處理 4 個浮點數，減少處理次數。

typedef float v4sf __attribute__((vector_size (16)));

for (unsigned long i=i_start; i<i_end; i++) {
    for (int k=0; k<n; k++) {
        float lik = inL[i*n+k];
        v4sf vlik = {lik,lik,lik,lik}; // process 4 instructions in 1 time

        for (int j=0; j<n; j+=4) {
            *(v4sf*)&out[i*n+j] += vlik * *(v4sf*)&inR[k*n + j];
        }
    }
}

使用 v8sf (256-bit 向量)

v8sf 允許處理更大的向量，進一步增加處理的浮點數數量。在支援 AVX 的硬體上，可以用 256-bit 向量一次處理 8 個浮點數。

typedef float v8sf __attribute__((vector_size (32)));

v8sf vlik = {lik,lik,lik,lik,
            lik,lik,lik,lik}; 

block matrix multiplication

與 loop tiling 類似，loop tiling 是針對迴圈的操作，而 block matrix 是將大矩陣劃分成小矩陣（block matrix），並對這些小矩陣分別進行運算，兩者皆是改善記憶體的局部性(locality)，減少了記憶體讀取的延遲，提高了計算效率，在實際應用中，可以根據具體的硬體配置選擇合適的分割大小 k，以達到最佳效能。

注意書寫規範：

• 程式碼註解不該出現中文，總是用美式英語書寫

for (i0=i_start; i0<i_end; i0+=block_size) {
    for (int j0=0; j0<n; j0+=block_size) {
        for (int k0=0; k0<n; k0+=block_size) {
            // cut into blocks size, then do block matrix multiplication.
            for (int bi=0; bi<block_size; bi++) {
                int i = i0+bi;
                for (int bk=0; bk<block_size; bk++) {
                    int k = k0+bk;
                    float lik = inL[i*n + k];
                    for (int bj=0; bj<block_size; bj++) {
                        int j = j0+bj;
                        out[i*n + j] += lik * inR[k*n + j]; // outer function
                    }
                }
            }
        }
    }
}

block matrix wiki

Other Intel SIMD Extensions

• AVX（Advanced Vector Extensions）

Intel AVX2

• FMA(Fused multiply-add)
  FMA（Fused Multiply-Add）是一種數學運算，可以在單一步驟中完成浮點數的乘法與加法運算。通過將兩個操作融合為一個，FMA 減少了指令數量，從而提高了計算速度。FMA 在硬體層面上同時執行
  $a+(b\times c)$ 的計算，而不是分別計算乘法和加法。這意味著 FMA 可以計算
  $b\times c$ 並將結果與
  $a$ 相加，然後將最終結果四捨五入到 N 個有效位。FMA 在一個操作中完成乘法和加法，減少了中間結果的舍入誤差。

Fused multiply–add wiki

測量不同實作方式的時間

TODO: Transformer 要解決什麼問題？

Transformer 跟過去的模型差在哪裡，解決什麼問題

Transformer 模型

Transformer（Attention is all you need）模型在 2017 年引入，為了解決機器翻譯中的長期依賴性問題，透過避免遞迴操作並使用平行化計算的方式減少訓練時間，提升訓練效率。其主要創新點如下：

1. 非序列處理 (Non Sequential)：
   • Transformer 將整句話同時處理，而不是逐字處理，允許並行處理，加快訓練速度。
2. 自注意力機制（Self-Attention）：
   • 自注意力機制能計算句中詞與詞之間的相似度分數，進而捕捉長距依賴性。
3. 位置嵌入（Positional Embeddings）：
   • 位置嵌入用來替代遞迴操作，編碼每個詞在句子中的位置訊息。

以下說明 Transformer 如何透過自注意力機制一步到位地取得全局關聯性。對於每個詞嵌入向量

$x_i$，計算 Query (

$Q$)、Key (

$K$) 和 Value (

$V$) 向量：

$$Q=x{W}_Q,K=x{W}_K,V=x{W}_V$$

其中

$W_Q$、

$W_K$ 和

$W_V$ 是學習的權重矩陣。

使用 Query 向量和 Key 向量計算 Attention scores：

$$\text{Attention}(Q,K,V)=\text{softmax}\left(\frac{Q{K}^T}{\sqrt{d_k}}\right)V$$

計算得到的 Attention score 矩陣是一個

$n\times n$ 的矩陣，其中

$n$ 是序列長度。這個矩陣的每個元素表示序列中兩個詞之間的關聯性。這些分數應用於 Value 向量，得到加權平均後的輸出向量：

$$\text{Output}=\text{softmax}\left(\frac{Q{K}^T}{\sqrt{d_k}}\right)V$$

在一次自注意力操作中，所有詞之間的關聯性都被計算出來。這意味著每個詞都能夠看到整個序列中的其他詞，從而捕捉到全局的關聯性。這與 CNN 逐層擴展視野的方法不同，Transformer 能在一步操作中完成全局關聯性的捕捉，且其中的計算可以被平行化。

RNN / LSTM 模型

RNN（遞迴神經網路）和 LSTM（長短期記憶網路）主要特性如下：

1. 序列處理(Sequential)：
   • RNN 和 LSTM 模型逐字處理輸入序列，每一步的輸出依賴於前一步的隱藏狀態 (基於 Markov property)，訊息傳遞受限於逐步更新，這種方式處裡句子只能逐字處理，無法並行計算。
2. 長期依賴性問題(Long-range Dependency)：
   • RNN 容易遇到梯度消失和爆炸 (gradient vanishing/exploding) 問題，難以捕捉長期依賴性。
   • LSTM 通過引入控制門 (Gate)，部分緩解這一問題，但仍受限於遞迴特性。

已知 RNN 模型主要是以序列處理，下一個輸入會是上一個輸出的結果，因此 RNN 的計算受限於時間，無法進行平行化。雖然透過局部平行化可以改善單個 RNN 節點的速度，但無法提升整體訓練速度。即使後續的 LSTM 和 GRU 透過門控機制解決長期依賴問題，它們依然是逐步序列處理，無法在根本上達到平行化的效益。

Recurrent neural network wiki

RNN 梯度消失與爆炸的原因

在 RNN 中，模型透過反向傳播更新參數。此過程中，梯度會反覆疊加並反向傳遞。當序列較長時，反向傳播的過程涉及將梯度多次相乘，尤其當 RNN 的隱藏狀態更新公式為：

$$h_t=\sigma (W_h\cdot h_{t-1}+W_x\cdot x_t)$$
這個更新公式中的權重矩陣

$W_h$ 反覆相乘，使得梯度在傳播過程中呈現指數增長或縮減。當權重絕對值大於 1 時，會引發梯度爆炸；相反，當權重小於 1 時，梯度會指數衰減而導致梯度消失。這些問題特別會影響較長序列的學習，使得 RNN 無法有效捕捉長期依賴。

Recurrent Neural Networks (RNNs) and Vanishing Gradients

LSTM/GRU 控制門 (Gate)

LSTM (Long Short-Term Memory) 和 GRU (Gated Recurrent Unit) 等架構透過控制門 (Gate) 的設計來解決梯度問題。具體來說，控制門透過引入額外的機制，讓模型能夠選擇性地保留或丟棄資訊，控制梯度的傳遞並保護長期記憶。以下是控制門的關鍵作用：

1. 輸入門 (Input Gate)：控制哪些新資訊需要被加入到記憶單元。輸入門會根據當前的輸入和先前的隱藏狀態生成一個範圍在 0 到 1 的控制信號，決定要寫入記憶的程度。

2. 遺忘門 (Forget Gate)：控制哪些資訊需要從記憶單元中移除。遺忘門會根據過去的狀態與當前輸入判斷應該遺忘的記憶比例。這樣可以避免無用資訊的累積，減少梯度爆炸的風險。

3. 輸出門 (Output Gate)：控制記憶單元中哪些資訊需要影響當前輸出的隱藏狀態。輸出門根據過去的記憶單元和當前的輸入決定應輸出的資訊，使得模型保留重要資訊。

CNN 模型

CNN（卷積神經網路）在 NLP 中也被廣泛應用，尤其適合短文本處理。其主要特性如下：

1. 卷積核（Kernels）：
   • CNN 依賴不同大小的卷積核來捕捉不同範圍的詞關係。例如，大小為 2 的卷積核可以學習詞對之間的關係，大小為 3 的卷積核可以捕捉三個詞之間的關係，以此類推。在這種方法模型可以有效捕捉局部特徵，非常適合處理短文本。
2. 多通道 CNN（Multichannel CNN）：
   • CNN 在文本分類中，通常會使用多種不同大小的卷積核來同時處理句子，這樣可以從不同的視角捕捉特徵。此外 CNN 的多通道(channel)能夠更全面地學習文本中的特徵，提升模型性能。
3. 限制：但隨著句子長度的增加，為了捕捉所有可能的依賴關係，所需的卷積核數量會呈指數級增長，這會導致計算資源需求巨大，實用性受限。

相較之下，CNN 的計算更容易平行化，因為卷積操作本質上可以平行處理。然而，CNN 也有其局限性。CNN 模型中兩個重要元件：卷積 (Convolution) 和池化 (Pooling)，都可以獨立操作，因而可以平行化。但隨著序列長度的增加，CNN 的運算成本也隨之倍增，主要有以下幾個原因：

1. 局部視野：卷積核一次只能看到固定大小的窗口（例如 3 個或 5 個詞），這意味著需要多層卷積層來擴大視野以捕捉更長距離的依賴性。

2. 層數增加：隨著序列長度增加，需要更多的卷積層來捕捉長距依賴性，這會導致模型變得更深，從而增加計算量和訓練時間。以下是層數增多的例子：

   1. 第一層：視野大小為 3，即每個卷積核看到 3 個詞。
   2. 第二層：如果第一層的輸出再經過第二層大小為 3 的卷積核，視野大小為 5，即
      $3+3-1$（相鄰重疊區域減 1）。
   3. 第三層：繼續重複這個過程，第三層的視野大小為 7，即
      $5+3-1$。

要捕捉長度為 100 的序列的長距離依賴性，需要多層卷積：

• 第 n 層的視野：
  $R_n=R_{n-1}+k-1$，其中
  $k$ 是卷積核的大小。
• 為了達到 100 個詞的視野，假設卷積核大小固定為 3，計算公式如下：
  
  $$R_n\approx 1+n(k-1)\text{（假設卷積層數 n 很大）}$$
  例如，為了達到視野 100，則需要大約
  $n\approx \frac{100-1}{3-1}=50$ 層。

具體來說，它解決了以下問題：

1. 長期依賴性問題：在傳統的RNN中，長期依賴性往往難以捕捉，因為它們在處理長序列時可能會出現 gradient vanish and explosion (梯度消失或爆炸)的問題。Transformer 通過引入 Attention mechanism，模型可以直接注意序列中的任何位置，從而更有效地捕捉長期資訊。
2. 並行化問題：RNN 的計算通常是順序進行的，因此難以實現有效的並行化。而 Transformer 中的自注意層允許所有位置之間的計算以並行方式進行，從而提高了訓練效率。
3. 計算效率問題：在傳統的 RNN 或 CNN 中，隨著序列長度的增加，計算成本可能急劇增加。而 Transformer 中的自注意層和位置編碼可以在不增加計算成本的情況下處理任意長度的序列。

參考：Why does the transformer do better than RNN and LSTM in long-range context dependencies?

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

找出其中效能瓶頸的地方，可能是矩陣運算，或是記憶體 (特別是 data model)

• 矩陣相乘：矩陣翻轉、矩陣相乘的平行化
• 模型參數大小：模型量化

矩陣相乘

llama.cpp 在導入 llamafile sgemm (commit 8cc91dc) 之前，使用 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

llamafile sgemm 之前，瓶頸出現在 ggml_vec_dot

Matrix vector multiplication is an operation where latency (not throughput) is the bottleneck, and the bloat of fancy libraries has a measurable impact.
LLaMA Now Goes Faster on CPUs

| model                          |       size |     params | backend    |    threads | test       |              t/s |
| ------------------------------ | ---------: | ---------: | ---------- | ---------: | ---------- | ---------------: |
| llama 1B Q4_0                  | 606.53 MiB |     1.10 B | CPU        |          6 | pp 512     |     88.24 ± 2.15 |
| llama 1B Q4_0                  | 606.53 MiB |     1.10 B | CPU        |          6 | tg 128     |     36.09 ± 0.18 |

build: dbceec87 (2684)
[ perf record: Woken up 525 times to write data ]
[ perf record: Captured and wrote 132.627 MB perf.data (2888859 samples) ]

# Samples: 1M of event 'cpu-cycles'
# Overhead  Command      Shared Object       Symbol                     
    72.52%  llama-bench  llama-bench         [.] ggml_vec_dot_q4_0_q8_0
    13.96%  llama-bench  llama-bench         [.] ggml_graph_compute_thread
     3.48%  llama-bench  llama-bench         [.] ggml_vec_dot_f16
     3.07%  llama-bench  llama-bench         [.] ggml_compute_forward_mul_mat
     2.26%  llama-bench  llama-bench         [.] ggml_vec_dot_q6_K_q8_K
     1.06%  llama-bench  llama-bench         [.] ggml_compute_forward_soft_max
     0.30%  llama-bench  llama-bench         [.] ggml_compute_forward_mul

llamafile sgemmg 之後

| model                          |       size |     params | backend    |    threads | test       |              t/s |
| ------------------------------ | ---------: | ---------: | ---------- | ---------: | ---------- | ---------------: |
| llama 1B Q4_0                  | 606.53 MiB |     1.10 B | CPU        |          6 | pp 512     |     85.39 ± 0.59 |
| llama 1B Q4_0                  | 606.53 MiB |     1.10 B | CPU        |          6 | tg 128     |     34.81 ± 0.18 |

build: dbceec87 (2684)
[ perf record: Woken up 547 times to write data ]
[ perf record: Captured and wrote 137.806 MB perf.data (3001985 samples) ]
# Samples: 1M of event 'cpu-cycles'
# Overhead  Command      Shared Object       Symbol
#
    41.44%  llama-bench  llama-bench         [.] (anonymous namespace)::tinyBLAS_Q0_AVX2<block_q4_0, block_q8_0, float>::gemm4x3(int, int, int, int)
    31.11%  llama-bench  llama-bench         [.] (anonymous namespace)::tinyBLAS_Q0_AVX2<block_q4_0, block_q8_0, float>::gemm4x1(int, int, int, int)
    14.65%  llama-bench  llama-bench         [.] ggml_graph_compute_thread
     4.13%  llama-bench  llama-bench         [.] (anonymous namespace)::tinyBLAS<8, float __vector(8), float __vector(8), unsigned short, float, float>::gemm3x4(int, int, int, int)
     2.97%  llama-bench  llama-bench         [.] ggml_vec_dot_q6_K_q8_K
     1.40%  llama-bench  llama-bench         [.] ggml_compute_forward_soft_max
     0.38%  llama-bench  llama-bench         [.] (anonymous namespace)::tinyBLAS<8, float __vector(8), float __vector(8), unsigned short, float, float>::gemm4x1(int, int, int, int)
     0.37%  llama-bench  llama-bench         [.] ggml_compute_forward_mul

模型量化

LLM通常使用全精度（float32）或半精度（float16）浮點數進行訓練。因此，量化過程尋找一種方法將 FP32 權重值的範圍（資料類型為 [min, max]）表示為較低精確度值，例如 FP16 甚至 INT4（整數 4 位元）資料類型。典型的情況是從 FP32 到 INT8。除了降低模型載入及儲存的空間，同時也因為精度下降，進而讓計算速度有所提升

然而精度的下降意味著模型預測的準確度會有所影響，可以藉由 Perplexity 指標確認量化模型後的品質

$$PPL(X)=exp\{-\frac{1}{t}\sum _i^{t}log{p}_{\theta }(x_i|x_{<i})\}$$
這表示模型在給定序列

$x_{<i}$ 的條件下預測下個詞

$x_i$ 的能力，而我們是希望在

$p_{\theta }(x_i|x_{<i})$ 的機率越大越好，條件機率的算法是

$\prod _i^t{p}_{\theta }(x_i|x_{<i})$，加上

$log$ 之後就可以方便計算其總和並取平均，得到一個數值。Perplexity 越低，表示模型預測越準確；值越高，表示模型預測越困難。
在 llama.cpp 中有探討不同大小的語言模型量化後的評估 k-quants #1684

參考：Quantization for Neural Networks

Quantization Mapping

將浮點數

$x\in [\alpha ,\beta ]$ 映射到整數

$x_q\in [{\alpha }_q,{\beta }_q]$

• de-quantization
  
  $$x=c(x_q+d)$$
• quantization
  
  $$x_q=round(\frac{1}{c}x-d)$$

計算極值的映射關係

$$\beta =c({\beta }_q+d),\alpha =c({\alpha }_q+d)$$

藉由極值

$\alpha ,\beta$ 計算

$c$,

$d$

$$c=\frac{\beta -\alpha }{{\beta }_q-{\alpha }_q},d=\frac{\alpha {\beta }_q-\beta {\alpha }_q}{\beta -\alpha }$$

最後 zero point z 表示原點，且 scale s 表示縮放數值

$x=s(x_q-z)$

$x_q=round(\frac{1}{s}x+z)$

透過極值

$\alpha ,\beta$ 取得 zero point 及 scale

$$s=\frac{\beta -\alpha }{{\beta }_q-{\alpha }_q},z=\frac{\beta {\alpha }_q-\alpha {\beta }_q}{\beta -\alpha }$$

實際使用中

$\alpha =min(x),\beta =max(x)$, 而量化極值可以是

${\alpha }_q=0,{\beta }_q=2^{n-bit}-1$, 通常這個數值範圍可以自訂

Quantized Matrix Multiplication

量化矩陣相乘過程，先知道矩陣相乘

$Y=XW+b$,

$X\in {\mathbb{R}}^{m\times p},W\in {\mathbb{R}}^{p\times n},b\in {\mathbb{R}}^n,Y\in {\mathbb{R}}^{m\times n}$
矩陣相乘數學表示式

$$Y_{i,j}=b_j+\sum _{k=1}^p{X}_{i,k}{W}_{k,j}$$

依照上面量化公式

$x=s(x_q-z)$代入

$$\begin{array}{rl}{Y}_{ij}& =b_j+\sum _{k=1}^p{X}_{i,k}{W}_{k,j}\\ & =s_b(b_{q,j}-z_b)+\sum _{k=1}^p{s}_X(X_{q,i,k}-z_X)s_W(W_{q,k,j}-z_W)\\ & =s_b(b_{q,j}-z_b)+s_X{s}_W\sum _{k=1}^p(X_{q,i,k}-z_X)(W_{q,k,j}-z_W)\\ & =s_b(b_{q,j}-z_b)+s_X{s}_W(\left(\sum _{k=1}^p{X}_{q,i,k}{W}_{q,k,j}\right)-\left(z_W\sum _{k=1}^p{X}_{q,i,k}\right)-\left(z_X\sum _{k=1}^p{W}_{q,k,j}\right)+p{z}_X{z}_W)\\ & =s_Y(Y_{q,i,j}-z_Y)\end{array}$$

從量化後的

$X_{q,i,k}$,

$W_{q,k,j}$,

$b_{q,j}$ 計算出量化結果

$Y_{q,i,j}$

$$\begin{array}{rl}{Y}_{q,i,j}=z_Y& +\frac{s_b}{s_Y}(b_{q,j}-z_b)\\ & +\frac{s_X{s}_W}{s_Y}(\left(\sum _{k=1}^p{X}_{q,i,k}{W}_{q,k,j}\right)-\left(z_W\sum _{k=1}^p{X}_{q,i,k}\right)-\left(z_X\sum _{k=1}^p{W}_{q,k,j}\right)+p{z}_X{z}_W)\end{array}$$

以下數值在模型推算時是常數，可以事先計算好

• $z_Y$
• $\frac{s_b}{s_Y}(b_{q,j}-z_b)$
• $z_X\sum _{k=1}^p{W}_{q,k,j}$
• $p{z}_X{z}_W$

而其中

$\sum _{k=1}^p{X}_{q,i,k}{W}_{q,k,j}$ 就是整數的矩陣相乘，可以使用硬體加速或是演算法處理

Quantized ReLU

• 經常用於深度模型的 ReLU
  

  $$\text{ReLU}(x,0,0,1)=\{\begin{array}{ll}0& \text{if }x<0\\ x& \text{if }x\ge 0\end{array}$$

• 用於表示量化 ReLU 且更通用的定義
  

  $$\text{ReLU}(x,z_x,z_y,k)=\{\begin{array}{ll}{z}_y& \text{if }x<z_x\\ z_y+k(x-z_x)& \text{if }x\ge z_x\end{array}$$

• 數學公式推導量化 ReLU

  • de-quantization
    
    $$\begin{array}{rl}y& =\text{ReLU}(x,0,0,1)\\ & =\{\begin{array}{ll}0& \text{if }x<0\\ x& \text{if }x\ge 0\end{array}\\ & =s_y(y_q-z_y)\\ & =\text{ReLU}(s_x(x_q-z_x),0,0,1)\\ & =\{\begin{array}{ll}0& \text{if }{s}_x(x_q-z_x)<0\\ s_x(x_q-z_x)& \text{if }{s}_x(x_q-z_x)\ge 0\end{array}\\ & =\{\begin{array}{ll}0& \text{if }{x}_q<z_x\\ s_x(x_q-z_x)& \text{if }{x}_q\ge z_x\end{array}\end{array}$$
    
    $$y=x\Rightarrow s_y(y_q-z_y)=s_x(x_q-z_x)$$
  • quantization
    
    $$\begin{array}{rl}{y}_q& =\{\begin{array}{ll}{z}_y& \text{if }{x}_q<z_x\\ z_y+\frac{s_x}{s_y}(x_q-z_x)& \text{if }{x}_q\ge z_x\end{array}\\ & =\text{ReLU}(x_q,z_x,z_y,\frac{s_x}{s_y})\end{array}$$

Quantization Implementation Unit Tests

驗證方法

1. 浮點數運算 (Floating Point Operations)：一般的神經網路浮點數運算。
2. 量化運算 (Quantized Operations)：量化參數運算，以 low bit-width 來加速推算。
3. 建立量化函數 (quantized function) 及反量化函數 (dequantized function)。
4. 建立浮點數 input tensor
   $(x,s_x,z_x)$ 和 output tensor
   $(y,s_y,z_y)$。
5. 執行浮點數運算，輸入為經過量化及反量化處理的 input tensor x，即
   $f_d(f_q(x),s_x,z_x)$，輸出為浮點數的 output tensor y，並計算
   $f_q(y)$ 及
   $f_d(f_q(y),s_y,z_y)$。
6. 執行量化運算，輸入為量化的 input tensor
   $x_q,s_x,z_x$，並提供 output tensor 的參數
   $s_y,z_y$，輸出為量化的 output tensor
   $y_q$，並計算
   $f_d(y_q,s_y,z_y)$。
7. 最後驗證
   $f_q(y)=y_q$ 的二進制相同，以及
   $f_d(f_q(y),s_y,z_y)=f_d(y_q,s_y,z_y)$ 的浮點數結果相同。

錯誤方式

1. 執行浮點數運算，取得浮點數的 output tensor y。
2. 執行量化運算，取得量化的 output tensor
   $y_q$，並計算
   $f_d(y_q,s_y,z_y)$。
3. 計算 y 與
   $f_d(y_q,s_y,z_y)$ 的誤差。

Neural Networks Integer Quantization Modes

• Dynamic Quantization
• Static Quantization
• Quantization Aware Training

Dynamic Quantization 動態量化

動態量化是在神經網路推理過程中使用，主要目的是使用整數操作來代替浮點數操作，以提高計算效率和降低資源消耗。在動態量化中，網路的權重在推理運行前被量化為整數，但網路並不知道 output tensor

$s_y$,

$z_y$ 及 activation tensors。而在執行階段中，只要一旦取得浮點數，即可以動態地計算出

$(\alpha ,\beta )$，進而算出中間層的 scale 及 zero point。

動態量化過程會直接算出

$Y_{i,j}$ 透過

$X_q$，前面有計算

$Y_{i,j}$ 的方式

優點：動態量化不需要任何預先校準資料，且是使用整數運算，動態量化仍能通過減少浮點運算來提升性能。
缺點：動態量化在計算 scale 和 zero points 時需要更多的運算，這可能在某些情況下增加計算負擔。

Static Quantization 靜態量化

靜態量化與動態量化的不同之處在於，靜態量化在推理前已經預先計算好所有 activation tensor 的 scales 和 zero points。因此，在推理期間不需要再計算這些參數。activation tensor 可以以整數格式儲存在記憶體中，而無需從浮點數格式轉換。

注意用語

• store 是「儲存」，不是「存儲」

計算所有 activation tensor 的 scales 和 zero points 的方法非常簡單。給定一個浮點數神經網路，我們只需使用一些代表性的無標籤資料運行該神經網路，收集所有 activation layer 的分布統計數據。然後，我們可以使用這些分布統計資料，按照前面的數學公式計算 scales 和 zero points。

在推理期間，由於所有計算都是使用整數運算，因此推理最快。唯一的缺點是我們必須準備代表性的無標籤資料。如果資料不具有代表性，計算出的 scales 和 zero points 可能無法反映到真實情況，從而影響推理準確性。

Quantization Aware Training 量化感知訓練

量化感知訓練在模型訓練期間進行模擬量化操作，使得模型在訓練時就能適應量化後的精度變化。

$$x=f_d(f_q(x),s_x,z_x)+{\mathrm{\Delta }}_x$$

${\mathrm{\Delta }}_x$ 是一個未知的小數值，如果

${\mathrm{\Delta }}_x=0$ 則表示量化模型的推理準確性將與浮點數模型完全一致，而這也是一個方式可以讓神經網路去學習之間的關係，並讓

${\mathrm{\Delta }}_x$ 的數值越小越好，並將之間的誤差降到最低。

這種量化技術不僅提高了模型的推理速度，還使得儲存需求更低，同時確保精度上有較小的損失。

探討其代價。