تحليل LU | LA 101

ما هو تحليل LU؟

تخيّل أنك تحتاج إلى حل المنظومة الخطية Ax = b ليس مرة واحدة، بل عشرات أو مئات المرات — في كل مرة مع متجه طرف أيمن b مختلف لكن مصفوفة المعاملات A ذاتها. هذا الوضع يظهر دائماً في الهندسة: محلّلات العناصر المنتهية تبني مصفوفة الصلابة K مرة واحدة ثم تحل Ku = f لعدد كبير من متجهات الحمل f؛ ومحاكيات الدوائر تبني مصفوفة التوصيلية G مرة ثم تحل Gv = i لمتجهات تيار كثيرة.

النهج الساذج — تشغيل حذف غاوس من الصفر في كل مرة — مضيّع للموارد. كل عمل اختزال A إلى الشكل السلّمي للصفوف متطابق بصرف النظر عن b. يلتقط تحليل LU ذلك العمل مرة واحدة للأبد بتحليل A إلى حاصل ضرب مصفوفتين مثلثيتين.

الفكرة المحورية هي أن عملية حذف غاوس على A يمكن تسجيلها كعاملين مثلثيين: مصفوفة مثلثية سفلية L بإدخالات قطرية تساوي 1، ومصفوفة مثلثية علوية U هي الشكل السلّمي للصفوف. بمجرد حساب A = LU، يُحل أي نظام Ax = b في خطوتين مثلثيتين رخيصتين.

العاملان L وU

المصفوفة المثلثية السفلية L هي مصفوفة تقع جميع مدخلاتها غير الصفرية على القطر الرئيسي أو تحته. الحالة الخاصة ذات الصلة بتحليل LU هي المصفوفة المثلثية السفلية الوحدوية — التي تكون إدخالاتها القطرية تساوي 1، وتحت القطر قيم اعتباطية، وفوقه أصفار. المصفوفة المثلثية العلوية U تحتوي جميع مدخلاتها غير الصفرية على القطر أو فوقه. يكتب تحليل LU:

تحليل LU

A = LU, \quad L = \begin{pmatrix} 1 & 0 & \cdots & 0 \\ l_{21} & 1 & \cdots & 0 \\ \vdots & \ddots & \ddots & \vdots \\ l_{n1} & l_{n2} & \cdots & 1 \end{pmatrix}, \quad U = \begin{pmatrix} u_{11} & u_{12} & \cdots & u_{1n} \\ 0 & u_{22} & \cdots & u_{2n} \\ \vdots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & u_{nn} \end{pmatrix}

تُعبَّر A على أنها حاصل ضرب المصفوفة المثلثية السفلية الوحدوية L والمصفوفة المثلثية العلوية U. تخزّن L المضاعفات المستخدمة في حذف غاوس؛ وU هي الشكل السلّمي للصفوف الناتج. يوجد التحليل متى أمكن اختزال A دون مبادلة صفوف — وبالمحاور الجزئية يوجد دائماً.

البنية أنيقة: L تُرمّز ما جرى فعله لاختزال A، وU تُسجّل النتيجة. معاً يحفظان كل حسابات حذف غاوس في شكل قابل للإعادة.

الصلة بحذف غاوس

يسير حذف غاوس باختيار محور في كل عمود وطرح مضاعفات صف المحور من جميع الصفوف تحته. المضاعف المستخدم لحذف المدخل في الصف i باستخدام المحور في الصف k هو:

m_ik = a_ik / a_kk

عادةً تُهمَل هذه المضاعفات بعد انتهاء الحذف. يحفظها تحليل LU: المضاعف m_ik يُخزَّن في الموضع (i, k) من L — بالضبط في الخانة التي أُصفِرت أثناء الحذف. كل خطوة حذف تحت القطر تملأ مدخلاً واحداً في L. المصفوفة المختزلة المتبقية هي U.

باختصار: U هي ما ينتجه حذف غاوس؛ L يتذكر كيف وصل إلى ذلك.

حل Ax = b عبر التعويض الأمامي والخلفي

بمجرد توفر التحليل A = LU، يصبح حل Ax = b عملية من خطوتين. تعويض A = LU في Ax = b يعطي LUx = b. بإدخال المتجه الوسيط y = Ux نحصل على منظومتين مثلثيتين:

التعويض الأمامي والخلفي

Ax = b \;\Rightarrow\; LUx = b \;\Rightarrow\; \begin{cases} Ly = b & \text{(forward substitution)} \\ Ux = y & \text{(back substitution)} \end{cases}

الخطوة 1 — التعويض الأمامي: حل Ly = b. لأن L مثلثية سفلية بإحاديات على القطر، فإن y₁ = b₁ مباشرةً، وكل yᵢ لاحقة تُحسب بتعويض القيم المحسوبة مسبقاً. الخطوة 2 — التعويض الخلفي: حل Ux = y. لأن U مثلثية علوية، فإن xₙ = yₙ/uₙₙ مباشرةً، وكل xᵢ سابقة تُحسب بالتعويض الخلفي.

التعويض الأمامي بالتفصيل

في التعويض الأمامي، تُحل المجاهيل من أعلى إلى أسفل. لأن L وحدوية مثلثية سفلية، فالمعادلة الأولى ببساطة y₁ = b₁. المعادلة الثانية هي l₂₁ y₁ + y₂ = b₂، فتعطي y₂ = b₂ − l₂₁ y₁. عموماً، في الصف i:

yᵢ = bᵢ − Σⱼ₌₁ⁱ⁻¹ lᵢⱼ yⱼ

يكلّف هذا O(n²) عملية إجمالاً. يسير التعويض الخلفي لـ Ux = y بشكل متماثل من أسفل إلى أعلى، بتكلفة O(n²) أيضاً.

لماذا LU؟ الكفاءة لأطراف أيمنة متعددة

الهدف الكامل من تحليل LU هو الفصل بين التكاليف. تحليل A = LU يتطلب تقريباً n³/3 ضرباً وجمعاً — استثمار أولي O(n³) لمرة واحدة. بمجرد تخزين L وU، يتطلب كل طرف أيمن b جديد تعويضاً أمامياً وخلفياً فقط، كلٌّ منهما بتكلفة O(n²).

التكلفة الحسابية

\underbrace{\frac{n^3}{3}\text{ flops}}_{{\text{factorize } A = LU}} + \underbrace{n^2\text{ flops}}_{{\text{each new } b}}

يكلّف التحليل ~n³/3 عملية فاصلة (flops) — مصروف لمرة واحدة. يكلّف كل حل ~n² عملية. لـ k طرف أيمن، التكلفة الإجمالية n³/3 + k·n² مقابل k·n³/3 لو أعدت تشغيل الحذف في كل مرة. التوفير يكبر بشكل هائل كلما زاد k أو n.

لنظام كبير بـ n = 10,000، يكلّف التحليل نحو 3.3 × 10¹¹ عملية — مكلف لكن لمرة واحدة. كل حل لاحق يكلّف 2 × 10⁸ عملية فقط، أرخص بنحو 1,600 مرة. هذا التباين هو سبب كون تحليل LU، لا حذف غاوس المتكرر، هو الخوارزمية المعتمدة في كل مكتبة حوسبة علمية جادة.

المحاور الجزئية للاستقرار العددي

يفترض تحليل LU الأساسي ألا يكون أي محور صفراً — ألا نضطر للقسمة على صفر أثناء الحذف. عملياً، حتى المحاور غير الصفرية يمكن أن تكون خطرة عددياً إن كانت صغيرة جداً: تصبح المضاعفات الناتجة كبيرة جداً مما يضخّم أخطاء التقريب في الحساب بالفاصلة العائمة.

المحاور الجزئية تتفادى ذلك بإعادة ترتيب الصفوف قبل كل خطوة حذف: عند العمود k، ابحث في الصفوف من k إلى n عن المدخل ذي القيمة المطلقة الأكبر، وانقل ذلك الصف إلى موضع المحور. يضمن هذا أن جميع المضاعفات تحقق |m_ik| ≤ 1 مما يحدّ نمو الأخطاء.

يُرمَّز إعادة ترتيب الصفوف بـمصفوفة التبديل P. بدلاً من تحليل A مباشرةً، تحسب المحاور الجزئية تحليل نسخة مبدَّلة من A:

LU مع المحاور الجزئية

PA = LU, \quad P \text{ permutation matrix encoding row swaps}

P هي مصفوفة تبديل تُرمّز جميع مبادلات الصفوف التي أُجريت أثناء الإمحاور. PA هي النسخة المبدَّلة (المعاد ترتيبها) من A، ويمكن دائماً تحقيق PA = LU. يصبح حل Ax = b: (١) تبديل b للحصول على Pb، (٢) تعويض أمامي Ly = Pb، (٣) تعويض خلفي Ux = y. في LAPACK وMAT LAB، هذا ما يجري داخل كل استدعاء لـ A\b.

المحاور الجزئية كافية في الغالب. المحاور الكاملة — التي تبدّل الأعمدة أيضاً — توفر ضمانات استقرار أقوى لكنها تتطلب تتبع مبادلات الأعمدة ونادراً ما تُستخدم عملياً.

مثال محلول 3×3

لنأخذ المصفوفة:

A = [[2, 1, 1], [4, 3, 3], [8, 7, 9]]

الخطوة 1 — حذف العمود 1. المحور هو a₁₁ = 2. المضاعفات هي m₂₁ = 4/2 = 2 وm₃₁ = 8/2 = 4. طرح 2 × الصف 1 من الصف 2 و4 × الصف 1 من الصف 3 يعطي مصفوفة جديدة بأصفار في العمود 1 تحت القطر، وتُخزَّن m₂₁ = 2 وm₃₁ = 4 في L.

الخطوة 2 — حذف العمود 2. المحور الجديد هو 1 (المدخل (2,2) بعد الخطوة 1). المضاعف هو m₃₂ = 1/1 = 1. طرح 1 × الصف 2 من الصف 3 يكمل U، ويُخزَّن m₃₂ = 1 في L.

مثال محلول 3×3

\begin{pmatrix}2&1&1\\4&3&3\\8&7&9\end{pmatrix} = \underbrace{\begin{pmatrix}1&0&0\\2&1&0\\4&1&1\end{pmatrix}}_{L} \underbrace{\begin{pmatrix}2&1&1\\0&1&1\\0&0&2\end{pmatrix}}_{U}

المصفوفة 3×3 هي A = [[2,1,1],[4,3,3],[8,7,9]] تتحلل إلى L (مثلثية سفلية بإحاديات على القطر، المضاعفات تحته) وU (الشكل السلّمي للصفوف). للتحقق: اضرب L × U وتأكد من الحصول على A. تخزّن L: m₂₁ = 2، m₃₁ = 4، m₃₂ = 1؛ وU هي المصفوفة بعد كلتا خطوتي الحذف.

لحل Ax = b مع، مثلاً، b = [1, 3, 9]ᵀ: يعطي التعويض الأمامي y₁ = 1، y₂ = 3 − 2·1 = 1، y₃ = 9 − 4·1 − 1·1 = 4. يعطي التعويض الخلفي x₃ = 4/2 = 2، x₂ = (1 − 1·2)/1 = −1، x₁ = (1 − 1·(−1) − 1·2)/2 = 0. إذن x = [0, −1, 2]ᵀ.

الأثر الهندسي

تحليل LU هو حصان العمل في الجبر الخطي العددي. مشغّل الشرطة المائلة العكسية في MATLAB (A\b)، ودالة numpy.linalg.solve في NumPy، وdgesv في LAPACK، كلها تحسب تحليل PA = LU في الكواليس. تبني محلّلات العناصر المنتهية مصفوفة الصلابة الكلية مرة واحدة ثم تحلها لعشرات حالات الحمل باستخدام العاملين L وU المخزَّنين. تبني محاكيات الدوائر مصفوفة القبولية العقدية المعدَّلة وتحلها في كل خطوة زمنية مع إعادة استخدام التحليل عندما لا تتغير الطبولوجيا. تحليل LU هو أيضاً أساس حساب معكوس المصفوفة (بحل AX = I عموداً عموداً) والمحدد (det A = حاصل ضرب الإدخالات القطرية لـ U، لأن det L = 1).

أبرز ما تعلمناه

يحلّل تحليل LU المصفوفة A = LU (أو PA = LU مع الإمحاور) بتسجيل مضاعفات حذف غاوس في L والشكل السلّمي للصفوف في U. يكلّف حل Ax = b بعد ذلك O(n²) فقط لكل طرف أيمن بعد استثمار أولي O(n³/3) لمرة واحدة. تضمن المحاور الجزئية — اختيار المحور ذي أكبر قيمة مطلقة في كل خطوة — الاستقرار العددي وهي معيارية في جميع التطبيقات الإنتاجية. البنية المثلثية لـ L وU هي ما يجعل الحل بخطوتين رخيصاً: يملأ التعويض الأمامي y من Ly = b من أعلى إلى أسفل، ويستعيد التعويض الخلفي x من Ux = y من أسفل إلى أعلى.