حل دستگاه‌های معادلات خطی و تجزیه ماتریس‌ها در پایتون

در یادگیری ماشین و علم داده، کارهای روزمره‌ای مثل آموزش رگرسیون خطی، برازش کمترین مربعات و حتی انتشار رو به جلو در شبکه‌های عصبی همگی به حل معادله‌ی پایه‌ای دستگاه‌های خطی وابسته‌اند. فهم سازوکار این دستگاه‌ها تفاوت بین یک کاربر کتابخانه‌های آماده و یک مهندس مسلط به بهینه‌سازی محاسباتی را مشخص می‌کند.

در این درس تمام روش‌های کلیدی حل معادلات خطی از جمله حذف گاوسی، تجزیه‌های LU، QR و چولسکی و همچنین روش‌های تکراری مانند گرادیان مزدوج را از صفر پیاده‌سازی می‌کنید. یاد می‌گیرید که چطور عدد حالت پایداری محاسبات ممیز شناور را رقم می‌زند و چگونه منظم‌سازی مانع از فروپاشی دقت الگوریتم‌ها می‌شود.

کدهای این آموزش متکی بر ابزارهای استاندارد پایتون مانند NumPy هستند و برای اجرا نیازی به سخت‌افزار خاص یا پردازش ابری ندارید؛ بنابراین به‌راحتی و بدون چالش روی رایانه‌ی شخصی قابل اجرا و آزمایش هستند.

حل معادله‌ی $Ax = b$ قدیمی‌ترین مسئله‌ی ریاضیات است که هنوز هم بار محاسباتی شبکه‌های عصبی شما را به دوش می‌کشد.

اهداف یادگیری

• حل $Ax = b$ با استفاده از حذف گاوسی به همراه محورگیری جزئی (partial pivoting) و جانشینی پس‌رو (back substitution)
• تجزیه‌ی ماتریس‌ها با روش‌های LU، QR و چولسکی (Cholesky) و توضیح کاربرد مناسب هر یک
• استخراج معادلات نرمال (normal equations) برای روش کمترین مربعات و پیوند دادن آن‌ها به رگرسیون خطی و رگرسیون ستیغی (ridge regression)
• عیب‌یابی دستگاه‌های بدحالت (ill-conditioned) با عدد حالت (condition number) و پایدارسازی آن‌ها به کمک منظم‌سازی (regularization)

صورت مسئله

هر بار که یک رگرسیون خطی را آموزش می‌دهید، در حال حل یک دستگاه خطی هستید. هر بار که برازش کمترین مربعات (least-squares fit) را محاسبه می‌کنید، یک دستگاه خطی حل می‌کنید. هر بار که لایه‌ای از یک شبکه عصبی عبارت y = Wx + b را محاسبه می‌کند، در واقع یک سمت از یک دستگاه خطی را ارزیابی می‌کند. وقتی منظم‌سازی اضافه می‌کنید، ساختار دستگاه را تغییر می‌دهید. وقتی از فرایندهای گاوسی استفاده می‌کنید، ماتریسی را تجزیه می‌کنید. وقتی ماتریس کوواریانس را برای محاسبه‌ی فاصله‌ی ماهالانوبیس (Mahalanobis distance) معکوس می‌کنید، باز هم دارید یک دستگاه خطی حل می‌کنید.

معادله‌ی $Ax = b$ همه‌جا به چشم می‌خورد. در این رابطه، $A$ ماتریس ضرایب معلوم، $b$ بردار خروجی‌های معلوم، و $x$ بردار مجهولاتی است که می‌خواهید بیابید. در رگرسیون خطی، $A$ ماتریس داده‌ها، $b$ بردار متغیر هدف و $x$ بردار وزن‌هاست. کل مدل در نهایت به همین مسئله خلاصه می‌شود: بردار $x$ را طوری پیدا کنید که $Ax$ تا جای ممکن به $b$ نزدیک باشد.

در این درس، تمام روش‌های اصلی حل این معادله را از صفر پیاده‌سازی می‌کنیم. متوجه خواهید شد چرا برخی روش‌ها سریع و برخی دیگر پایدارند، چرا دسته‌ای فقط برای دستگاه‌های مربعی جواب می‌دهند و دسته‌ای دیگر دستگاه‌های فرامعین (overdetermined) را هم مدیریت می‌کنند، و چرا عدد حالت ماتریس مشخص می‌کند که آیا اصلاً می‌توان به جواب حاصل اعتماد کرد یا نه.

مفاهیم پایه

معنای هندسی Ax = b

هر دستگاه معادلات خطی تفسیری هندسی دارد. هر معادله معرف یک ابرصفحه (hyperplane) است. پاسخ مسئله، نقطه (یا مجموعه‌نقاطی) است که همه‌ی این ابرصفحه‌ها در آن یکدیگر را قطع می‌کنند.

2x + y = 5          Two lines in 2D.
x - y  = 1          They intersect at x=2, y=1.

سه حالت ممکن است رخ دهد:

در قالب ماتریسی، «یک جواب» یعنی ماتریس $A$ وارون‌پذیر است. «بدون جواب» یعنی دستگاه ناسازگار است. «بینهایت جواب» یعنی ماتریس $A$ فضای پوچ (null space) دارد. اغلب مسائل یادگیری ماشین در دسته‌ی «بدون جواب دقیق» جای می‌گیرند؛ زیرا تعداد معادلات (نقاط داده) بیشتر از مجهولات (پارامترها) است. اینجاست که روش کمترین مربعات وارد میدان می‌شود.

دیدگاه سطری در برابر دیدگاه ستونی

معادله‌ی $Ax = b$ را می‌توان به دو شکل تفسیر کرد:

دیدگاه سطری. هر سطر $A$ یک معادله می‌سازد. هر معادله نشان‌دهنده‌ی یک ابرصفحه است. جواب مسئله نقطه‌ی تلاقی همه‌ی این ابرصفحه‌هاست.

دیدگاه ستونی. هر ستون $A$ یک بردار است. در این حالت مسئله این است: چه ترکیب خطی‌ای از ستون‌های $A$ بردار $b$ را می‌سازد؟

A = | 2  1 |    b = | 5 |
    | 1 -1 |        | 1 |

Row picture: solve 2x + y = 5 and x - y = 1 simultaneously.

Column picture: find x1, x2 such that:
  x1 * [2, 1] + x2 * [1, -1] = [5, 1]
  2 * [2, 1] + 1 * [1, -1] = [4+1, 2-1] = [5, 1]   check.

دیدگاه ستونی بنیادی‌تر است. اگر $b$ در فضای ستونیِ $A$ جای بگیرد، دستگاه جواب دارد. اگر $b$ درون این فضا نباشد، نزدیک‌ترین نقطه در فضای ستونی را پیدا می‌کنید. آن نزدیک‌ترین نقطه، همان جواب کمترین مربعات است.

حذف گاوسی

روش حذف گاوسی، معادله‌ی $Ax = b$ را به یک دستگاه بالامثلثی $Ux = c$ تبدیل می‌کند که می‌توان آن را با جانشینی پس‌رو حل کرد. این مستقیم‌ترین شیوه‌ی حل مسئله است.

مراحل الگوریتم:

1. For each column k (the pivot column):
   a. Find the largest entry in column k at or below row k (partial pivoting).
   b. Swap that row with row k.
   c. For each row i below k:
      - Compute multiplier m = A[i][k] / A[k][k]
      - Subtract m times row k from row i.
2. Back substitute: solve from the last equation upward.

مثال:

Original:
| 2  1  1 | 8 |       R2 = R2 - (2)R1     | 2  1   1 |  8 |
| 4  3  3 |20 |  -->  R3 = R3 - (1)R1 --> | 0  1   1 |  4 |
| 2  3  1 |12 |                            | 0  2   0 |  4 |

                       R3 = R3 - (2)R2     | 2  1   1 |  8 |
                                       --> | 0  1   1 |  4 |
                                           | 0  0  -2 | -4 |

Back substitute:
  -2 * x3 = -4    -->  x3 = 2
  x2 + 2  = 4     -->  x2 = 2
  2*x1 + 2 + 2 = 8 --> x1 = 2

هزینه‌ی محاسباتی حذف گاوسی برابر با $O(n^3)$ عملیات است. برای دستگاهی با ابعاد $1000 \times 1000$، این مقدار معادل حدود یک میلیارد عملیات ممیز شناور است؛ مقداری سریع، اما اگر نیاز داشته باشید چندین دستگاه با ماتریس $A$ یکسان را حل کنید، راهکارهای کارآمدتری هم وجود دارد.

محورگیری جزئی: چرا اهمیت دارد؟

بدون محورگیری (pivoting)، حذف گاوسی ممکن است شکست بخورد یا خروجی کاملاً بی‌ارزشی تولید کند. اگر عنصر لولا صفر باشد، با خطای تقسیم بر صفر مواجه می‌شوید؛ اگر خیلی کوچک باشد، خطاهای گرد کردن به‌شدت تقویت می‌شوند.

Bad pivot:                       With partial pivoting:
| 0.001  1 | 1.001 |            Swap rows first:
| 1      1 | 2     |            | 1      1 | 2     |
                                 | 0.001  1 | 1.001 |
m = 1/0.001 = 1000              m = 0.001/1 = 0.001
R2 = R2 - 1000*R1               R2 = R2 - 0.001*R1
| 0.001  1     | 1.001   |      | 1      1     | 2     |
| 0     -999   | -999.0  |      | 0      0.999 | 0.999 |

x2 = 1.000 (correct)            x2 = 1.000 (correct)
x1 = (1.001 - 1)/0.001          x1 = (2 - 1)/1 = 1.000 (correct)
   = 0.001/0.001 = 1.000        Stable because the multiplier is small.

در محاسبات ممیز شناور با دقت محدود، نسخه‌ی بدون محورگیری ممکن است ارقام بااهمیت خود را از دست بدهد. محورگیری جزئی همواره بزرگ‌ترین عنصر در دسترس را به عنوان لولا انتخاب می‌کند تا تقویت خطا به حداقل برسد.

تجزیه‌ی LU

تجزیه‌ی LU ماتریس $A$ را به حاصل‌ضرب یک ماتریس پایین‌مثلثی $L$ و یک ماتریس بالامثلثی $U$ تجزیه می‌کند: $A = LU$. ماتریس $L$ ضرایب حاصل از حذف گاوسی را در خود ذخیره می‌کند و ماتریس $U$ نتیجه‌ی عملیات حذف است.

A = L @ U

| 2  1  1 |   | 1  0  0 |   | 2  1   1 |
| 4  3  3 | = | 2  1  0 | @ | 0  1   1 |
| 2  3  1 |   | 1  2  1 |   | 0  0  -2 |

چرا به‌جای حذف مستقیم، سراغ فاکتورگیری و تجزیه برویم؟ چون پس از به‌دست آوردن $L$ و $U$، حل معادله‌ی $Ax = b$ برای هر بردار جدیدِ $b$ فقط هزینه‌ای معادل $O(n^2)$ خواهد داشت:

$$\begin{aligned} Ax &= b \\ LUx &= b \\ &\text{Let } y = Ux\text{:} \\ Ly &= b \quad \text{(forward substitution, } O(n^2)\text{)} \\ Ux &= y \quad \text{(back substitution, } O(n^2)\text{)} \end{aligned}$$

هزینه‌ی سنگین $O(n^3)$ فقط یک بار هنگام تجزیه پرداخت می‌شود و هر حلِ پس از آن فقط به اندازه‌ی $O(n^2)$ هزینه دارد. اگر لازم باشد 1000 دستگاه را با ماتریس $A$ یکسان اما بردارهای $b$ متفاوت حل کنید، تجزیه‌ی LU بار کل محاسبات را به نسبت 1000/3 کاهش می‌دهد.

اگر محورگیری جزئی را هم اعمال کنید، رابطه به صورت $PA = LU$ درمی‌آید که در آن $P$ ماتریس جایگشتی است که جابه‌جایی سطرها را ثبت می‌کند.

تجزیه‌ی QR

تجزیه‌ی QR ماتریس $A$ را به یک ماتریس متعامد (orthogonal) به نام $Q$ و یک ماتریس بالامثلثی $R$ تجزیه می‌کند: $A = QR$.

ماتریس متعامد ویژگی $Q^T Q = I$ را دارد. ستون‌های آن بردار‌هایی یکه‌متعامد (orthonormal) هستند. ضرب در $Q$ طول بردارها و زاویه‌ی بین آن‌ها را حفظ می‌کند.

A = Q @ R

Q has orthonormal columns: Q^T Q = I
R is upper triangular

To solve Ax = b:
  QRx = b
  Rx = Q^T b    (just multiply by Q^T, no inversion needed)
  Back substitute to get x.

تجزیه‌ی QR برای حل مسائل کمترین مربعات، از نظر عددی پایدارتر از LU عمل می‌کند. فرایند گرام-اشمیت (Gram-Schmidt) ماتریس $Q$ را ستون‌به‌ستون می‌سازد:

Given columns a1, a2, ... of A:

q1 = a1 / ||a1||

q2 = a2 - (a2 . q1) * q1        (subtract projection onto q1)
q2 = q2 / ||q2||                (normalize)

q3 = a3 - (a3 . q1) * q1 - (a3 . q2) * q2
q3 = q3 / ||q3||

R[i][j] = qi . aj    for i <= j

در هر مرحله، مؤلفه‌ی بردار در امتداد تمام بردارهای قبلیِ $q$ حذف می‌شود و تنها جهت متعامد جدید باقی می‌ماند.

تجزیه‌ی چولسکی

وقتی $A$ متقارن ($A = A^T$) و معین مثبت باشد (یعنی تمام مقادیر ویژه‌ی آن مثبت باشند)، می‌توان آن را به‌صورت $A = L L^T$ تجزیه کرد که در آن $L$ یک ماتریس پایین‌مثلثی است. به این فرایند، تجزیه‌ی چولسکی می‌گویند.

A = L @ L^T

| 4  2 |   | 2  0 |   | 2  1 |
| 2  5 | = | 1  2 | @ | 0  2 |

L[i][i] = sqrt(A[i][i] - sum(L[i][k]^2 for k < i))
L[i][j] = (A[i][j] - sum(L[i][k]*L[j][k] for k < j)) / L[j][j]    for i > j

چولسکی دو برابر سریع‌تر از LU است و نصف آن فضا می‌گیرد. این روش فقط برای ماتریس‌های متقارنِ معین مثبت کاربرد دارد، اما این ماتریس‌ها مدام در یادگیری ماشین ظاهر می‌شوند:

• ماتریس‌های کوواریانس متقارن و نیمه‌معین مثبت هستند (و با افزودن منظم‌سازی، معین مثبت می‌شوند).
• ماتریس کرنل (kernel matrix) در فرایندهای گاوسی متقارن و معین مثبت است.
• ماتریس هسیان (Hessian) یک تابع محدب در نقطه‌ی کمینه، متقارن و معین مثبت است.
• عبارت $A^T A$ همیشه متقارن و نیمه‌معین مثبت است.

در فرایندهای گاوسی، ماتریس کرنل $K$ را با روش چولسکی تجزیه می‌کنید و سپس برای یافتن میانگین پیش‌بین، دستگاه $K \alpha = y$ را حل می‌نمایید. فاکتور چولسکی همچنین مقدار لگاریتم دترمینان مورد نیاز برای محاسبه‌ی درست‌نمایی حاشیه‌ای را به دست می‌دهد: $\log \det(K) = 2 \sum \log(\mathop{\mathrm{diag}}(L))$.

کمترین مربعات: وقتی Ax = b جواب دقیقی ندارد

اگر $A$ ماتریسی $m \times n$ با شرط $m > n$ باشد (تعداد معادلات بیشتر از مجهولات)، با دستگاهی فرامعین روبه‌رو هستیم و هیچ جواب دقیقی وجود ندارد. در چنین وضعیتی، به‌جای حل دقیق، مجموع مربعات خطا را کمینه می‌کنیم:

minimize ||Ax - b||^2

This is the sum of squared residuals:
  sum((A[i,:] @ x - b[i])^2 for i in range(m))

نقطه‌ی کمینه در معادلات نرمال صدق می‌کند:

$$A^T A x = A^T b$$

اثبات: عبارت $\lVert Ax - b \rVert^2 = (Ax - b)^T (Ax - b) = x^T A^T A x - 2 x^T A^T b + b^T b$ را بسط دهید، گرادیان آن را نسبت به $x$ حساب کنید و برابر صفر بگذارید: $2 A^T A x - 2 A^T b = 0$.

Original system (overdetermined, 4 equations, 2 unknowns):
| 1  1 |         | 3 |
| 1  2 | x     = | 5 |       No exact x satisfies all 4 equations.
| 1  3 |         | 6 |
| 1  4 |         | 8 |

Normal equations:
A^T A = | 4  10 |    A^T b = | 22 |
        | 10 30 |            | 63 |

Solve: x = [1.5, 1.7]

This is linear regression. x[0] is the intercept, x[1] is the slope.

معادلات نرمال = رگرسیون خطی

این تطابق کاملاً دقیق و مستقیم است. در رگرسیون خطی، ماتریس داده‌های شما یعنی $X$ به‌ازای هر نمونه یک سطر و به‌ازای هر ویژگی یک ستون دارد. بردار متغیر هدف $y$ نیز برای هر نمونه یک درایه دارد. بردار وزن‌های $w$ در معادله‌ی زیر صدق می‌کند:

$$\begin{aligned} X^T X w &= X^T y \\ w &= (X^T X)^{-1} X^T y \end{aligned}$$

این همان جواب فرم بسته (closed-form solution) برای رگرسیون خطی است. هر بار که متد fit() را در sklearn.linear_model.LinearRegression اجرا می‌کنید، دقیقاً همین رابطه (یا معادل آن از طریق QR یا SVD) محاسبه می‌شود.

اگر یک جمله‌ی منظم‌سازیِ $\lambda I$ به این ماتریس بیفزایید، به رگرسیون ستیغی می‌رسید:

$$\begin{aligned} (X^T X + \lambda I) w &= X^T y \\ w &= (X^T X + \lambda I)^{-1} X^T y \end{aligned}$$

منظم‌سازی وضعیت ماتریس را بهبود می‌دهد (معکوس کردن دقیق آن را آسان‌تر می‌کند) و با جمع کردن وزن‌ها به سمت صفر، جلوی بیش‌برازش (overfitting) را می‌گیرد. ماتریس $X^T X + \lambda I$ به‌ازای $\lambda > 0$ همیشه متقارن و معین مثبت است؛ در نتیجه می‌توانید برای حل آن از تجزیه‌ی چولسکی استفاده کنید.

شبه‌وارون (مور-پنروز)

شبه‌وارون (pseudoinverse) که با $A^+$ نمایش داده می‌شود، مفهوم معکوس ماتریس را برای ماتریس‌های غیرمربعی و منفرد (singular) تعمیم می‌دهد. برای هر ماتریس $A$:

$$\begin{aligned} x &= A^+ b \\[4pt] \text{where } A^+ &= V \Sigma^+ U^T \quad \text{(computed via SVD)} \end{aligned}$$

ماتریس $\Sigma^+$ با وارون کردن تک‌تک مقادیر منفرد غیرصفر و سپس ترانهاده کردن حاصل به دست می‌آید. اگر $A = U \Sigma V^T$ باشد، آنگاه $A^+ = V \Sigma^+ U^T$.

A = U Sigma V^T        (SVD)

Sigma = | 5  0 |       Sigma+ = | 1/5  0  0 |
        | 0  2 |                | 0  1/2  0 |
        | 0  0 |

A+ = V Sigma+ U^T

شبه‌وارون، جوابی از مسئله‌ی کمترین مربعات را تحویل می‌دهد که کمترین نُرم را دارد. اگر دستگاه:

• یک جواب داشته باشد: $A^+ b$ همان جواب را به دست می‌دهد.
• جوابی نداشته باشد: $A^+ b$ جواب کمترین مربعات را تولید می‌کند.
• بینهایت جواب داشته باشد: $A^+ b$ جوابی را تعیین می‌کند که کمترین $\lVert x \rVert$ را دارد.

توابع np.linalg.lstsq و np.linalg.pinv در NumPy در پشت صحنه از تجزیه‌ی مقادیر منفرد (SVD) استفاده می‌کنند.

عدد حالت

عدد حالت میزان حساسیت جواب دستگاه به تغییرات ناچیز در ورودی را نشان می‌دهد. برای ماتریس $A$، عدد حالت به‌صورت زیر تعریف می‌شود:

$$\kappa(A) = \lVert A \rVert \cdot \lVert A^{-1} \rVert = \frac{\sigma_{\max}}{\sigma_{\min}}$$

که در آن $\sigma_{\max}$ و $\sigma_{\min}$ بزرگ‌ترین و کوچک‌ترین مقادیر منفرد هستند.

Well-conditioned (kappa ~ 1):        Ill-conditioned (kappa ~ 10^15):
Small change in b -->                Small change in b -->
small change in x                    huge change in x

| 2  0 |   kappa = 2/1 = 2          | 1   1          |   kappa ~ 10^15
| 0  1 |   safe to solve            | 1   1+10^(-15) |   solution is garbage

قواعد سرانگشتی:

• $\kappa < 100$: دستگاه پایدار است و جواب دقت مطلوبی دارد.
• $\kappa \approx 10^k$: حدود $k$ رقم از دقت محاسبات ممیز شناور خود را از دست می‌دهید.
• $\kappa \approx 10^{16}$ (برای داده‌های float64): جواب کاملاً بی‌معنی است و ماتریس عملاً منفرد به شمار می‌رود.

در یادگیری ماشین، بدحالتی معمولاً زمانی پیش می‌آید که ویژگی‌ها تقریباً هم‌خط (collinear) باشند. منظم‌سازی (افزودن $\lambda I$) عدد حالت را از $\sigma_{\max} / \sigma_{\min}$ به $(\sigma_{\max} + \lambda) / (\sigma_{\min} + \lambda)$ ارتقا می‌دهد.

روش‌های تکراری: گرادیان مزدوج

برای دستگاه‌های خلوت (sparse) و بسیار بزرگ با میلیون‌ها مجهول، استفاده از روش‌های مستقیم مانند LU یا چولسکی از نظر بار محاسباتی بسیار سنگین است. روش‌های تکراری، با اصلاح و بهبود تدریجی یک حدس اولیه طی چندین تکرار، جواب را تقریب می‌زنند.

روش گرادیان مزدوج (conjugate gradient یا CG) معادله‌ی $Ax = b$ را در وضعیتی حل می‌کند که $A$ متقارن و معین مثبت باشد. این روش در محاسبات دقیق تئوری، حداکثر ظرف $n$ تکرار به جواب دقیق می‌رسد؛ اما در عمل، اگر مقادیر ویژه‌ی $A$ به صورت خوشه‌ای و نزدیک به هم باشند، بسیار سریع‌تر همگرا خواهد شد.

Algorithm sketch:
  x0 = initial guess (often zero)
  r0 = b - A x0           (residual)
  p0 = r0                 (search direction)

  For k = 0, 1, 2, ...:
    alpha = (rk . rk) / (pk . A pk)
    x_{k+1} = xk + alpha * pk
    r_{k+1} = rk - alpha * A pk
    beta = (r_{k+1} . r_{k+1}) / (rk . rk)
    p_{k+1} = r_{k+1} + beta * pk
    if ||r_{k+1}|| < tolerance: stop

گرادیان مزدوج در این موارد کاربرد دارد:

• بهینه‌سازی در مقیاس بزرگ (روش Newton-CG)
• حل معادلات دیفرانسیل با مشتقات جزئی (PDE) گسسته‌سازی‌شده
• روش‌های کرنل که ماتریس کرنل در آن‌ها بزرگ‌تر از توان تجزیه شدن است
• پیش‌شرط‌سازی (preconditioning) برای سایر حل‌کننده‌های تکراری

سرعت همگرایی وابسته به عدد حالت است. دستگاه‌های خوش‌حالت‌تر سریع‌تر همگرا می‌شوند؛ این نکته دلیل دیگری است بر اینکه چرا منظم‌سازی به کار می‌آید.

جمع‌بندی کلی: کدام روش برای چه کاری مناسب است؟

ارتباط با یادگیری ماشین

تمام روش‌های این درس به شکل ملموس در سیستم‌های واقعی یادگیری ماشین حضور دارند:

رگرسیون خطی. جواب فرم بسته‌ی این مدل، معادلات نرمال $X^T X w = X^T y$ را حل می‌کند. این حل بر حسب شرایط با استفاده از روش چولسکی (اگر $n$ کوچک باشد)، QR (اگر پایداری عددی حیاتی باشد) یا SVD (اگر احتمال نقص رتبه در ماتریس وجود داشته باشد) پیاده‌سازی می‌شود.

رگرسیون ستیغی. مقدار $\lambda I$ را به $X^T X$ اضافه می‌کند. دستگاه منظم‌شده‌ی $(X^T X + \lambda I) w = X^T y$ همیشه با روش چولسکی حل‌شدنی است؛ زیرا $X^T X + \lambda I$ به‌ازای $\lambda > 0$ همواره ماتریسی متقارن و معین مثبت خواهد بود.

فرایندهای گاوسی. محاسبه‌ی میانگین پیش‌بین نیازمند حل دستگاه $K \alpha = y$ است که در آن $K$ ماتریس کرنل به شمار می‌رود. رویکرد استاندارد در این زمینه، تجزیه‌ی چولسکیِ $K$ است. محاسبه‌ی لگاریتم درست‌نمایی حاشیه‌ای نیز از $\log \det(K) = 2 \sum \log(\mathop{\mathrm{diag}}(L))$ استفاده می‌کند.

مقداردهی اولیه‌ی شبکه‌های عصبی. روش مقداردهی اولیه‌ی متعامد (orthogonal initialization) از تجزیه‌ی QR کمک می‌گیرد تا ماتریس‌های وزنی با ستون‌های یکه‌متعامد بسازد. این کار از افت شدید سیگنال در شبکه‌های عمیق جلوگیری می‌کند.

پیش‌شرط‌سازی. بهینه‌سازهای مقیاس‌بزرگ، از چولسکی ناقص یا LU ناقص به عنوان پیش‌شرط‌ساز در حل‌کننده‌های گرادیان مزدوج بهره می‌برند.

مهندسی ویژگی. عدد حالتِ $X^T X$ به شما می‌گوید که آیا ویژگی‌ها هم‌خط هستند یا خیر. اگر مقدار $\kappa$ بزرگ باشد، باید برخی ویژگی‌ها را حذف کنید یا سراغ منظم‌سازی بروید.

پیاده‌سازی از صفر

گام 1: حذف گاوسی با محورگیری جزئی

import numpy as np

def gaussian_elimination(A, b):
    n = len(b)
    Ab = np.hstack([A.astype(float), b.reshape(-1, 1).astype(float)])

    for k in range(n):
        max_row = k + np.argmax(np.abs(Ab[k:, k]))
        Ab[[k, max_row]] = Ab[[max_row, k]]

        if abs(Ab[k, k]) < 1e-12:
            raise ValueError(f"Matrix is singular or nearly singular at pivot {k}")

        for i in range(k + 1, n):
            m = Ab[i, k] / Ab[k, k]
            Ab[i, k:] -= m * Ab[k, k:]

    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (Ab[i, -1] - Ab[i, i+1:n] @ x[i+1:n]) / Ab[i, i]

    return x

گام 2: تجزیه‌ی LU

def lu_decompose(A):
    n = A.shape[0]
    L = np.eye(n)
    U = A.astype(float).copy()
    P = np.eye(n)

    for k in range(n):
        max_row = k + np.argmax(np.abs(U[k:, k]))
        if max_row != k:
            U[[k, max_row]] = U[[max_row, k]]
            P[[k, max_row]] = P[[max_row, k]]
            if k > 0:
                L[[k, max_row], :k] = L[[max_row, k], :k]

        for i in range(k + 1, n):
            L[i, k] = U[i, k] / U[k, k]
            U[i, k:] -= L[i, k] * U[k, k:]

    return P, L, U

def lu_solve(P, L, U, b):
    n = len(b)
    Pb = P @ b.astype(float)

    y = np.zeros(n)
    for i in range(n):
        y[i] = Pb[i] - L[i, :i] @ y[:i]

    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - U[i, i+1:] @ x[i+1:]) / U[i, i]

    return x

گام 3: تجزیه‌ی چولسکی

def cholesky(A):
    n = A.shape[0]
    L = np.zeros_like(A, dtype=float)

    for i in range(n):
        for j in range(i + 1):
            s = A[i, j] - L[i, :j] @ L[j, :j]
            if i == j:
                if s <= 0:
                    raise ValueError("Matrix is not positive definite")
                L[i, j] = np.sqrt(s)
            else:
                L[i, j] = s / L[j, j]

    return L

گام 4: کمترین مربعات از طریق معادلات نرمال

def least_squares_normal(A, b):
    AtA = A.T @ A
    Atb = A.T @ b
    return gaussian_elimination(AtA, Atb)

def ridge_regression(A, b, lam):
    n = A.shape[1]
    AtA = A.T @ A + lam * np.eye(n)
    Atb = A.T @ b
    L = cholesky(AtA)
    y = np.zeros(n)
    for i in range(n):
        y[i] = (Atb[i] - L[i, :i] @ y[:i]) / L[i, i]
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - L.T[i, i+1:] @ x[i+1:]) / L.T[i, i]
    return x

گام 5: عدد حالت

def condition_number(A):
    U, S, Vt = np.linalg.svd(A)
    return S[0] / S[-1]

کاربرد عملی

ترکیب این بخش‌ها برای اجرای رگرسیون خطی و رگرسیون ستیغی روی داده‌های واقعی:

np.random.seed(42)
X_raw = np.random.randn(100, 3)
w_true = np.array([2.0, -1.0, 0.5])
y = X_raw @ w_true + np.random.randn(100) * 0.1

X = np.column_stack([np.ones(100), X_raw])

w_ols = least_squares_normal(X, y)
print(f"OLS weights (ours):    {w_ols}")

w_np = np.linalg.lstsq(X, y, rcond=None)[0]
print(f"OLS weights (numpy):   {w_np}")
print(f"Max difference: {np.max(np.abs(w_ols - w_np)):.2e}")

w_ridge = ridge_regression(X, y, lam=1.0)
print(f"Ridge weights (ours):  {w_ridge}")

from sklearn.linear_model import Ridge
ridge_sk = Ridge(alpha=1.0, fit_intercept=False)
ridge_sk.fit(X, y)
print(f"Ridge weights (sklearn): {ridge_sk.coef_}")

خروجی این درس

این درس نتایج زیر را ارائه می‌دهد:

• فایل code/linear_systems.py شامل پیاده‌سازی پایه‌ای و از صفرِ روش‌های حذف گاوسی، تجزیه‌ی LU، تجزیه‌ی چولسکی، کمترین مربعات و رگرسیون ستیغی
• یک نمایش عملی از اینکه معادلات نرمال و تابع LinearRegression در sklearn وزن‌های کاملاً یکسانی به دست می‌دهند

تمرین‌ها

1. دستگاه [[1,2,3],[4,5,6],[7,8,10]] x = [6, 15, 27] را با تابع حذف گاوسی خود، تابع حل‌کننده‌ی LU خود و تابع np.linalg.solve حل کنید. بررسی کنید که هر سه روش با در نظر گرفتن خطای مجاز ممیز شناور، جواب یکسانی تولید کنند.

2. یک ماتریس تصادفی X با ابعاد 50 در 5 و بردار هدف y = X @ w_true + noise تولید کنید. بردار w را با معادلات نرمال، تجزیه‌ی QR (از طریق np.linalg.qr)، تجزیه‌ی SVD (از طریق np.linalg.svd) و تابع np.linalg.lstsq به دست آورید. هر چهار پاسخ را با هم بسنجید. عدد حالت $X^T X$ را محاسبه کنید و توضیح دهید این عدد چگونه بر میزان اطمینان شما به هر روش اثر می‌گذارد.

3. با ساختن دو ستون تقریباً یکسان، ماتریسی نزدیک به منفرد بسازید (برای نمونه، ستون دوم = ستون اول + 1e-10 * نویز). عدد حالت آن را حساب کنید. دستگاه $Ax = b$ را یک بار بدون منظم‌سازی و یک بار با منظم‌سازی (افزودن $0.01 I$) حل کنید. جواب‌ها و نُرم باقیمانده‌ها را با هم مقایسه کرده و شرح دهید چرا منظم‌سازی به حل مسئله کمک می‌کند.

4. الگوریتم گرادیان مزدوج را برای یک ماتریس تصادفیِ متقارن معین مثبت با ابعاد 100 در 100 پیاده‌سازی کنید. بشمارید که رسیدن به تلرانس 1e-8 به چند تکرار نیاز دارد. نتیجه را با سقف تئوریِ $n$ تکرار مقایسه کنید.

5. زمان اجرای حل‌کننده‌ی چولسکی خود را با حل‌کننده‌ی LU خود و تابع np.linalg.solve روی ماتریس‌های متقارن معین مثبت در ابعاد 10، 50، 200 و 500 مقایسه کنید. نمودار نتایج را رسم کرده و بررسی کنید که آیا چولسکی عملاً حدود 2 برابر سریع‌تر از LU است یا خیر.

راهنمایی‌ها

تمرین 1. ابتدا gaussian_elimination(A, b)، سپس lu_decompose(A) و به دنبال آن lu_solve(P, L, U, b)، و در نهایت np.linalg.solve(A, b) را فراخوانی کنید. هر پاسخ را با دستور np.allclose(A @ x, b) بیازمایید و هر سه نتیجه را با np.allclose با یکدیگر مقایسه کنید.

تمرین 2. از تابع least_squares_normal برای معادلات نرمال و از سه تابع نام‌برده در NumPy برای سایر روش‌ها استفاده کنید. هر چهار جواب باید نزدیک به w_true باشند. مقدار condition_number(X.T @ X) را با condition_number(X) مقایسه کنید: اولی مجذور دومی است؛ بنابراین معادلات نرمال بیشترین افت دقت محاسباتی را خواهند داشت.

تمرین 3. انتظار می‌رود عدد حالت فراتر از $10^{9}$ باشد. دستگاه را یک بار با np.linalg.solve و یک بار با ridge_regression(A, b, 0.01) حل کنید و در هر حالت، مقدار np.linalg.norm(x) و نُرم باقیمانده را مقایسه کنید. باقیمانده‌ی هر دو روش ممکن است کوچک باشد، اما جواب حاصل از روشِ بدون منظم‌سازی معمولاً درایه‌هایی به شدت بزرگ دارد.

تمرین 4. ماتریس متقارن معین مثبت را با دستور B @ B.T + n * np.eye(n) بسازید تا تمام مقادیر ویژه‌ی آن مثبت باشند. یک متغیر شمارنده به حلقه‌ی تابع conjugate_gradient اضافه کنید، سپس بررسی کنید که نُرم باقیمانده‌ی نهایی کمتر از 1e-8 بوده و شمارنده‌ی تکرار از 100 بیشتر نشود.

تمرین 5. برای هر سه حل‌کننده از ماتریس‌های یکسانِ B @ B.T + n * np.eye(n) استفاده کنید و زمان هر حل را با time.perf_counter بسنجید و بهترین زمان را از میان چند اجرا برگزینید. چولسکی تقریباً به نیمی از حجم محاسبات LU نیاز دارد؛ بنابراین این اختلاف سرعت در ابعاد بزرگ‌تر خودش را بهتر نشان می‌دهد.

اصطلاحات کلیدی

برای مطالعه‌ی بیشتر

• MIT 18.06: Linear Algebra (Gilbert Strang) -- دوره‌ی مرجع و معتبر پیرامون دستگاه‌های خطی و تجزیه‌ی ماتریس‌ها
• Numerical Linear Algebra (Trefethen & Bau) -- مرجع استاندارد برای درک عمیق پایداری عددی، عدد حالت و علل شکست الگوریتم‌ها
• Matrix Computations (Golub & Van Loan) -- کتابخانه‌ای جامع و دانشنامه‌ای درباره‌ی تمام الگوریتم‌های ماتریسی
• 3Blue1Brown: Inverse Matrices -- ارائه‌ی شهود بصری از معنای هندسی حل معادله‌ی $Ax = b$