লিনিয়ার ইকুয়েশন সিস্টেম সমাধানের গাউস মেথড#
$m$ টা অজানা রাশি বিশিষ্ট $n$ টা লিনিয়ার বীজগাণিতিক ইকুয়েশনের একটা সিস্টেম (SLAE) দেওয়া আছে। তোমাকে সিস্টেমটা সমাধান করতে হবে: এর কোনো সমাধান নেই, ঠিক একটা সমাধান আছে, নাকি অসীম সংখ্যক সমাধান আছে তা নির্ণয় করতে হবে। এবং যদি কমপক্ষে একটা সমাধান থাকে, তাহলে যেকোনো একটা বের করতে হবে।
ফর্মালি, সমস্যাটা এরকম: সিস্টেমটা সমাধান করো:
$$\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$ এর সমান।
গাউসীয় এলিমিনেশন দুটি সরল রূপান্তরের উপর ভিত্তি করে:
- দুটি ইকুয়েশন পরস্পর বিনিময় করা সম্ভব
- যেকোনো ইকুয়েশনকে সেই রোর (row) (অশূন্য সহগ সহ) এবং অন্য কিছু রোর (যথেচ্ছ সহগ সহ) একটি লিনিয়ার কম্বিনেশন দ্বারা প্রতিস্থাপন করা যায়।
প্রথম ধাপে, গাউস-জর্ডান অ্যালগরিদম প্রথম রোকে $a_{11}$ দ্বারা ভাগ করে। তারপর, অ্যালগরিদম প্রথম রোটি অবশিষ্ট রোগুলোতে এমনভাবে যোগ করে যাতে প্রথম কলামের সকল সহগ শূন্য হয়ে যায়। এটি অর্জন করতে, i-তম রোতে, আমাদের $- a_{i1}$ দ্বারা গুণিত প্রথম রো যোগ করতে হবে। লক্ষ্য করুন যে, এই অপারেশনটি ভেক্টর $b$-তেও করতে হবে। একটি অর্থে, এটি এমনভাবে আচরণ করে যেন ভেক্টর $b$ হলো ম্যাট্রিক্স $A$ এর $m+1$-তম কলাম।
তাহলে, প্রথম ধাপের পরে, ম্যাট্রিক্স $A$ এর প্রথম কলামে প্রথম রোতে $1$ এবং অন্যান্য রোতে $0$ থাকবে।
একইভাবে, আমরা অ্যালগরিদমের দ্বিতীয় ধাপ সম্পাদন করি, যেখানে আমরা দ্বিতীয় রোর দ্বিতীয় কলাম কনসিডার করি। প্রথমে, রোটিকে $a_{22}$ দ্বারা ভাগ করা হয়, তারপর অন্যান্য রো থেকে এটি বিয়োগ করা হয় যাতে সমস্ত দ্বিতীয় কলাম $0$ হয়ে যায় (দ্বিতীয় রো ব্যতীত)।
আমরা ম্যাট্রিক্স $A$ এর সকল কলামের জন্য এই প্রক্রিয়া চালিয়ে যাই। যদি $n = m$ হয়, তাহলে $A$ আইডেন্টিটি ম্যাট্রিক্সে পরিণত হবে।
পিভটিং এলিমেন্ট সার্চ#
উপরের স্কিমটা অনেক বিশদ বিবরণ বাদ দিয়েছে। $i$-তম ধাপে, যদি $a_{ii}$ শূন্য হয়, তাহলে আমরা সরাসরি উপরের মেথড প্রয়োগ করতে পারি না। পরিবর্তে, আমাদের প্রথমে একটা পিভটিং রো নির্বাচন করতে হবে: ম্যাট্রিক্সের এমন একটা রো খুঁজে বের করতে হবে যেখানে $i$-তম কলাম অশূন্য, এবং তারপর দুটো রো অদলবদল করতে হবে।
লক্ষ্য করো যে, এখানে আমরা রো অদলবদল করি কিন্তু কলাম নয়। কারণ তুমি যদি কলাম অদলবদল করো, তাহলে সমাধান পাওয়ার পর তোমাকে সঠিক স্থানে ফিরিয়ে আনার কথা মনে রাখতে হবে। তাই রো অদলবদল করা অনেক সহজ।
অনেক ইমপ্লিমেন্টেশনে, যখন $a_{ii} \neq 0$, তখনও $i$-তম রোকে কোনো পিভটিং রোর সাথে অদলবদল করা হয়, কিছু হিউরিস্টিক ব্যবহার করে যেমন $a_{ji}$ এর সর্বোচ্চ পরম মান বিশিষ্ট পিভটিং রো বেছে নেওয়া। এই হিউরিস্টিকটা পরবর্তী ধাপগুলোতে ম্যাট্রিক্সের মানের পরিসর কমাতে ব্যবহৃত হয়। এই হিউরিস্টিক ছাড়া, প্রায় $20$ আকারের ম্যাট্রিক্সের জন্যও ত্রুটি অনেক বেশি হবে এবং C++ এর ফ্লোটিং পয়েন্ট ডেটা টাইপে ওভারফ্লো হতে পারে।
ডিজেনারেট কেস#
যে ক্ষেত্রে $m = n$ এবং সিস্টেমটি নন-ডিজেনারেট (অর্থাৎ এর অশূন্য ডিটারমিন্যান্ট আছে, এবং অনন্য সমাধান আছে), উপরে বর্ণিত অ্যালগরিদম $A$ কে আইডেন্টিটি ম্যাট্রিক্সে রূপান্তর করবে।
এখন আমরা সাধারণ ক্ষেত্র কনসিডার করি, যেখানে $n$ এবং $m$ অগত্যা সমান নয়, এবং সিস্টেমটা ডিজেনারেট হতে পারে। এই ক্ষেত্রগুলোতে, $i$-তম ধাপে পিভটিং এলিমেন্ট পাওয়া নাও যেতে পারে। এর মানে হলো $i$-তম কলামে, বর্তমান রো থেকে শুরু করে, সবগুলোতে শূন্য আছে। এই ক্ষেত্রে, হয় ভ্যারিয়েবল $x_i$ এর কোনো সম্ভাব্য মান নেই (মানে SLAE এর কোনো সমাধান নেই), অথবা $x_i$ একটা ফ্রি ভ্যারিয়েবল (free variable) এবং যথেচ্ছ মান নিতে পারে। গাউস-জর্ডান ইমপ্লিমেন্ট করার সময়, তোমার পরবর্তী ভ্যারিয়েবলগুলোর জন্য কাজ চালিয়ে যাওয়া উচিত এবং $i$-তম কলাম এড়িয়ে যাওয়া উচিত (এটা ম্যাট্রিক্সের $i$-তম কলাম অপসারণের সমতুল্য)।
তাহলে, প্রক্রিয়ায় কিছু ভ্যারিয়েবল ফ্রি হিসেবে পাওয়া যেতে পারে। যখন ভ্যারিয়েবলের সংখ্যা $m$, ইকুয়েশনের সংখ্যা $n$ এর চেয়ে বেশি, তখন কমপক্ষে $m - n$ টা ফ্রি ভ্যারিয়েবল পাওয়া যাবে।
সাধারণত, যদি তুমি কমপক্ষে একটা ফ্রি ভ্যারিয়েবল পাও, এটা যেকোনো যথেচ্ছ মান নিতে পারে, যেখানে অন্যান্য (ডিপেন্ডেন্ট) ভ্যারিয়েবলগুলো এর মাধ্যমে প্রকাশিত হয়। এর মানে হলো যখন আমরা রিয়েল নাম্বারের ক্ষেত্রে কাজ করি, সিস্টেমের সম্ভাব্যভাবে অসীম সংখ্যক সমাধান থাকতে পারে। কিন্তু মনে রেখো যে ফ্রি ভ্যারিয়েবল থাকলেও SLAE এর কোনো সমাধান না-ও থাকতে পারে। এটা ঘটে যখন অবশিষ্ট অপ্রসেসড ইকুয়েশনে কমপক্ষে একটা অশূন্য ধ্রুবক পদ থাকে। তুমি সকল ফ্রি ভ্যারিয়েবলে শূন্য বসিয়ে, অন্যান্য ভ্যারিয়েবল বের করে, এবং তারপর মূল SLAE-তে বসিয়ে যাচাই করতে পারো।
ইমপ্লিমেন্টেশন#
নিচে গাউস-জর্ডানের একটি ইমপ্লিমেন্টেশন দেওয়া হলো। পিভট রো নির্বাচন হিউরিস্টিক দিয়ে করা হয়: বর্তমান কলামে সর্বোচ্চ মান নির্বাচন।
ফাংশন gauss এর ইনপুট হলো সিস্টেমের ম্যাট্রিক্স $a$। এই ম্যাট্রিক্সের শেষ কলামটি হলো ভেক্টর $b$।
ফাংশনটা সিস্টেমের সমাধানের সংখ্যা $(0, 1,\textrm{or } \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)$।
লক্ষ্য করো যে যখন SLAE রিয়েল নাম্বারে নয়, বরং মডুলো ২-তে হয়, তখন সিস্টেমটা অনেক দ্রুত সমাধান করা যায়, যেটা নিচে বলা হয়েছে।
অ্যালগরিদমের ত্বরণ#
আগের ইমপ্লিমেন্টেশনকে দুই গুণ দ্রুত করা যায়, অ্যালগরিদমকে দুটি ফেজে ভাগ করে: ফরওয়ার্ড এবং রিভার্স:
- ফরওয়ার্ড ফেজ: আগের ইমপ্লিমেন্টেশনের অনুরূপ, কিন্তু বর্তমান রো শুধুমাত্র এর পরের রোগুলোতে যোগ করা হয়। তাহলে, আমরা ডায়াগোনালের পরিবর্তে একটা ট্রায়াঙ্গুলার ম্যাট্রিক্স পাই।
- রিভার্স ফেজ: যখন ম্যাট্রিক্স ত্রিভুজাকার, আমরা প্রথমে শেষ চলকের মান গণনা করি। তারপর পরবর্তী চলকের মান বের করতে এই মান বসাই। তারপর এই দুটি মান বসিয়ে পরবর্তী চলক বের করি…
রিভার্স ফেজে শুধু $O(nm)$ সময় লাগে, যা ফরওয়ার্ড ফেজের চেয়ে অনেক দ্রুত। ফরওয়ার্ড ফেজে, আমরা অপারেশনের সংখ্যা অর্ধেক করি, ফলে ইমপ্লিমেন্টেশনের রানিং টাইম কমে।
মডুলার SLAE সমাধান#
কোনো মডুলাসে SLAE সমাধানের জন্য, আমরা এখনও বর্ণিত অ্যালগরিদম ব্যবহার করতে পারি। তবে, যদি মডুলাস ২ এর সমান হয়, আমরা বিটওয়াইজ অপারেশন এবং C++ bitset ডেটা টাইপ ব্যবহার করে গাউস-জর্ডান এলিমিনেশন অনেক বেশি কার্যকরভাবে করতে পারি:
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;
}
// The rest of implementation is the same as above
}যেহেতু আমরা বিট কম্প্রেস ব্যবহার করি, ইমপ্লিমেন্টেশন শুধু ছোটই নয়, ৩২ গুণ দ্রুতও।
পিভটিং রো নির্বাচনের বিভিন্ন হিউরিস্টিক সম্পর্কে একটি ছোট নোট#
কোন হিউরিস্টিক ব্যবহার করতে হবে তার কোনো সাধারণ নিয়ম নেই।
আগের ইমপ্লিমেন্টেশনে ব্যবহৃত হিউরিস্টিকটি বাস্তবে বেশ ভালো কাজ করে। এটি “ফুল পিভটিং”-এর প্রায় একই উত্তর দেয় - যেখানে সম্পূর্ণ সাবম্যাট্রিক্সের (বর্তমান রো এবং বর্তমান কলাম থেকে) সকল উপাদানের মধ্যে পিভটিং রো সার্চ করা হয়।
তবে, তোমার লক্ষ্য করা উচিত যে উভয় হিউরিস্টিকই মূল ইকুয়েশনগুলো কতটুকু স্কেল করা হয়েছে তার উপর নির্ভরশীল। উদাহরণস্বরূপ, যদি একটা ইকুয়েশনকে $10^6$ দ্বারা গুণ করা হয়, তাহলে এই ইকুয়েশনটা প্রথম ধাপে পিভট হিসেবে সিলেক্টেড হওয়া প্রায় নিশ্চিত। এটা বেশ অদ্ভুত মনে হয়, তাই ইমপ্লিসিট পিভটিং নামক আরও জটিল হিউরিস্টিকে পরিবর্তন করা যুক্তিসঙ্গত মনে হয়।
ইমপ্লিসিট পিভটিং উপাদানগুলোকে এমনভাবে তুলনা করে যেন উভয় লাইন নরমালাইজড, যাতে সর্বোচ্চ উপাদান একক হয়। এই কৌশল ইমপ্লিমেন্ট করতে, প্রতিটি রোতে সর্বোচ্চ বজায় রাখতে হবে (অথবা প্রতিটি লাইন এমনভাবে বজায় রাখতে হবে যাতে সর্বোচ্চ একক হয়, কিন্তু এটি সঞ্চিত ত্রুটি বাড়াতে পারে)।
সমাধানের উন্নতি#
বিভিন্ন হিউরিস্টিক থাকা সত্ত্বেও, গাউস-জর্ডান অ্যালগরিদম $50 - 100$ আকারের বিশেষ ম্যাট্রিক্সেও বড় ত্রুটি তৈরি করতে পারে।
অতএব, গাউস-জর্ডান থেকে প্রাপ্ত সমাধানকে কখনো কখনো একটি সরল সংখ্যাসূচক মেথড প্রয়োগ করে উন্নত করতে হয় - উদাহরণস্বরূপ, সিম্পল ইটারেশন মেথড।
এইভাবে, সমাধানটি দুই-ধাপের হয়ে যায়: প্রথমে, গাউস-জর্ডান অ্যালগরিদম প্রয়োগ করা হয়, এবং তারপর একটি সংখ্যাসূচক মেথড প্রথম ধাপের সমাধানকে প্রাথমিক সমাধান হিসেবে নেয়।