• برنامج تعليمي حول نظرية الاحتمالات (مقالة تمهيدية ، اختيارية)
  
• مقدمة في أنظمة المعادلات الخطية
  
• طرق العناصر المحدودة
  
• التقليل من الأشكال التربيعية وأمثلة لمشاكل OLS
  
• من المربعات الصغرى إلى الشبكات العصبية
  

لذا ، مقال آخر من دورة "الرياضيات على الأصابع ". اليوم نواصل المناقشة حول طرق المربعات الصغرى ، ولكن هذه المرة من وجهة نظر المبرمج. هذه مقالة أخرى في السلسلة ، لكنها متميزة ، لأنها لا تتطلب عمومًا أي معرفة بالرياضيات. تم تصور المقالة كمقدمة للنظرية ، لذلك من المهارات الأساسية تتطلب القدرة على تشغيل الكمبيوتر وكتابة خمسة أسطر من التعليمات البرمجية. بالطبع ، لن أسهب في الحديث عن هذه المقالة ، وفي المستقبل القريب سأقوم بنشر تكملة. إذا تمكنت من إيجاد الوقت الكافي ، فسوف أكتب كتابًا من هذه المادة. الجمهور المستهدف مبرمجون ، لذلك فإن هبر هي المكان المناسب للاقتحام. بشكل عام ، لا أحب كتابة الصيغ ، وأحب حقًا أن أتعلم من الأمثلة ، يبدو لي أن هذا مهم جدًا - ليس فقط النظر إلى المتاعب على لوحة المدرسة ، ولكن جرب كل شيء على السن.

لذلك دعونا نبدأ. دعونا نتخيل أن لدي سطحًا مثلثيًا مع مسح وجهي (في الصورة على اليسار). ما الذي أحتاجه لتحسين الميزات ، وتحويل هذا السطح إلى قناع غريب؟



في هذه الحالة بالذات ، أحل معادلة تفاضلية إهليلجية تسمى Simeon Demi Poisson . زملائي المبرمجين ، دعنا نلعب لعبة: خمن عدد الأسطر الموجودة في كود C ++ التي تقررها؟ لا يمكن استدعاء مكتبات الطرف الثالث ، لدينا مترجم فقط تحت تصرفنا. الجواب تحت القطع.

في الواقع ، يكفي عشرين سطرًا من التعليمات البرمجية للحل. إذا كنت تحسب مع كل شيء ، كل شيء ، بما في ذلك المحلل اللغوي لملف النموذج ثلاثي الأبعاد ، فاحتفظ به في مائتي سطر - فقط ابصق.

مثال 1: تنعيم البيانات

لنتحدث عن كيفية عملها. لنبدأ من بعيد ، تخيل أن لدينا مصفوفة منتظمة f ، على سبيل المثال ، من 32 عنصرًا ، تمت تهيئتها على النحو التالي:



ثم سنقوم بتنفيذ الإجراء التالي ألف مرة: لكل خلية f [i] نكتب متوسط ​​قيمة الخلايا المجاورة فيها: f [i] = (f [i-1] + f [i + 1]) / 2. لتوضيح الأمر ، إليك الرمز الكامل:

import matplotlib.pyplot as plt f = [0.] * 32 f[0] = f[-1] = 1. f[18] = 2. plt.plot(f, drawstyle='steps-mid') for iter in range(1000): f[0] = f[1] for i in range(1, len(f)-1): f[i] = (f[i-1]+f[i+1])/2. f[-1] = f[-2] plt.plot(f, drawstyle='steps-mid') plt.show() 

كل تكرار سوف يعمل على تنعيم بيانات الصفيف ، وبعد ألف تكرار سنحصل على قيمة ثابتة في جميع الخلايا. إليك رسم متحرك لأول مائة وخمسين تكرارًا:



إذا لم يكن من الواضح لك سبب حدوث التعرج ، فتوقف الآن ، وامسك بقلم وحاول استخلاص الأمثلة ، وإلا فلن يكون من المنطقي قراءة المزيد. السطح المثلث هو في الأساس نفس هذا المثال. تخيل أننا في كل قمة نجد جيرانها ، نحسب مركز كتلتهم ، وننقل رؤوسنا إلى مركز الكتلة ، وهكذا عشر مرات. ستكون النتيجة هكذا:



بالطبع ، إذا لم تتوقف عند عشرة تكرارات ، فبعد فترة سيتم ضغط السطح بأكمله إلى نقطة واحدة بالطريقة نفسها تمامًا كما في المثال السابق ، تم ملء المصفوفة بالكامل بنفس القيمة.

مثال 2: تحسين / تخفيف الميزات

الكود الكامل متاح على github ، ولكن هنا سأعطي الجزء الأكثر أهمية ، مع حذف قراءة وكتابة النماذج ثلاثية الأبعاد فقط. لذا ، فإن النموذج المثلث بالنسبة لي يتمثل في مصفوفتين: علامات ووجوه. إن مجموعة verts هي مجرد مجموعة من النقاط ثلاثية الأبعاد ، وهي رؤوس شبكة المضلع. صفيف الوجوه عبارة عن مجموعة من المثلثات (عدد المثلثات يساوي faces.size ()) ، لكل مثلث يتم تخزين المؤشرات من صفيف قمة الرأس في الصفيف. وصفت تنسيق البيانات والعمل معها بالتفصيل في مسار محاضراتي على رسومات الكمبيوتر. هناك أيضًا مصفوفة ثالثة ، أرويها من الأولين (بشكل أكثر دقة ، فقط من مصفوفة الوجوه) - vvadj. هذا هو الصفيف الذي لكل رأس (الفهرس الأول لصفيف ثنائي الأبعاد) يخزن مؤشرات جيرانه (الفهرس الثاني).

 std::vector<Vec3f> verts; std::vector<std::vector<int> > faces; std::vector<std::vector<int> > vvadj; 

أول شيء أفعله هو بالنسبة لكل قمة من سطحي ، فأنا أعتبر متجه الانحناء. دعونا نوضح: بالنسبة للقمة الحالية v ، أقوم بالتكرار على جميع جيرانها n1-n4 ؛ ثم أحسب مركز كتلتها b = (n1 + n2 + n3 + n4) / 4. حسنًا ، يمكن حساب متجه الانحناء النهائي على أنه c = vb ، فهو ليس مثل الاختلافات المحدودة العادية للمشتقة الثانية .



في الشفرة مباشرة ، يبدو كما يلي:

  std::vector<Vec3f> curvature(verts.size(), Vec3f(0,0,0)); for (int i=0; i<(int)verts.size(); i++) { for (int j=0; j<(int)vvadj[i].size(); j++) { curvature[i] = curvature[i] - verts[vvadj[i][j]]; } curvature[i] = verts[i] + curvature[i] / (float)vvadj[i].size(); } 

حسنًا ، ثم نقوم بالشيء التالي عدة مرات (انظر الصورة السابقة): ننتقل قمة v إلى v: = b + const * c. يرجى ملاحظة أنه إذا كان الثابت يساوي واحدًا ، فلن يتحرك رأسنا إلى أي مكان! إذا كان الثابت صفرًا ، فسيتم استبدال الرأس بمركز كتلة القمم المجاورة ، والتي ستنعّم سطحنا. إذا كان الثابت أكبر من الوحدة (تم عمل صورة العنوان باستخدام const = 2.1) ، فإن الرأس سيتحول في اتجاه متجه انحناء السطح ، مما يعزز التفاصيل. هكذا تبدو في الكود:

  for (int it=0; it<100; it++) { for (int i=0; i<(int)verts.size(); i++) { Vec3f bary(0,0,0); for (int j=0; j<(int)vvadj[i].size(); j++) { bary = bary + verts[vvadj[i][j]]; } bary = bary / (float)vvadj[i].size(); verts[i] = bary + curvature[i]*2.1; // play with the coeff here } } 

بالمناسبة ، إذا كانت أقل من الوحدة ، فإن التفاصيل ستضعف على العكس (const = 0.5) ، لكن هذا لن يكون مكافئًا للتنعيم البسيط ، سيبقى "تباين الصورة":



يرجى ملاحظة أن الكود الخاص بي يولد ملف نموذج ثلاثي الأبعاد بتنسيق Wavefront .obj ، لقد قمت بتقديمه في برنامج تابع لجهة خارجية. يمكنك مشاهدة النموذج الناتج ، على سبيل المثال ، في العارض عبر الإنترنت . إذا كنت مهتمًا بطرق العرض ، بدلاً من إنشاء النموذج ، فقم بقراءة الدورة التدريبية الخاصة بي على رسومات الكمبيوتر .

مثال 3: إضافة قيود

دعنا نرجع إلى المثال الأول ، ونفعل نفس الشيء بالضبط ، ولكن لن نعيد كتابة عناصر المصفوفة تحت الأرقام 0 و 18 و 31:

 import matplotlib.pyplot as plt x = [0.] * 32 x[0] = x[31] = 1. x[18] = 2. plt.plot(x, drawstyle='steps-mid') for iter in range(1000): for i in range(len(x)): if i in [0,18,31]: continue x[i] = (x[i-1]+x[i+1])/2. plt.plot(x, drawstyle='steps-mid') plt.show() 

لقد قمت بتهيئة العناصر المتبقية "المجانية" من المصفوفة بالأصفار ، وما زلت استبدليها بشكل متكرر بمتوسط ​​قيمة العناصر المجاورة. هذه هي الطريقة التي يبدو بها تطور الصفيف في أول مائة وخمسين تكرارًا:



من الواضح تمامًا أن الحل هذه المرة لن يتقارب إلى عنصر ثابت يملأ المصفوفة ، بل إلى منحدرين خطيين. بالمناسبة ، هل هذا واضح للجميع؟ إذا لم يكن الأمر كذلك ، فقم بتجربة هذا الرمز ، وأعطيك أمثلة على وجه التحديد برمز قصير جدًا حتى تتمكن من فهم ما يحدث بدقة.

الاستطراد الغنائي: الحل العددي لأنظمة المعادلات الخطية.

دعونا نعطي النظام المعتاد للمعادلات الخطية:



يمكن إعادة كتابته تاركًا في كل من المعادلات على جانب واحد من علامة المساواة x:



دعونا نعطي ناقلات تعسفية تقريب حل نظام المعادلات (على سبيل المثال ، صفر).

بعد ذلك ، عند التمسك به في الجانب الأيمن من نظامنا ، يمكننا الحصول على متجه تقريبي محدث للحلول .

لجعل الأمر أكثر وضوحًا ، يتم الحصول على x1 من x0 على النحو التالي:



تكرار العملية k مرات ، سيتم تقريب الحل بواسطة المتجه

دعنا فقط في حالة ، اكتب صيغة العودية:



تحت بعض الافتراضات حول معاملات النظام (على سبيل المثال ، من الواضح أن العناصر القطرية لا يجب أن تكون صفراً ، لأننا نقسمها) ، يتحول هذا الإجراء إلى الحل الحقيقي. هذه الجمباز معروفة في الأدب باسم طريقة جاكوبي . بالطبع ، هناك طرق أخرى لحل أنظمة المعادلات الخطية عدديًا ، وأكثر قوة بكثير ، على سبيل المثال ، طريقة التدرج المترافق ، ولكن ، ربما ، تعد طريقة جاكوبي واحدة من أبسط الطرق.

المثال 3 مرة أخرى ، ولكن من ناحية أخرى

والآن دعونا نلقي نظرة فاحصة على الحلقة الرئيسية من المثال 3:

 for iter in range(1000): for i in range(len(x)): if i in [0,18,31]: continue x[i] = (x[i-1]+x[i+1])/2. 

لقد بدأت ببعض المتجه الأولي x ، وقمت بتحديثه ألف مرة ، وإجراء التحديث مشابه بشكل مثير للريبة لطريقة Jacobi! دعونا نكتب نظام المعادلات هذا بشكل صريح:



خذ بعض الوقت ، تأكد من أن كل تكرار في رمز Python الخاص بي هو تحديث واحد فقط لطريقة Jacobi لنظام المعادلات هذا. القيم x [0] و x [18] و x [31] ثابتة بالنسبة لي ، على التوالي ، لا يتم تضمينها في مجموعة المتغيرات ، وبالتالي يتم نقلها إلى الجانب الأيمن.

في المجموع ، تبدو جميع المعادلات في نظامنا - x [i-1] + 2 x [i] - x [i + 1] = 0. هذا التعبير ليس أكثر من اختلافات محدودة عادية للمشتق الثاني . أي أن نظام المعادلات الخاص بنا ينص ببساطة على أن المشتق الثاني يجب أن يكون مساوياً للصفر في كل مكان (حسناً ، باستثناء النقطة x [18]). تذكر ، قلت أنه من الواضح أن التكرارات يجب أن تتلاقى إلى منحدرات خطية؟ هذا هو بالضبط السبب في أن المشتق الخطي للمشتق الثاني هو صفر.

هل تعلم أننا قمنا للتو بحل مشكلة Dirichlet لمعادلة لابلاس ؟

بالمناسبة ، كان يجب أن يلاحظ القارئ المتنبه أنه ، بشكل صارم ، في شفرتي ، يتم حل أنظمة المعادلات الخطية ليس من خلال طريقة جاكوبي ، ولكن من خلال طريقة غاوس-سيدل ، وهو نوع من تحسين طريقة جاكوبي:

مثال 4: معادلة بواسون

ودعنا نغير المثال الثالث قليلاً: كل خلية لا توضع فقط في مركز كتلة الخلايا المجاورة ، ولكن في وسط الكتلة بالإضافة إلى بعض الثابت (التعسفي):

 import matplotlib.pyplot as plt x = [0.] * 32 x[0] = x[31] = 1. x[18] = 2. plt.plot(x, drawstyle='steps-mid') for iter in range(1000): for i in range(len(x)): if i in [0,18,31]: continue x[i] = (x[i-1]+x[i+1] +11./32**2 )/2. plt.plot(x, drawstyle='steps-mid') plt.show() 

في المثال السابق ، وجدنا أن وضعه في مركز الكتلة هو تقدير لعامل لابلاس (في حالتنا ، المشتق الثاني). أي أننا الآن بصدد حل نظام من المعادلات التي تنص على أنه يجب أن يكون لإشارتنا بعض المشتقات الثانية الثابتة. المشتق الثاني هو انحناء السطح. وبالتالي ، ينبغي أن تكون الوظيفة التربيعية التدريجية هي الحل لنظامنا. دعونا نتحقق من أخذ العينات في 32 عينة:



مع صفيف طوله 32 عنصرًا ، يتقارب نظامنا إلى حل في بضع مئات من التكرارات. ولكن ماذا لو جربنا مجموعة من 128 عنصرًا؟ كل شيء أكثر حزنًا هنا ، يجب أن يحسب عدد التكرار بالفعل بالآلاف:



طريقة Gauss-Seidel بسيطة للغاية للبرمجة ، ولكنها لا تنطبق على أنظمة المعادلات الكبيرة. يمكنك محاولة تسريعها ، باستخدام ، على سبيل المثال ، طرق متعددة . بالكلمات ، قد يبدو هذا مرهقًا ، لكن الفكرة بدائية للغاية: إذا أردنا حلًا يحتوي على ألف عنصر ، فيمكننا أولاً حلها بعشرة عناصر ، والحصول على تقريب تقريبي ، ثم مضاعفة الدقة ، وحلها مرة أخرى ، وما إلى ذلك ، حتى نصل إلى النتيجة المرجوة. في الممارسة العملية ، يبدو هذا:



ولا يمكنك التباهي ، واستخدام المحلل الحقيقي لأنظمة المعادلات. لذا ، أحل نفس المثال عن طريق إنشاء المصفوفة A والعمود b ، ثم حل معادلة المصفوفة Ax = b:

 import numpy as np import matplotlib.pyplot as plt n=1000 x = [0.] * n x[0] = x[-1] = 1. m = n*57//100 x[m] = 2. A = np.matrix(np.zeros((n, n))) for i in range(1,n-2): A[i, i-1] = -1. A[i, i] = 2. A[i, i+1] = -1. A = A[1:-2,1:-2] A[m-2,m-1] = 0 A[m-1,m-2] = 0 b = np.matrix(np.zeros((n-3, 1))) b[0,0] = x[0] b[m-2,0] = x[m] b[m-1,0] = x[m] b[-1,0] = x[-1] for i in range(n-3): b[i,0] += 11./n**2 x2 = ((np.linalg.inv(A)*b).transpose().tolist()[0]) x2.insert(0, x[0]) x2.insert(m, x[m]) x2.append(x[-1]) plt.plot(x2, drawstyle='steps-mid') plt.show() 

وهنا نتيجة عمل هذا البرنامج ، لاحظ أن الحل تحول على الفور:



وهكذا ، في الواقع ، وظيفتنا هي تربيعية مجزأة (المشتق الثاني ثابت). في المثال الأول ، قمنا بتعيين المشتق الثاني الثاني على صفر ، والثالث على غير الصفر ، ولكن نفس الشيء في كل مكان. وماذا كان في المثال الثاني؟ قمنا بحل معادلة بواسون المنفصلة بتحديد انحناء السطح. دعني أذكرك بما حدث: لقد حسبنا انحناء السطح الوارد. إذا قمنا بحل مشكلة بواسون عن طريق تعيين انحناء السطح عند الإخراج يساوي انحناء السطح عند الإدخال (const = 1) ، فلن يتغير شيء. يحدث تقوية ملامح الوجه عندما نقوم ببساطة بزيادة الانحناء (الثابت = 2.1). وإذا كان const <1 ، فإن انحناء السطح الناتج ينخفض.

تحديث: لعبة أخرى مثل الواجبات المنزلية

تطوير الفكرة التي اقترحها SquareRootOfZero ، العب مع هذا الرمز:



هذه هي النتيجة الافتراضية ، لينين الأحمر هو البيانات الأولية ، المنحنى الأزرق هو تطورهم ، في اللانهاية تتقارب النتيجة إلى نقطة:



وإليكم النتيجة بالمعامل 2.:



الواجب المنزلي: لماذا في الحالة الثانية ، يتحول لينين أولاً إلى دزيرجينسكي ، ثم يتحول مرة أخرى إلى لينين ، ولكن بحجم أكبر؟

الخلاصة

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

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



ترقبوا المزيد!