A GPU-resident Memory-Aware Algorithm for Accelerating Bidiagonalization of Banded Matrices
Ringoot, Alomairy, Edelman
The reduction of a banded matrix to a bidiagonal form is a crucial step in the Singular Value Decomposition (SVD), a cornerstone of scientific computing and AI. Despite being a highly parallel algorithm, it was previously believed to be unsuitable for GPU computation because it is memory bandwidth-bound. Recent developments in GPU hardware, including larger L1 memory per Streaming Multiprocessor/Compute Unit, have changed that. We present the first GPU algorithm for reducing a banded matrix to bidiagonal form as part of the NextLA.jl open-source software package. Our algorithm is based on previous CPU-based multicore parallel cache-efficient bulge chasing algorithms and adapted to optimize for GPU throughput. We leverage Julia Language's Array abstractions and KernelAbstractions to implement a single hardware- and data precision-agnostic function on NVIDIA, AMD, Intel, and Apple Metal GPUs for half, single, and double precision, and examine performance optimization across hardware architectures and data precision. We also develop a hardware-aware performance model and identify key hyperparameters, such as inner tilewidth and block concurrency, that govern optimal GPU execution for bandwidth-bound workloads. We demonstrate highly parallel bandwidth-bound algorithm on the GPU can outperform CPU-based implementations: the GPU algorithm outperforms multithreaded CPU High-Performance libraries PLASMA and SLATE as of matrix size 1024 x 1024 and by a factor over 100 for matrices of 32k x 32k. In addition, the performance of the algorithm increases linearly with matrix bandwidth size, making faster reduction of larger matrix bandwidths now also possible. With this work, we break memory bandwidth barriers, as well as matrix bandwidth barriers, resulting in orders-of-magnitude faster algorithms for the reduction of banded matrices to bidiagonal form on the GPU.
academic
خوارزمية موجودة في وحدة معالجة الرسومات وتدرك الذاكرة لتسريع ثنائي الأقطار للمصفوفات ذات النطاق
تقدم هذه الورقة أول خوارزمية موجودة في وحدة معالجة الرسومات وتدرك الذاكرة لتسريع اختزال المصفوفات ذات النطاق إلى مصفوفات ثنائية الأقطار، وهي خطوة حاسمة في تحليل القيم الذاتية (SVD). على الرغم من أن هذه الخوارزمية تتمتع بدرجة عالية من التوازي، إلا أنها كانت تُعتبر سابقاً غير مناسبة لحسابات وحدات معالجة الرسومات بسبب قيودها على نطاق ذاكرة التخزين المؤقت. مع تطور أجهزة وحدات معالجة الرسومات، خاصة مع وجود ذاكرة L1 أكبر لكل معالج متعدد التدفقات، تغير هذا الوضع. يعتمد المؤلفون على خوارزمية bulge-chasing السابقة الفعالة من حيث الذاكرة المؤقتة متعددة النوى على وحدة المعالجة المركزية، ويحسّنونها لإنتاجية وحدة معالجة الرسومات. تحقق الخوارزمية دالة واحدة مستقلة عن الأجهزة والدقة على وحدات معالجة الرسومات من NVIDIA و AMD و Intel و Apple Metal، مع دعم الحسابات بدقة نصف وواحد وثنائي. تُظهر التجارب أن خوارزمية وحدة معالجة الرسومات تتفوق على مكتبات وحدة المعالجة المركزية عالية الأداء متعددة الخيوط PLASMA و SLATE بدءاً من حجم مصفوفة 1024×1024، مع تحسن في الأداء يتجاوز 100 مرة على مصفوفات 32k×32k.
تحليل القيم الذاتية (SVD) هو أداة عددية أساسية في الحوسبة العلمية وتعلم الآلة وتحليل البيانات، مع تطبيقات واسعة في تحليل المكونات الرئيسية والفهرسة الدلالية الكامنة والتقريب منخفض الرتبة وإكمال المصفوفة. عادة ما يستخدم SVD على الأجهزة الكبيرة الحديثة عملية من ثلاث مراحل:
اختزال المصفوفة الكثيفة إلى مصفوفة ذات نطاق
اختزال المصفوفة ذات النطاق إلى مصفوفة ثنائية الأقطار (bulge-chasing)
بينما تمت دراسة تطبيقات وحدة معالجة الرسومات للمرحلة الأولى والثالثة على نطاق واسع، لم تُستكشف المرحلة الثانية بشكل كافٍ على وحدات معالجة الرسومات الحديثة. أشار Dongarra وآخرون في عام 2014 إلى أن "المسرعات تؤدي أداءً سيئاً في التعامل مع مهام الحوسبة الدقيقة المقيدة بالذاكرة (مثل bulge chasing)، مما يحد من الفوائد المحتملة لتطبيقات وحدة معالجة الرسومات في المرحلة الثانية".
أول خوارزمية bulge-chasing موجودة في وحدة معالجة الرسومات وتدرك الذاكرة: تقدم أول خوارزمية أصلية لوحدة معالجة الرسومات لاختزال المصفوفات ذات النطاق إلى مصفوفات ثنائية الأقطار، مما يستفيد بالكامل من خصائص الذاكرة الممتازة لأحدث أجيال وحدات معالجة الرسومات عالية الأداء، مع تحسن في الأداء بمعامل 10-100 مرة مقارنة بتطبيقات وحدة المعالجة المركزية المحسّنة.
استراتيجية bulge-chasing فعالة من حيث الذاكرة المؤقتة للمصفوفات ذات النطاق الكبير: من خلال استراتيجية جديدة فعالة من حيث الذاكرة المؤقتة لتقليل النطاق التدريجي، يمكن لخوارزمية وحدة معالجة الرسومات التعامل مع نطاقات أكبر من ذي قبل، مما يغير التوازن الأمثل للنطاق بين المرحلة الأولى والثانية في SVD الكبير.
تطبيق مفتوح المصدر مستقل عن الأجهزة والدقة: مكتبة مفتوحة المصدر بناءً على تجريدات لغة Julia، مع دالة واحدة تدعم الحسابات بدقة واحدة ونصف وثنائية على وحدات معالجة الرسومات من NVIDIA و AMD و Intel و Apple.
بالنظر إلى مصفوفة ذات نطاق A ∈ ℝⁿˣⁿ بنطاق BW، الهدف هو اختزال المصفوفة إلى مصفوفة ثنائية الأقطار B من خلال تحويلات متعامدة، بحيث A = UBVᵀ، حيث U و V مصفوفات متعامدة.
الخوارزمية 1: اختزال المصفوفة ذات النطاق إلى مصفوفة ثنائية الأقطار باستخدام متجهات Householder
الإدخال: النطاق BW، عرض البلاطة الداخلية TW، حجم المصفوفة n
1: لتقليل النطاق i = (BW-1)/TW → 1 افعل
2: النطاق المستهدف TBW = 1 + i·TW
3: لكل صف بالتوازي R = 1→n افعل
4: صف bulge k = R
5: لكل صف bulge: j = 0, j+=1 افعل
6: إذا كان 3(R-1) < j و k ≤ n ثم
7: احسب متجه HH للصف k لحذف عناصر TW وطبقه على الصفوف أدناه
8: احسب متجه HH لأعمدة bulge الأيسر الأكثر وطبقه على الأعمدة اليمنى
9: حدّث قيمة k
10: نهاية إذا
11: نهاية لكل
12: نهاية لكل
13: نهاية لكل
الخوارزمية 2: نواة داخلية لوحدة معالجة الرسومات تدرك الذاكرة، مع TPB خيط لكل كتلة
1: ذاكرة الخيط: Ai (TW+1)
2: ذاكرة الكتلة: X (TW+1)
3: تعاون جميع الخيوط في الكتلة: X ← A[k,..]
4: مزامنة الخيوط
5: تعاون جميع الخيوط في الكتلة: HH(X)
6: تعاون جميع الخيوط في الكتلة: A[k,..]← X
7: مزامنة الخيوط
8: لـ l: 0→(CBW + TW)/TPB - 1 افعل
9: الخيط i: احسب الصف r = k + l·CPB + i
10: Ai ← A[r, ...]
11: HH(X, Ai)
12: A[r, ...]← Ai
13: نهاية لـ
إدخال معامل Max blocks لتحديد عدد كتل الخيوط المتزامنة لكل وحدة تنفيذ في وحدة معالجة الرسومات. عندما يتجاوز عدد كتل bulge-chasing المطلوبة بواسطة الخوارزمية الحد، يتم استخدام فك الحلقة على مستوى البرنامج، حيث يتم تعيين كتلة واحدة مهام متعددة وتنفيذها بالتسلسل في إطلاق نواة واحد.
لضمان الصحة وتفعيل التوازي، تفرض الخوارزمية فصل ثلاث دورات بين عمليات المسح المتتالية للصفوف. أي أنه بعد إكمال ثلاثة صفوف bulge، يمكن للصف التالي أن يبدأ المسح دون تداخل الوصول إلى البيانات.
بخلاف محللات SVD، هناك عدد كبير من الأعمال على محللات القيم الذاتية للنطاق المتماثل ثنائي المرحلة على الأنظمة الهجينة CPU-GPU، لكن البحث في اختزال النطاق غير المتماثل إلى ثنائي الأقطار متخلف.
تم اقتراحها في الأصل كخوارزمية متوازية، لكن البحث في حالة ثنائي الأقطار على وحدة المعالجة المركزية نادر جداً. في المقابل، البحث في bulge-chasing للنطاق المتماثل إلى ثلاثي الأقطار والشكل Hessenberg العلوي أكثر انتشاراً.
تجاوز حاجز نطاق ذاكرة التخزين المؤقت: إثبات أن خوارزميات الجبر الخطي المقيدة بالذاكرة ليست فقط ممكنة على وحدات معالجة الرسومات الحديثة، بل يمكنها أن تتفوق بشكل كبير على أداء وحدة المعالجة المركزية
أهمية خصائص الأجهزة: نسبة تأخير L1/L2 إلى الحجم أكثر أهمية من الحجم نفسه، مما يؤكد أهمية الوصول السريع للذاكرة منخفضة التأخير
مبادئ تصميم الخوارزمية: الأداء الأمثل يأتي من تصميم الخوارزمية الأساسي والمعاملات التي تدرك الذاكرة المؤقتة والبنية الأساسية للمترجم المحمول
نموذج احتلال وحدة معالجة الرسومات: الخوارزمية لديها عدد كتل وحدة معالجة الرسومات منخفض بشكل كبير في المراحل الأولية مقارنة بوحدات التنفيذ المتاحة، مما يتطلب مصفوفات كبيرة بما يكفي للاستفادة الكاملة من موارد وحدة معالجة الرسومات
تجاوز السجلات: قد يؤدي عرض البلاطة الداخلية العالي إلى تجاوز السجلات، مما يؤثر على الأداء
الاعتماد على الأجهزة: الفروقات الكبيرة في الأداء بين معماريات وحدة معالجة الرسومات المختلفة تتطلب ضبط معاملات محدد للأجهزة
تستشهد الورقة بـ 86 مرجعاً ذا صلة، تغطي مجالات متعددة بما في ذلك حوسبة وحدات معالجة الرسومات وخوارزميات الجبر الخطي والنظام البيئي للغة Julia، مما يوفر أساساً نظرياً متيناً للبحث.