เหนือกว่าการคูณเมทริกซ์ของ NumPy ด้วยโค้ด C 150 บรรทัด
(salykova.github.io)- การคูณเมทริกซ์ของ NumPy พึ่งพา ไลบรารี BLAS ภายนอก แต่การติดตั้งใช้งานนี้ตั้งเป้าดันประสิทธิภาพแบบเธรดเดียวและหลายเธรดให้ถึงระดับ BLAS ด้วย C ล้วนและ FMA3·AVX2 เท่านั้น
- แกนหลักของประสิทธิภาพคือโครงสร้างที่แบ่ง $C$ ออกเป็นบล็อกเล็ก ๆ แล้วให้ ไมโครเคอร์เนล 16×6 ทำ rank-1 update ซ้ำภายในรีจิสเตอร์
YMMเพื่อลดการเข้าถึงหน่วยความจำ - สำหรับเมทริกซ์ขนาดใด ๆ การจัดการขอบเขตมักกลายเป็นคอขวดได้ง่าย จึงผสมผสาน mask store กับบัฟเฟอร์ padding ด้วย 0 เพื่อหลีกเลี่ยงประสิทธิภาพที่ลดลงจาก mask load
- การใช้แคชซ้ำทำได้ด้วยการ blocking
k_c,m_c,n_cและประสิทธิภาพสูงสุดจริงขึ้นอยู่กับการจูน จำนวนเธรด·ขนาดเคอร์เนล·ขนาดไทล์ อย่างมาก - ตัด AVX-512 ออกเพื่อรองรับ CPU ได้กว้างขึ้น ดังนั้นบน CPU ที่มี AVX-512 นั้น BLAS อาจเร็วกว่า และการเทียบกับ OpenBLAS ก็ทำภายใต้เงื่อนไขที่ปิด AVX-512
เป้าหมายการติดตั้งใช้งานและสิ่งที่ใช้เปรียบเทียบ
- โค้ดติดตั้งใช้งานเผยแพร่ไว้ที่ sgemm.c และปรับแต่ง การคูณเมทริกซ์ FP32 แบบหลายเธรด บนโปรเซสเซอร์รุ่นใหม่
- NumPy พึ่งพาไลบรารี BLAS ภายนอกสำหรับงานพีชคณิตเชิงเส้น เช่น การคูณเมทริกซ์
- ตัวอย่างได้แก่ Intel MKL, Accelerate, BLIS, GotoBLAS, OpenBLAS
- OpenBLAS, GotoBLAS, BLIS เขียนด้วย C/FORTRAN/Assembly และมีการติดตั้งใช้งานการคูณเมทริกซ์ที่ปรับแต่งด้วยมือสำหรับแต่ละ microarchitecture ของ CPU
- เป้าหมายคือการติดตั้งใช้งานการคูณเมทริกซ์ที่เขียนด้วย C ล้วน โดยไม่ใช้แอสเซมบลีระดับต่ำ แต่ยังตอบโจทย์เงื่อนไขต่อไปนี้
- ทำงานได้กับขนาดเมทริกซ์ใด ๆ
- รันบนโปรเซสเซอร์ x86-64 รุ่นใหม่
- แข่งขันกับไลบรารี BLAS ที่มีอยู่ได้
- โค้ดเรียบง่ายและขยายต่อได้ง่าย
- เอกสารอ้างอิงคือ Fast Multidimensional Matrix Multiplication on CPU from Scratch ของ Simon Boehm, Matrix Multiplication ของ Sergey Slotin, Can you multiply a matrix? ของ Geohot และงานวิจัยที่เกี่ยวข้องกับ GotoBLAS·BLIS
เงื่อนไขเบนช์มาร์กและการคำนวณ FLOPS
- สภาพแวดล้อมทดสอบคือ AMD Ryzen 7 9700X, 32GB DDR5 6000 MHz CL36, OpenBLAS 0.3.26, GCC 13.3, Ubuntu 24.04.1 LTS
- ใช้แฟล็กคอมไพล์
-O3 -march=native -mno-avx512f -fopenmp - เพื่อให้เปรียบเทียบอย่างเป็นธรรม เมื่อติดตั้ง OpenBLAS ต้องตั้งค่า
TARGETที่เหมาะสมและ ปิดใช้คำสั่ง AVX-512- โปรเซสเซอร์ Zen4/5 คอมไพล์ด้วย
make TARGET=ZEN - มิฉะนั้น OpenBLAS จะใช้คำสั่ง AVX-512 เป็นค่าเริ่มต้น
- โปรเซสเซอร์ Zen4/5 คอมไพล์ด้วย
- การคูณเมทริกซ์ FP32 ของ OpenBLAS รันผ่าน API
cblas_sgemm - เบนช์มาร์กทำกับเมทริกซ์จัตุรัส
- ประเมินตั้งแต่
m=n=k=200ถึงm=n=k=10000โดยเพิ่มทีละ200 - ทำการคูณเมทริกซ์ซ้ำ
n_iterครั้ง และใช้เวลารันค่ามัธยฐานในการวัดประสิทธิภาพ
- ประเมินตั้งแต่
- เมื่อคูณเมทริกซ์ $A$ ขนาด $M \times K$ กับเมทริกซ์ $B$ ขนาด $K \times N$ ปริมาณงานคำนวณรวมคือ $2MNK$ FLOP
- คำนวณประสิทธิภาพด้วย
FLOPS=(2*m*n*k)/exec_time
- คำนวณประสิทธิภาพด้วย
ขีดจำกัดเชิงทฤษฎีและพื้นฐาน SIMD
- CPU x86-64 รุ่นใหม่ประมวลผลข้อมูลหลายตัวพร้อมกันด้วยส่วนขยาย SIMD
- คำสั่งหลักคือ AVX2 และ FMA
- ทั้งคู่ใช้รีจิสเตอร์
YMMขนาด 256 บิต - รีจิสเตอร์
YMMแต่ละตัวเก็บ float 32 บิตได้ 8 ค่า
- ทั้งคู่ใช้รีจิสเตอร์
- คำสั่ง FMA
VFMADD231PSทำ packed single operation ในรูปYMM1 = YMM2 * YMM3 + YMM1 - บน Ryzen 9700X throughput ของ fused multiply-add คือ 0.5 cycles/instruction หรือ 2 คำสั่งต่อไซเคิล
- ตามทฤษฎี Ryzen 9700X สามารถทำงานได้ 32 FLOP ต่อไซเคิลบนคอร์เดียว
- สูตรคำนวณคือ
8 floats × 2(add+mul) × 2(1/TP) - หากสมมติคล็อกต่อเนื่อง 4.7GHz บน 8 คอร์ peak เชิงทฤษฎีแบบหลายเธรดจะประมาณ 1203 FLOPS
- สูตรคำนวณคือ
การติดตั้งใช้งานพื้นฐานและไมโครเคอร์เนล
- เมทริกซ์จัดเก็บแบบ column-major
A[row][col]เข้าถึงจากพอยน์เตอร์ C ด้วยptr[col*M + row]
- การติดตั้งใช้งานที่ง่ายที่สุดคือวนผ่านทุกแถวและคอลัมน์ของ $C$ แล้วคำนวณ dot product ของแถวใน $A$ กับคอลัมน์ใน $B$ สำหรับแต่ละองค์ประกอบ
- หัวใจของการติดตั้งใช้งานประสิทธิภาพสูงคือ ไมโครเคอร์เนล ที่แบ่ง $C$ ออกเป็นเมทริกซ์ย่อยขนาด $m_R \times n_R$ แล้วคำนวณเมทริกซ์ย่อยแต่ละส่วนอย่างมีประสิทธิภาพ
- เคอร์เนลเริ่มจากตั้งค่า $\bar{C}$ ในรีจิสเตอร์เป็น 0 แล้ววนซ้ำตามมิติ $K$
- นำเวกเตอร์คอลัมน์ของ $\bar{A}$ และเวกเตอร์แถวของ $\bar{B}$ เข้ารีจิสเตอร์
- คำนวณ outer product ของเวกเตอร์ทั้งสองแล้วบวกเข้า accumulator ของ $\bar{C}$
- แต่ละขั้นคือ rank-1 update
- วิธีนี้ลดจำนวนองค์ประกอบที่ต้องนำเข้ารีจิสเตอร์เหลือ $(m_R+n_R)K$ เมื่อเทียบกับปริมาณการเข้าถึงหน่วยความจำ $2K m_R n_R$ ของวิธี naive
- CPU ที่มี AVX มี รีจิสเตอร์ YMM 16 ตัว ดังนั้นขนาดเคอร์เนลต้องเป็นไปตามข้อจำกัดต่อไปนี้
- $(m_R/8) \cdot n_R + m_R/8 + 1 \le 16$
- $m_R$ ต้องเป็นพหุคูณของ 8
- ตามทฤษฎี ยิ่ง $m_R$ และ $n_R$ มีขนาดใหญ่และใกล้เคียงกัน การลดการเข้าถึงหน่วยความจำก็ยิ่งมาก แต่บน Ryzen 9700X จริง ๆ แล้ว เคอร์เนล 16×6 ให้ประสิทธิภาพดีที่สุด
- การติดตั้งใช้งานใช้ intrinsic จาก
immintrin.h__m256เป็นชนิดเวกเตอร์ 256 บิต และแทนเนื้อหาในรีจิสเตอร์YMM- โหลดเวกเตอร์คอลัมน์ของ
Aด้วย_mm256_loadu_ps - broadcast ค่าสเกลาร์ของ
Bเป็นเวกเตอร์ float 8 ค่า ด้วย_mm256_broadcast_ss - อัปเดต accumulator ด้วย
_mm256_fmadd_ps - เก็บผลลัพธ์ลงหน่วยความจำด้วย
_mm256_storeu_ps
- แอสเซมบลีที่สร้างขึ้นมีคำสั่ง SIMD FMA เช่น
vfmadd231psและvbroadcastss
Padding สำหรับเมทริกซ์ขนาดใด ๆ
- เคอร์เนล 16×6 พื้นฐานทำงานได้ทันทีเมื่อ $M$ และ $N$ เป็นพหุคูณของ 16 และ 6 ตามลำดับ
- ในบริเวณขอบ หากจำนวนคอลัมน์ $n$ น้อยกว่า 6 จะทำลูปจัดเก็บถึงแค่
j < n - เมื่อจำนวนแถว $m$ น้อยกว่า 16 เนื่องจาก
_mm256_storeu_psเก็บองค์ประกอบครั้งละ 8 ค่า จึงต้องใช้ mask store_mm256_maskstore_psจะเก็บลงหน่วยความจำเฉพาะองค์ประกอบที่บิต mask เปิดอยู่- mask สร้างตามจำนวนแถว $m$ ที่ทับซ้อน
- หากจัดการโหลดบริเวณขอบด้วย
_mm256_maskload_psด้วย ประสิทธิภาพของเคอร์เนลอาจตกลงมาก- คำสั่งเพิ่มเติมสำหรับคำนวณ mask ทำให้เกิด overhead
- เนื่องจาก $n$ ไม่ใช่ค่าคงที่ตอนคอมไพล์ คอมไพเลอร์จึง unroll ลูปอย่างมีประสิทธิภาพได้ยาก
- แทนที่จะทำเช่นนั้น หาก $m \neq m_R$ จะคัดลอก $\bar{A}$ ไปยังบัฟเฟอร์แล้ว padding ด้วย 0 และหาก $n \neq n_R$ ก็จะคัดลอก $\bar{B}$ ไปยังบัฟเฟอร์แล้วเติม 0 เช่นกัน
- การติดตั้งใช้งานที่เกี่ยวข้องอยู่ใน matmul_pad.h
Cache blocking และการใช้ข้อมูลซ้ำ
- ระหว่างรีจิสเตอร์กับ DRAM มีลำดับชั้นของ CPU cache และ CPU เดสก์ท็อปรุ่นใหม่โดยทั่วไปใช้แคช L1, L2, L3
- แคชเร็วกว่า DRAM แต่มีความจุจำกัด จึงเป็นไปไม่ได้ที่จะเก็บ $A$, $B$, $C$ ทั้งหมดไว้ในแคช
- วิธีแบ่งเมทริกซ์เป็นบล็อกเล็ก ๆ นำขึ้นแคช แล้วใช้ข้อมูลเดิมซ้ำใน rank-1 update หลายครั้ง เรียกว่า cache blocking หรือ tiling
- cache blocking แบบเธรดเดียวมีรูปแบบลูป 5 ชั้น คล้ายโครงสร้าง BLIS
- ลูปนอกสุดสร้างบล็อก $C_j$ และ $B_j$ ตามมิติ $N$
- ลูปถัดไปสร้างบล็อก $A_j$ และ $B_p$ ตามมิติ $K$
- $B_p$ ถูก pack เป็น $\tilde{B}_p$ และหากจำเป็นจะ padding ด้วย 0 เพื่อหวังใช้ซ้ำในแคช L3
- ลูปถัดไปสร้างบล็อก $C_i$ และ $A_j$ ตามมิติ $M$ และ $A_j$ ถูก pack เป็น $\tilde{A}_j$
- สองลูปสุดท้ายแบ่งบล็อกแคชเป็น panel ขนาด $m_R \times k_c$, $k_c \times n_R$ แล้วส่งต่อให้เคอร์เนล
- $\tilde{A}_j$ และ $\tilde{B}_p$ ที่ pack แล้วถูกจัดเก็บต่างกัน
- panel ภายใน $\tilde{A}_j$ ถูกเก็บแบบ column-major
- panel ภายใน $\tilde{B}_p$ ถูกเก็บแบบ row-major
- พารามิเตอร์ cache blocking ต้องปรับตามขนาดแคชของแต่ละรุ่น CPU
- $k_c \times n_c$ เป็นจุดเริ่มต้นสำหรับเติมแคช L3
- $m_c \times k_c$ เป็นจุดเริ่มต้นสำหรับเติมแคช L2
- $k_c \times n_R$ เป็นจุดเริ่มต้นสำหรับเติมแคช L1
- ในทางปฏิบัติ ค่าที่ใหญ่กว่าค่าทฤษฎีมักให้ประสิทธิภาพดีกว่า และเนื่องจาก CPU จัดการการวางข้อมูลในแคชโดยอัตโนมัติ จึงต้องออกแบบลูปและรูปแบบการเข้าถึงในระดับอัลกอริทึม
- การติดตั้งใช้งานอยู่ใน matmul_cache.h
การปรับแต่งละเอียดของเคอร์เนล
- แทนที่จะนิยาม accumulator เป็นอาร์เรย์ เช่น
__m256 C_buffer[6][2]ให้ประกาศตัวแปร accumulator แบบคลี่ออกอย่างชัดเจน - วิธีนี้ช่วยให้ GCC ปรับแต่งโค้ดได้ดีขึ้นและช่วยหลีกเลี่ยง register spilling
- การคำนวณ mask ก็เปลี่ยนให้ใช้คำสั่งเวกเตอร์ด้วย
- มีอาร์เรย์ static
mask[32]แล้วใช้_mm256_cvtepi8_epi32และ_mm_loadu_si64
- มีอาร์เรย์ static
- การติดตั้งใช้งานดังกล่าวอยู่ใน matmul_micro.h
กลยุทธ์มัลติเธรดดิ้ง
- สิ่งที่ทำ parallelization คือทั้งการคำนวณเชิงเลขและการ packing
- ลูปที่ 5, 4, 3 นอกไมโครเคอร์เนลวนซ้ำตามขนาดบล็อกแคช
- เพื่อให้ทุกเธรดทำงานตลอด จำนวนรอบวนซ้ำต้องไม่น้อยกว่าจำนวนเธรด
- มิติของเมทริกซ์อินพุตควรมีขนาดประมาณ
จำนวนเธรด × ขนาดบล็อกแคชขึ้นไป
- ขนาดบล็อกแคชที่ให้ประสิทธิภาพดีบน Ryzen 9700X แบบเธรดเดียวคือ $n_c=1535$, $m_c=1024$
- หากต้องการใช้ครบ 8 คอร์ ต้องมีมิติขนาดอย่างน้อย $\max(m_c,n_c) \times 8 = 1535 \times 8 = 12280$
- ในทางกลับกัน สองลูปสุดท้ายวนซ้ำบนบล็อกขนาดเล็ก $m_R$, $n_R$ จึงเหมาะกับการ parallelize
- โดยทั่วไป $m_R$, $n_R$ น้อยกว่า 20
- หากเลือก $m_c$, $n_c$ เป็นพหุคูณของจำนวนคอร์ ก็แบ่งงานได้อย่างสม่ำเสมอ
- บน Ryzen 9700X วิธีที่ให้ประสิทธิภาพดีที่สุดคือ parallelize ลูปภายในสองชั้นร่วมกันด้วย
#pragma omp parallel for collapse(2) num_threads(NTHREADS) - สำหรับโปรเซสเซอร์ที่มีคอร์จำนวนมาก โดยเฉพาะสภาพแวดล้อมที่เกิน 16 คอร์ อาจพิจารณา nested parallelism และการ parallelize ลูป 2–3 ชั้น
- การ packing $\tilde{A}$ และ $\tilde{B}$ ก็ parallelize ด้วย OpenMP เช่นกัน
pack_blockAparallelize โดยวนmcทีละMRpack_blockBparallelize โดยวนncทีละNR
- พารามิเตอร์ที่ให้ประสิทธิภาพดีบน Ryzen 9700X ในการติดตั้งใช้งานแบบหลายเธรดมีดังนี้
- $m_c = m_R \times \text{number of threads} \times 5$
- $n_c = n_R \times \text{number of threads} \times 50$
- การติดตั้งใช้งานแบบหลายเธรดขั้นสุดท้ายอยู่ใน matmul_parallel.h
1 ความคิดเห็น
ความคิดเห็นจาก Hacker News
ถ้าประเด็นของบทความนี้คือโดยทั่วไปแล้วยังมี ช่องว่างด้านประสิทธิภาพเหลืออยู่ ก็ถือว่าประเมินขนาดของการปรับปรุงต่ำไปด้วยซ้ำ ทั้งที่ความพยายามที่ใส่ลงไปในไลบรารีคูณเมทริกซ์นั้นมากกว่าซอฟต์แวร์ส่วนใหญ่หลายเท่า
ถ้าไม่ใช่โค้ดที่ถูกปรับแต่งอย่างหนักอยู่แล้ว การปรับปรุงจากโค้ดเดิมได้ 10~1000 เท่าขึ้นไป โดยไม่ต้องลงแรงมากนักเป็นเรื่องที่พบได้บ่อย หากเรียงตามความสำคัญคร่าว ๆ สิ่งที่สำคัญที่สุดคือการเลือกอัลกอริทึมให้เหมาะสมหรือไม่ และสามารถตัดงานนั้นทิ้งไปเลยได้หรือไม่ รวมถึงการลดงานหนัก ๆ อย่างการไป-กลับระหว่างเคอร์เนลหรือ
mallocก็มีผลมากเช่นกันการทำเวกเตอร์ไรซ์ด้วย intrinsic แบบชัดเจนก็เป็นเรื่องดี แต่บ่อยครั้งแค่จัดโครงสร้างข้อมูลใหม่จากอาร์เรย์ของสตรักต์ ไปเป็นอาร์เรย์ของอาร์เรย์/สตรักต์ ก็ได้ภาษาเครื่องแบบเดียวกันแล้ว ประสิทธิภาพของแคชก็สำคัญ และในโค้ดแบบขนานจะยิ่งซับซ้อนขึ้นเมื่อข้อมูลของแต่ละเธรดไม่ได้ถูกแยกจากกัน เช่น false sharing สุดท้ายยังมีการปรับแต่งเฉพาะฮาร์ดแวร์ได้ด้วย เช่น intrinsic หรือแอสเซมบลีที่เขียนเอง
พอเปลี่ยนคิวรีให้ join เกิดขึ้นบนเซิร์ฟเวอร์ระยะไกล และส่งผ่านเครือข่ายมาแค่ 5~10 แถว ก็เร็วขึ้นทันที แม้จะมีโอเวอร์เฮดคงที่และ latency อยู่เสมอ แต่ถ้าส่งข้อมูลผ่านการเชื่อมต่อเครือข่ายมากเกินกว่าที่จำเป็นมาก ๆ สุดท้ายประสิทธิภาพก็พังอยู่ดี บทความ “It's the latency, stupid” ที่พูดถึงผลของ latency ก็น่าอ่าน: http://www.stuartcheshire.org/rants/latency.html
โดยรวมแล้วเห็นด้วยกับประเด็นที่กล่าวมาข้างต้นและลำดับความสำคัญคร่าว ๆ
การสัมภาษณ์จำนวนมาก แทนที่จะดูวิธีอนุมานว่า implementation ช้าเพราะอะไร ทำ benchmark และแก้ไขอย่างไร กลับกลายเป็นควิซท่องจำอัลกอริทึม obscure แบบ “เพราะ Google ทำแบบนั้น”
แพตเทิร์นการเขียนโค้ดที่พบบ่อยไม่ได้เฉพาะเจาะจงกับฮาร์ดแวร์มากพอ จึง ทิ้งประสิทธิภาพไว้มาก บทความนี้เป็นตัวอย่างที่น่าสนใจ และยังมีการสาธิตคลาสสิกอีกอย่างคือ “There's plenty of room at the top”
https://www.science.org/doi/10.1126/science.aam9744
ถ้าจะทำความเข้าใจเรื่องนี้ บทความวิจัยในคลัง BLIS ถือว่าใกล้เคียงกับเอกสารอ้างอิงหลักที่สุด ผมไม่เข้าใจว่าทำไมถึงคิดว่า BLAS ที่ปรับแต่งแล้วจะทำประสิทธิภาพได้ไม่ดี และถ้าเมทริกซ์ใหญ่พอ ก็ควรคาดหวังได้มากกว่า 90% ของพีก CPU
ครั้งล่าสุดที่ดู OpenBLAS แบบ serial โดยรวมก็ใกล้เคียงกับ MKL และ BLAS ในฐานะบล็อกพื้นฐานของพีชคณิตเชิงเส้นนั้นไม่ได้ implement
matmulแต่ implement GEMM โดยปกติผมก็ไม่เข้าใจว่าทำไมใช้ numpy แทนเฟรมเวิร์ก benchmark และบน Zen ก็ควรเทียบกับ BLAS ของ AMD หรือ implementation ที่อิง BLIS มากกว่า แต่ก่อน BLIS มีเรื่องราวด้านการ parallelization ดีกว่า OpenBLAS และ AMD BLIS ยังมีการสลับ implementation สำหรับมิติที่ “เล็ก” ด้วย แต่ไม่รู้ว่าปัจจุบัน OpenBLAS มีหรือยังการทำเวกเตอร์ไรซ์ micro-kernel ไม่จำเป็นต้องใช้ SIMD intrinsic เสมอไป และคอมไพเลอร์ C ที่ดีจะทำเวกเตอร์ไรซ์ได้เต็มที่และ unroll ลูปให้ด้วย micro-kernel แบบ C ล้วนของ BLIS ที่ใช้ขนาดบล็อกเหมาะสม ทำประสิทธิภาพได้มากกว่า 80% เมื่อเทียบกับ implementation ที่ปรับแต่งด้วยมือบน Haswell ความต่างน่าจะมาจาก prefetch แต่ผมก็ไม่ได้เข้าใจอย่างแม่นยำ
แม้แต่คอมไพเลอร์รุ่นใหม่ก็ยังไม่สามารถทำ auto-vectorization และ unroll ลูปได้ถูกต้องด้วยอัตราสำเร็จ 100%
บทความและการนำไปทำดูดี แต่ก็สงสัยว่า “เคล็ดลับ” คืออะไร OpenBLAS ถูก ปรับแต่งด้วยแอสเซมบลี+C มานานหลายสิบปีเพื่อโจทย์แบบนี้โดยเฉพาะ แล้วมันเอาชนะได้อย่างไร
บทความลงรายละเอียดเรื่องแคช ฯลฯ เลยสงสัยว่า BLAS ไม่ได้ใช้ประโยชน์จากสิ่งเหล่านี้หรือเปล่า หรือว่าแค่ปรับให้เข้ากับโปรเซสเซอร์บางรุ่นได้ดีกว่า
ความต่างของประสิทธิภาพเห็นชัดที่ปริมาณงานสูงสุดมากกว่าปริมาณงานเฉลี่ยมาก แต่แทบไม่มีแอปพลิเคชันไหนที่ปริมาณงานสูงสุดสำคัญ โค้ดเบนช์มาร์กที่แสดงดูเหมือนว่าฝั่ง numpy จะผ่านตัวจัดสรรหน่วยความจำของ Python ส่วนการทำงานใน C ไม่ได้ผ่านตัวจัดสรร ดังนั้นนี่เป็นจุดแรกที่ควรตรวจสอบข้อผิดพลาดหรือความไม่สอดคล้องของไมโครเบนช์มาร์ก รูทีนของ numpy จำนวนมากรองรับการคำนวณแบบ in-place ดังนั้นน่าจะต้องดูเบนช์มาร์กเวอร์ชัน in-place ของทั้งสองฝั่งอย่างชัดเจน
numpy ยังมีการตรวจขอบเขตและการจัดการข้อผิดพลาดที่ทำงานโดยไม่ขึ้นกับ implementation ชั้นล่างด้วย ซึ่งเป็นเหตุผลว่าทำไมเมทริกซ์เล็ก ๆ ถึงช้ากว่าแม้กระทั่งลิสต์ Python ล้วน ๆ มาก เมื่อบวกโอเวอร์เฮดล้วน ๆ ระดับไม่กี่พันไซเคิลเข้าไป ก็ทำให้เร็วได้ยาก
implementation นี้เป็นแนวทางที่ค่อนข้างมีหลักการเพื่อทำให้แคชที่เกี่ยวข้องถูกใช้อย่างเต็มที่ และในแง่หนึ่งก็เป็นเรื่องชัดเจนอยู่แล้ว แต่การปรับปรุงทางวิศวกรรมที่ชัดเจนแบบนี้ควรถูกเน้นในการถกเถียงเช่นนี้ OpenBLAS ก็ใช้กำลังคนไปมาก แต่ไม่น่าจะคิดครอบคลุมทุกอย่างได้ทั้งหมด หากจะอธิบายให้ถูกต้องจริง ๆ ต้องวิเคราะห์โค้ดทั้งสองฝั่งอย่างลึกซึ้ง
สำหรับ implementation ใน C++ และ C ให้ดูแนวทางเมตาโปรแกรมมิง [2], [3] สิ่งที่น่าประหลาดใจจริง ๆ คือภาษาสมัยใหม่จำนวนมากอย่าง Matlab, Julia, Mojo ยังพึ่งพา OpenBLAS อยู่ ซึ่งแน่นอนว่าแต่ละภาษาคงมีเหตุผลของตัวเอง
[1] Numeric age for D: Mir GLAS is faster than OpenBLAS and Eigen (2016):
http://blog.mir.dlang.io/glas/benchmark/openblas/2016/09/23/...
[2] Vastly outperforming LAPACK with C++ metaprogramming (2018):
https://wordsandbuttons.online/vastly_outperforming_lapack_w...
[3] Outperforming LAPACK with C metaprogramming (2018):
https://wordsandbuttons.online/outperforming_lapack_with_c_m...
-march=nativeอาจได้เปรียบ เพราะคอมไพล์ให้ตรงกับรุ่น CPU ที่แน่นอน ส่วน numpy มีโอกาสสูงว่าจะถูกคอมไพล์สำหรับเป้าหมาย x86-64 ที่เก่ากว่าและทั่วไปกว่าบน CPU Ryzen นั้น
-march=nativeน่าจะใช้ v4 ส่วน numpy น่าจะตั้งเป้าไว้ที่ v1 หรือ v2https://en.wikipedia.org/wiki/X86-64#Microarchitecture_level...
บทความดี และยอดเยี่ยมที่ทำให้ทำซ้ำเบนช์มาร์กได้ง่าย บน 16-core Xeon W-2245 3.90GHz ของผม
matmul.cทำการคูณเมทริกซ์ 8192x8192 ได้ใน 1.41 วินาทีด้วยgcc -O3และ 1.47 วินาทีด้วยclang -O2ส่วน NumPy ใช้เวลา 1.07 วินาทีผมคิดว่าถ้าเป็น เคอร์เนล AVX-512 จะเร็วกว่านี้มาก อีกเหตุผลหนึ่งที่ประสิทธิภาพยังไม่ดีอาจเป็น OpenMP จากประสบการณ์ หากจัดการเธรดพูลเองอย่างชัดเจนด้วย
pthreadsจะลดโอเวอร์เฮดได้ และควรใช้sysconf(_SC_NPROCESSORS_ONLN)แทนการฮาร์ดโค้ดจำนวน CPUไม่มีเหตุผลที่จะทำให้ภาระต่างกันโดยฝั่งหนึ่งเป็น Python อีกฝั่งเป็น C ควรเขียนทั้งสองฝั่งด้วย C แล้วให้ฝั่งหนึ่งเรียก ไลบรารี BLAS อีกฝั่งเรียก implementation นี้ แบบนั้นถึงจะเปรียบเทียบกันแบบแอปเปิลกับแอปเปิลได้
โอเวอร์เฮดไม่ได้ใหญ่มาก แต่ตามที่มีคนพูดในที่อื่นของเธรดนี้ สิ่งสำคัญคือต้องเรียกใช้อย่างถูกต้อง การเอาโค้ด numpy แบบตรงไปตรงมามาชนกับโค้ด C ที่ปรับแต่งแล้ว แน่นอนว่าไม่ใช่การเปรียบเทียบที่ยุติธรรม
แม้จะไม่ใช่ hot path แต่ความไม่มีประสิทธิภาพในการสร้าง mask หรือก็คือการใช้
bit_maskมันขัดตา วิธีที่มีประสิทธิภาพกว่าน่าจะเป็นการสร้างอาร์เรย์ค่าคงที่แบบ global ในรูป{-1,-1,...,0,0,...}แล้วโหลดจาก offset ของ element ที่16-m,8-mหรือใช้วิธีเปรียบเทียบเวกเตอร์ค่าคงที่{0,1,2,3,4,...}กับmและm-8ที่ถูก broadcastแต่เพราะมันเกี่ยวข้องกับแค่คอลัมน์เดียวของเมทริกซ์ และลูป
maskload/maskstoreที่ตามมากินเวลานานกว่ามาก นี่จึงเป็นการจิกกัดเล็กน้อยจริง ๆ โดยเฉพาะฝั่ง store ยังช้าอยู่แม้แต่บน Zen 4[1] และคำสั่ง AVX-512 เร็วกว่า 6 เท่า ทั้งที่ต่างกันแค่รับ mask จาก mask register เท่านั้น ส่วน clang ยังไงก็ auto-vectorize การ shift อยู่แล้ว จึงน่าจะช้ากว่าข้อเสนอของผมราว 2–3 เท่า[1]: https://uops.info/table.html?search=vmaskmovps&cb_lat=on&cb_...
ขอบคุณมากสำหรับ feedback ที่ให้มุมมองใหม่ ๆ วิธี “สร้างอาร์เรย์ global ค่าคงที่แล้วโหลด” เท่าที่จำได้ ตอนทดสอบมันเหมือนจะช้ากว่าการ shift bit mask เล็กน้อย แต่จะทดสอบใหม่อีกครั้งให้แน่ใจ ส่วนวิธี “เปรียบเทียบเวกเตอร์ค่าคงที่
{0, 1, 2, 3, 4, ...}กับmและm-8ที่ถูก broadcast” เป็นไอเดียที่ดี จะลองทำดูint8_tแล้วตอนโหลดค่อย sign-extend byte เป็นint32_tได้ ชุด_mm_loadu_si64/_mm256_cvtepi8_epi32น่าจะ compile เป็นคำสั่งvpmovsxbdเดียวที่ใช้ memory operandแบบนี้ถ้า align ให้ถูกต้องด้วย
alignas(32)อาร์เรย์ค่าคงที่ทั้งชุดจะเข้าไปอยู่ใน cache line เดียวพอดี ใน use case ของต้นฉบับต้องใช้ mask สองตัว ดังนั้นคำสั่งvpmovsxbdตัวที่สองจะเป็น L1D cache hit แน่นอน จึงค่อนข้างเหมาะแล้ว tinyBLAS ของ jart ล่ะ
https://hacks.mozilla.org/2024/04/llamafiles-progress-four-m...
และ https://justine.lol/matmul/
นอกจาก benchmark แล้ว เหตุผลที่ต้องทำให้การคูณเมทริกซ์เองเป็นแบบ multi-threaded คืออะไร ในทางปฏิบัติ เอา multi-threading ไปใช้ในฝั่งอัลกอริทึมที่เรียกใช้การคูณจะได้ประโยชน์มากกว่าไม่ใช่หรือ
แต่โค้ด HPC ส่วนใหญ่ GEMM ไม่ได้เป็นคอขวด
ยังแค่อ่านผ่าน ๆ แต่บทความนี้มีรายละเอียดและคำอธิบายเยอะ ดูเหมือนเป็นบทความที่อธิบายได้ค่อนข้างดีว่า การคูณเมทริกซ์ที่เร็วถูก implement โดยคำนึงถึง ข้อพิจารณาด้านสถาปัตยกรรม อย่างไร เลยเก็บไว้ในรายการอ่านแล้ว