Day 15：換一個題目，這次我們要算 Hadamard

前面幾天我們一直在處理 GEMM。

那條主線的問題很清楚：我們想用 INT8 Tensor Core 做高精度矩陣乘法，所以一路補上 CRT、residue GEMM、DoubleDouble、Float128，以及 IEEE binary128 的轉換邊界。

今天開始題目換掉(換個主題水)。

這次我們不先問「怎麼把 GEMM 做得更準」，而是問另一個 GPU 上很常見、但性質完全不同的問題：

給一個長度 2^N 的 state vector x，
計算 y = H^{\otimes N} x

這是 quantum final state simulation 裡會碰到的 Hadamard transform。

如果只看符號，H^{\otimes N} 是一個 2^N x 2^N 的矩陣，元素只有 +1 和 -1。所以第一眼很容易想：

這不就是矩陣乘法嗎？
而且 H 只有 +1 / -1，看起來很適合 INT8 Tensor Core。

但這次我們不能這麼快跳到 Tensor Core。

因為 Hadamard matrix 不是一般 dense matrix。它有非常強的結構，可以用 Fast Walsh-Hadamard Transform，也就是 FWT，用 butterfly 直接做掉。

這一系列接下來要寫的，就是我們怎麼從這個數學結構一路拆到 GPU kernel，最後再回頭驗證：Tensor Core 在這題上到底幫不幫得上忙。

這幾天要怎麼切

我打算把這個 case 拆成 5 到 6 天(或是更多嘿嘿嘿)。

今天先講數學題目和計算模型。

也就是：

H^{\otimes N} 是什麼
為什麼它不用真的存成矩陣
FWT butterfly 怎麼把 O(2^N * 2^N) 變成 O(N * 2^N)
這個運算在 GPU 上真正吃的是什麼

Day 16 開始切第一段：single-block / intra-block FWT。

這天會進到 kernel 裡面看：

thread-local butterfly
float4 register blocking
warp shuffle
shared memory swizzle
N <= 14 / N = 15 的邊界

Day 17 切第二段：large-N multi-pass。

當 2^N 太大，單一 block 放不下時，就不能再幻想一個 kernel 吃完整個 vector。這時候會用 Kronecker decomposition，把整個 transform 拆成 row transform 和 column transform。

這天要看的問題會變成：

怎麼切 N
Pass 1 / Pass 2 分別在做什麼
為什麼一跨 block，global memory traffic 就變成主瓶頸

Day 18 切第三段：Tensor Core 嘗試。

這裡會回到前面幾天熟悉的東西，但不要急著套答案。我們會看：

Hadamard block 怎麼塞進 MMA
TF32 path 怎麼做
INT8 / Ozaki path 怎麼做
為什麼 layout conversion 和 transpose 會把理論算力吃掉

Day 19 講比較標的和 profile。

這天會把比較對象講清楚，包括 Dao-AILab 的 fast-hadamard-transform，還有我們後來在實驗室那邊看到的實作脈絡。這裡不是只看誰快，而是要看：

同一個 N、同一個 batch、同一個 normalization
output 是否對齊
latency 差在哪
NCU 看到的是 memory wall、occupancy、還是 instruction issue

Day 20 再做一次比較激進的版本。

如果前面的分析說明原本寫法輸在哪裡，那最後就應該用比較像 Razor / FlashAttention 那種精神，把路徑重新切乾淨：

不要 materialize 不必要的矩陣
不要做多餘 transpose
把資料盡量留在 register / shared memory
能 fuse 的 pass 就 fuse
用 profile 反推 kernel 邊界

這不是保證一定會贏，而是用工程方式把問題逼到硬體極限，最後誠實看結果。

數學題目：H 的 Kronecker 結構

單一 qubit 的 Hadamard 是：

H = 1/sqrt(2) * [ 1   1
                  1  -1 ]

如果有 N 個 qubit，就會得到：

H_N = H^{\otimes N}

輸入是一個 state vector：

x = [x_0, x_1, ..., x_{2^N-1}]

輸出是：

y = H_N x

如果直接把 H_N 展開成矩陣，再做矩陣乘法，計算量會是：

O(2^N * 2^N)

而且 H_N 本身就需要 2^N * 2^N 個元素。這對 N 稍微大一點就不合理。

例如 N=20：

vector length = 2^20 = 1,048,576
dense H size  = 2^20 * 2^20

這不是 kernel 寫好一點就能解決的問題。矩陣本身就不該被 materialize。

FWT：真正要算的是 butterfly

Hadamard 的結構讓我們可以改用 FWT。

FWT 每一層只做成對的加減：

a' = a + b
b' = a - b

每個 stage 的 stride 不同。

假設 stage 是 s，stride 就是：

stride = 2^s

每個 index i 的 partner 是：

j = i xor (1 << s)

所以 FWT 的整體形式是：

for s in 0 .. N-1:
    for every pair (i, i xor 2^s):
        a = x[i]
        b = x[j]
        x[i] = a + b
        x[j] = a - b

這樣計算量變成：

O(N * 2^N)

這就是整個 case_3 的核心。

我們不是要拿一個通用 GEMM 去乘一個巨大 Hadamard matrix，而是要把這個 butterfly network 映射到 GPU。

GPU 上真正的問題

FWT 的每次運算非常便宜。

一組 butterfly 只有：

1 個 add
1 個 subtract

所以它不像 GEMM 那樣可以用大量 multiply-add 把 memory cost 攤掉。

這會導致一個很重要的判斷：

FWT 很容易 memory-bound。

如果每個 stage 都從 global memory 讀一次、寫一次，那資料會被來回搬很多遍。

對 N 個 stage 來說，最壞的直覺版本會很接近：

每個 stage：
    read whole vector
    write whole vector

總 traffic 約 N * 2 * sizeof(vector)

所以這題的 kernel 設計，不是先想怎麼把 ALU 或 Tensor Core 塞滿，而是先問：

資料能不能留在 register？
warp 內能不能直接交換？
需要 shared memory 的時候能不能避免 bank conflict？
真的需要跨 block 時，global memory 要搬幾次？

這也是為什麼 case_3 裡面的 native path 會長成 register-heavy 的 butterfly kernel。

這跟前面的 GEMM 不一樣

前面做 Ozaki / CRT GEMM 時，我們願意付出 residue conversion、CRT reconstruction、甚至多個 modulus GEMM，是因為 GEMM 本身有足夠多的計算量。

對一個 M x K 乘 K x N 的 GEMM，每個 output element 會累積 K 次乘加。資料載入後可以被重複使用很多次。

但 FWT 不一樣。

每個 element 在每個 stage 只參與一次簡單加減。它沒有 GEMM 那種高 reuse 的 dense accumulation。

所以如果我們為了使用 Tensor Core，把 FWT 改寫成：

reshape
block GEMM
transpose
block GEMM
transpose back

那很可能還沒吃到 Tensor Core 的好處，就先被 layout conversion 和 memory traffic 打爆。

這不是理論問題，是後面幾天要用實作和 profile 去驗證的問題。

先看這個 project 的三條路

case_3 裡面大致有三條實作路線。

第一條是 Native FWT。

核心檔案是：

src/ImplicitHadamardNative.cu

這條路直接照 FWT 的 butterfly 結構寫 kernel。小 stride 用 thread-local 計算，warp 內用 __shfl_xor_sync，跨 warp 才用 shared memory。

第二條是 Tensor Core FP32 / TF32 path。

核心檔案是：

src/ImplicitHadamardTensorCoreFP32.cu

這條路嘗試把 Hadamard block 轉成 Tensor Core 可以吃的形狀。這很有吸引力，因為 Tensor Core 算力很高，但它會引入 layout 和 transpose 的成本。

第三條是 Ozaki / INT8 path。

核心檔案是：

src/ImplicitHadamardOzaki.cu

這條路把前面 GEMM 系列的想法帶進來：Hadamard matrix 的 +1/-1 很適合 INT8，input 則透過 scaling / residue 轉成 INT8，再用 MMA 做。

但這條路更需要小心，因為它除了 Tensor Core layout 問題，還多了 conversion 和 reconstruction。

Day 15 先停在這裡

今天先把題目定義清楚：

y = H^{\otimes N} x

它看起來像矩陣乘法，但真正該算的是 butterfly。

FWT 的成本不是 dense GEMM 的成本；它的硬體瓶頸也不是 Tensor Core TFLOPS，而是資料交換和 memory traffic。

接下來幾天我們就照這個切分往下走：

Day 16：single-block native butterfly kernel
Day 17：large-N multi-pass / Kronecker split
Day 18：Tensor Core / Ozaki 嘗試
Day 19：baseline、profile、為什麼原本路線不夠好
Day 20：用更激進的方式重切 kernel

這一輪的目標不是先選邊站，而是把問題拆到不能再拆，最後看硬體數據說話。