پرش به محتویات

روش گاوس برای حل دستگاه معادلات خطی

با یک دستگاه $n$ معادله جبری خطی (SLAE) با $m$ مجهول، از شما خواسته می‌شود که دستگاه را حل کنید: یعنی مشخص کنید که آیا جوابی ندارد، دقیقاً یک جواب دارد یا بی‌نهایت جواب دارد. و در صورتی که حداقل یک جواب وجود داشته باشد، یکی از آن‌ها را پیدا کنید.

به طور رسمی، مسئله به این صورت فرموله می‌شود: دستگاه زیر را حل کنید:

$$\begin{align} a_{11} x_1 + a_{12} x_2 + &\dots + a_{1m} x_m = b_1 \\ a_{21} x_1 + a_{22} x_2 + &\dots + a_{2m} x_m = b_2\\ &\vdots \\ a_{n1} x_1 + a_{n2} x_2 + &\dots + a_{nm} x_m = b_n \end{align}$$

که در آن ضرایب $a_{ij}$ (برای $i$ از ۱ تا $n$ و $j$ از ۱ تا $m$) و $b_i$ (برای $i$ از ۱ تا $n$) معلوم و متغیرهای $x_i$ (برای $i$ از ۱ تا $m$) مجهول هستند.

این مسئله یک نمایش ماتریسی ساده نیز دارد:

$$Ax = b,$$

که در آن $A$ یک ماتریس با اندازه $n \times m$ از ضرایب $a_{ij}$ و $b$ بردار ستونی با اندازه $n$ است.

شایان ذکر است که روش ارائه شده در این مقاله می‌تواند برای حل معادلات به پیمانه هر عدد p نیز استفاده شود، یعنی:

$$\begin{align} a_{11} x_1 + a_{12} x_2 + &\dots + a_{1m} x_m \equiv b_1 \pmod p \\ a_{21} x_1 + a_{22} x_2 + &\dots + a_{2m} x_m \equiv b_2 \pmod p \\ &\vdots \\ a_{n1} x_1 + a_{n2} x_2 + &\dots + a_{nm} x_m \equiv b_n \pmod p \end{align}$$

گاوس

به بیان دقیق، روشی که در ادامه شرح داده می‌شود باید «گاوس-جردن» یا حذف گاوس-جردن نامیده شود، زیرا این روش گونه‌ای از روش گاوس است که توسط جردن در سال ۱۸۸۷ توصیف شد.

مرور کلی

این الگوریتم یک حذف متوالی متغیرها در هر معادله است، تا زمانی که هر معادله تنها یک متغیر باقی‌مانده داشته باشد. اگر $n=m$ باشد، می‌توانید این فرآیند را به عنوان تبدیل ماتریس $A$ به ماتریس همانی در نظر بگیرید و معادله را در این حالت بدیهی حل کنید، که در آن جواب یکتا و برابر با ضریب $b_i$ است.

حذف گاوسی بر دو تبدیل ساده استوار است:

  • می‌توان دو معادله را با هم جابجا کرد
  • هر معادله را می‌توان با ترکیب خطی آن سطر (با ضریب غیرصفر) و چند سطر دیگر (با ضرایب دلخواه) جایگزین کرد.

در گام اول، الگوریتم گاوس-جردن سطر اول را بر $a_{11}$ تقسیم می‌کند. سپس، الگوریتم سطر اول را به سطرهای باقی‌مانده اضافه می‌کند به طوری که ضرایب در ستون اول همگی صفر شوند. برای رسیدن به این هدف، در سطر i-ام، باید سطر اول را که در $-a_{i1}$ ضرب شده است اضافه کنیم. توجه داشته باشید که این عملیات باید روی بردار $b$ نیز انجام شود. به نوعی، این فرآیند طوری عمل می‌کند که گویی بردار $b$ ستون $(m+1)$-ام ماتریس $A$ است.

در نتیجه، پس از گام اول، ستون اول ماتریس $A$ شامل ۱ در سطر اول و ۰ در سطرهای دیگر خواهد بود.

به طور مشابه، گام دوم الگوریتم را اجرا می‌کنیم، که در آن ستون دوم از سطر دوم را در نظر می‌گیریم. ابتدا، سطر بر $a_{22}$ تقسیم می‌شود، سپس از سطرهای دیگر کم می‌شود تا تمام ستون دوم (به جز سطر دوم) صفر شود.

این فرآیند را برای تمام ستون‌های ماتریس $A$ ادامه می‌دهیم. اگر $n = m$ باشد، آنگاه $A$ به ماتریس همانی تبدیل خواهد شد.

جستجو برای عنصر محوری (pivot)

طرح توصیف‌شده جزئیات زیادی را نادیده گرفته است. در گام i-ام، اگر $a_{ii}$ صفر باشد، نمی‌توانیم روش توصیف‌شده را مستقیماً به کار ببریم. به جای آن، ابتدا باید یک سطر محوری انتخاب کنیم: یعنی سطری از ماتریس را پیدا کنیم که در آن ستون i-ام غیرصفر باشد، و سپس آن دو سطر را با هم جابجا کنیم.

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

در بسیاری از پیاده‌سازی‌ها، حتی وقتی $a_{ii} \neq 0$ است، مشاهده می‌شود که افراد همچنان سطر i-ام را با یک سطر محوری دیگر جابجا می‌کنند، با استفاده از هیوریستیک‌هایی مانند انتخاب سطر محوری با بیشترین قدر مطلق $a_{ji}$. این هیوریستیک برای کاهش دامنه مقادیر ماتریس در گام‌های بعدی استفاده می‌شود. بدون این هیوریستیک، حتی برای ماتریس‌هایی با اندازه حدود ۲۰، خطا بسیار بزرگ شده و می‌تواند باعث سرریز (overflow) برای انواع داده ممیز شناور در C++ شود.

حالت‌های وا رفته (degenerate)

در حالتی که $m=n$ و دستگاه غیر وا رفته باشد (یعنی دترمینان غیرصفر داشته باشد و دارای جواب یکتا باشد)، الگوریتم توصیف‌شده ماتریس $A$ را به ماتریس همانی تبدیل خواهد کرد.

اکنون حالت عمومی را در نظر می‌گیریم، که در آن $n$ و $m$ لزوماً برابر نیستند و دستگاه می‌تواند وا رفته باشد. در این موارد، ممکن است عنصر محوری در گام i-ام پیدا نشود. این به این معنی است که در ستون i-ام، از سطر فعلی به بعد، همگی صفر هستند. در این حالت، یا هیچ مقدار ممکنی برای متغیر $x_i$ وجود ندارد (یعنی دستگاه معادلات خطی جوابی ندارد)، یا $x_i$ یک متغیر مستقل است و می‌تواند هر مقدار دلخواهی بگیرد. هنگام پیاده‌سازی گاوس-جردن، باید کار را برای متغیرهای بعدی ادامه دهید و فقط از ستون i-ام صرف نظر کنید (این معادل حذف ستون i-ام از ماتریس است).

بنابراین، برخی از متغیرها در این فرآیند ممکن است مستقل تشخیص داده شوند. وقتی تعداد متغیرها، $m$، از تعداد معادلات، $n$، بیشتر باشد، حداقل $m - n$ متغیر مستقل پیدا خواهد شد.

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

پیاده‌سازی

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

ورودی تابع gauss ماتریس دستگاه $a$ است. ستون آخر این ماتریس بردار $b$ است.

این تابع تعداد جواب‌های دستگاه را برمی‌گرداند $(0, 1, \text{یا } \infty)$. اگر حداقل یک جواب وجود داشته باشد، در بردار $ans$ برگردانده می‌شود.

const double EPS = 1e-9;
const int INF = 2; // it doesn't actually have to be infinity or a big number

int gauss (vector < vector<double> > a, vector<double> & ans) {
    int n = (int) a.size();
    int m = (int) a[0].size() - 1;

    vector<int> where (m, -1);
    for (int col=0, row=0; col<m && row<n; ++col) {
        int sel = row;
        for (int i=row; i<n; ++i)
            if (abs (a[i][col]) > abs (a[sel][col]))
                sel = i;
        if (abs (a[sel][col]) < EPS)
            continue;
        for (int i=col; i<=m; ++i)
            swap (a[sel][i], a[row][i]);
        where[col] = row;

        for (int i=0; i<n; ++i)
            if (i != row) {
                double c = a[i][col] / a[row][col];
                for (int j=col; j<=m; ++j)
                    a[i][j] -= a[row][j] * c;
            }
        ++row;
    }

    ans.assign (m, 0);
    for (int i=0; i<m; ++i)
        if (where[i] != -1)
            ans[i] = a[where[i]][m] / a[where[i]][i];
    for (int i=0; i<n; ++i) {
        double sum = 0;
        for (int j=0; j<m; ++j)
            sum += ans[j] * a[i][j];
        if (abs (sum - a[i][m]) > EPS)
            return 0;
    }

    for (int i=0; i<m; ++i)
        if (where[i] == -1)
            return INF;
    return 1;
}

نکات پیاده‌سازی:

  • این تابع از دو اشاره‌گر استفاده می‌کند - ستون فعلی $col$ و سطر فعلی $row$.
  • برای هر متغیر $x_i$، مقدار $where(i)$ سطری است که این ستون در آن صفر نیست. این بردار به این دلیل لازم است که برخی متغیرها می‌توانند مستقل باشند.
  • در این پیاده‌سازی، سطر i-ام فعلی بر $a_{ii}$ تقسیم نمی‌شود، همانطور که در بالا توضیح داده شد، بنابراین در پایان ماتریس به ماتریس همانی تبدیل نمی‌شود (هرچند ظاهراً تقسیم سطر i-ام می‌تواند به کاهش خطاها کمک کند).
  • پس از یافتن یک جواب، آن را در ماتریس جایگذاری می‌کنیم تا بررسی کنیم که آیا دستگاه حداقل یک جواب دارد یا نه. اگر جواب آزمایشی موفقیت‌آمیز بود، تابع بسته به اینکه آیا حداقل یک متغیر مستقل وجود دارد یا نه، ۱ یا $\inf$ را برمی‌گرداند.

پیچیدگی

اکنون باید پیچیدگی این الگوریتم را تخمین بزنیم. الگوریتم از $m$ فاز تشکیل شده است که در هر فاز:

  • جستجو و جابجایی سطر محوری. این کار با استفاده از هیوریستیک ذکر شده، $O(n + m)$ زمان می‌برد.
  • اگر عنصر محوری در ستون فعلی پیدا شود - آنگاه باید این معادله را به تمام معادلات دیگر اضافه کنیم که $O(nm)$ زمان می‌برد.

بنابراین، پیچیدگی نهایی الگوریتم $O(\min (n, m) . nm)$ است. در حالتی که $n=m$ باشد، پیچیدگی به سادگی $O(n^3)$ خواهد بود.

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

تسریع الگوریتم

پیاده‌سازی قبلی را می‌توان با تقسیم الگوریتم به دو فاز پیشرو و پسرو، تا دو برابر سریع‌تر کرد:

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

فاز پسرو تنها $O(nm)$ زمان می‌برد که بسیار سریع‌تر از فاز پیشرو است. در فاز پیشرو، تعداد عملیات را به نصف کاهش می‌دهیم، و در نتیجه زمان اجرای پیاده‌سازی را کم می‌کنیم.

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

برای حل دستگاه معادلات خطی در یک پیمانه خاص، همچنان می‌توانیم از الگوریتم توصیف‌شده استفاده کنیم. با این حال، در صورتی که پیمانه برابر با دو باشد، می‌توانیم حذف گاوس-جردن را با استفاده از عملیات بیتی و نوع داده bitset در C++ بسیار مؤثرتر انجام دهیم:

int gauss (vector < bitset<N> > a, int n, int m, bitset<N> & ans) {
    vector<int> where (m, -1);
    for (int col=0, row=0; col<m && row<n; ++col) {
        for (int i=row; i<n; ++i)
            if (a[i][col]) {
                swap (a[i], a[row]);
                break;
            }
        if (! a[row][col])
            continue;
        where[col] = row;

        for (int i=0; i<n; ++i)
            if (i != row && a[i][col])
                a[i] ^= a[row];
        ++row;
    }
        // باقی پیاده‌سازی مانند بالا است
}

از آنجایی که از فشرده‌سازی بیتی استفاده می‌کنیم، پیاده‌سازی نه تنها کوتاه‌تر، بلکه ۳۲ برابر سریع‌تر است.

نکته‌ای کوچک در مورد هیوریستیک‌های مختلف انتخاب سطر محوری

هیچ قاعده کلی برای اینکه از کدام هیوریستیک استفاده کنیم وجود ندارد.

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

با این حال، باید توجه داشت که هر دو هیوریستیک به میزان مقیاس‌دهی معادلات اصلی بستگی دارند. به عنوان مثال، اگر یکی از معادلات در $10^6$ ضرب شده باشد، آن معادله تقریباً به طور قطع در گام اول به عنوان محوری انتخاب خواهد شد. این موضوع کمی عجیب به نظر می‌رسد، بنابراین منطقی است که به یک هیوریستیک پیچیده‌تر به نام implicit pivoting (محورگیری ضمنی) روی بیاوریم.

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

بهبود جواب

با وجود هیوریستیک‌های مختلف، الگوریتم گاوس-جردن همچنان می‌تواند در ماتریس‌های خاص حتی با اندازه ۵۰ تا ۱۰۰ به خطاهای بزرگی منجر شود.

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

بنابراین، راه‌حل به یک فرآیند دو مرحله‌ای تبدیل می‌شود: ابتدا، الگوریتم گاوس-جردن اعمال می‌شود و سپس یک روش عددی که جواب مرحله اول را به عنوان جواب اولیه می‌گیرد، به کار می‌رود.

مسائل تمرینی