Шугаман тэгшитгэлийн системийг бодох Гауссын арга¶

$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$ нь 1-ээс $n$, $j$ нь 1-ээс $m$ хүртэл) ба $b_i$ ($i$ нь 1-ээс $n$ хүртэл) мэдэгдэж байгаа бөгөөд $x_i$ хувьсагчид ($i$ нь 1-ээс $m$ хүртэл) нь үл мэдэгдэгчид юм.

Энэ бодлого мөн энгийн матрицан илэрхийлэлтэй:

$$Ax = b,$$

энд $A$ нь $a_{ij}$ коэффициентүүдийн $n \times m$ хэмжээтэй матриц, $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}$$

Гаусс¶

Хатуухан хэлбэл доор тайлбарласан аргыг "Гаусс-Жордан" буюу Гаусс-Жорданы арга гэж нэрлэх ёстой, учир нь энэ нь 1887 онд Жорданы тайлбарласан Гауссын аргын нэг хувилбар юм.

Тойм¶

Алгоритм нь тэгшитгэл бүр дэх хувьсагчдыг тэгшитгэл бүрт зөвхөн нэг хувьсагч үлдэх хүртэл дараалан арилгах явдал юм. Хэрэв $n = m$ бол үүнийг $A$ матрицыг нэгж матриц болгон хувиргаж, шийд нь цор ганц бөгөөд $b_i$ коэффициенттэй тэнцүү байх энэ илэрхий тохиолдолд тэгшитгэлийг бодож байгаа гэж бодож болно.

Гауссын арга нь хоёр энгийн хувиргалт дээр суурилдаг:

Хоёр тэгшитгэлийг солих боломжтой
Дурын тэгшитгэлийг тухайн мөрийн (тэг биш коэффициенттэй) ба зарим бусад мөрийн (дурын коэффициенттэй) шугаман комбинацаар солиж болно.

Эхний алхамд Гаусс-Жорданы алгоритм эхний мөрийг $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++-ийн хөвөгч цэгтэй өгөгдлийн төрлүүдэд халилт үүсгэж болно.

Degenerate тохиолдлууд¶

$m = n$ бөгөөд систем degenerate биш (өөрөөр хэлбэл тэг биш тодорхойлогчтой, цор ганц шийдтэй) тохиолдолд дээр тайлбарласан алгоритм $A$-г нэгж матриц болгон хувиргана.

Одоо бид $n$ ба $m$ заавал тэнцүү биш бөгөөд систем degenerate байж болох ерөнхий тохиолдол-ыг авч үзнэ. Эдгээр тохиолдолд $i$ дугаар алхмын тулгуур элемент олдохгүй байж болно. Энэ нь $i$ дугаар баганад одоогийн мөрөөс эхлэн бүгд тэг агуулж байна гэсэн үг. Энэ тохиолдолд $x_i$ хувьсагчийн боломжит утга байхгүй (SLAE шийдгүй гэсэн үг), эсвэл $x_i$ нь хамааралгүй хувьсагч бөгөөд дурын утга авч болно. Гаусс-Жорданыг хэрэгжүүлэхдээ та дараагийн хувьсагчдын хувьд ажлаа үргэлжлүүлж, $i$ дугаар баганыг зүгээр л алгасах ёстой (энэ нь матрицын $i$ дугаар баганыг устгахтай эквивалент).

Ингэснээр процессын явцад зарим хувьсагчид хамааралгүй болох нь олдож болно. $m$ хувьсагчийн тоо нь $n$ тэгшитгэлийн тооноос их бол дор хаяж $m - n$ хамааралгүй хувьсагч олдоно.

Ерөнхийдөө хэрэв та дор хаяж нэг хамааралгүй хувьсагч олбол энэ нь дурын утга авч болох ба бусад (хамааралтай) хувьсагчид түүгээр дамжуулан илэрхийлэгдэнэ. Энэ нь бид бодит тооны талбарт ажиллаж байх үед систем хязгааргүй олон шийдтэй байж болно гэсэн үг. Гэвч хамааралгүй хувьсагч байх үед SLAE огт шийдгүй байж болохыг санах ёстой. Энэ нь үлдсэн боловсруулагдаагүй тэгшитгэлүүд дор хаяж нэг тэг биш тогтмол гишүүнтэй байх үед тохиолдоно. Үүнийг бүх хамааралгүй хувьсагчид тэг оноож, бусад хувьсагчийг тооцоолоод, дараа нь анхны SLAE-д орлуулж, тэдгээрийг хангаж байгаа эсэхийг шалгаж болно.

Implementation¶

Following is an implementation of Gauss-Jordan. Choosing the pivot row is done with heuristic: choosing maximum value in the current column.

The input to the function gauss is the system matrix $a$. The last column of this matrix is vector $b$.

The function returns the number of solutions of the system $(0, 1,\textrm{or } \infty)$. If at least one solution exists, then it is returned in the vector $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;
}

Implementation notes:

The function uses two pointers - the current column $col$ and the current row $row$.
For each variable $x_i$, the value $where(i)$ is the line where this column is not zero. This vector is needed because some variables can be independent.
In this implementation, the current $i$th line is not divided by $a_{ii}$ as described above, so in the end the matrix is not identity matrix (though apparently dividing the $i$th line can help reducing errors).
After finding a solution, it is inserted back into the matrix - to check whether the system has at least one solution or not. If the test solution is successful, then the function returns 1 or $\inf$, depending on whether there is at least one independent variable.

Complexity¶

Одоо бид энэ алгоритмын complexity-г үнэлэх ёстой. Алгоритм нь $m$ үе шатаас тогтох ба үе шат бүрд:

Тулгуур мөрийг хайж дахин байрлуулна. Дээр дурдсан эвристикийг ашиглах үед энэ нь $O(n + m)$ болно.
Хэрэв одоогийн баганад тулгуур элемент олдвол бид энэ тэгшитгэлийг бусад бүх тэгшитгэлд нэмэх ёстой ба энэ нь $O(nm)$ хугацаа авна.

Тиймээс алгоритмын эцсийн complexity нь $O(\min (n, m) . nm)$ болно. $n = m$ тохиолдолд complexity нь зүгээр л $O(n^3)$ болно.

SLAE нь бодит тоон дээр биш, харин хоёрын модулиар байх үед системийг хамаагүй хурдан бодож болохыг анхаарна уу, үүнийг доор тайлбарлав.

Алгоритмыг хурдасгах¶

Өмнөх хэрэгжүүлэлтийг алгоритмыг шууд ба урвуу гэсэн хоёр үе шатанд хуваах замаар хоёр дахин хурдасгаж болно:

Шууд үе шат: Өмнөх хэрэгжүүлэлттэй төстэй боловч одоогийн мөрийг зөвхөн түүнээс хойших мөрүүдэд нэмнэ. Үр дүнд нь бид диагональ матрицын оронд гурвалжин матриц олж авна.
Урвуу үе шат: Матриц гурвалжин болсон үед бид эхлээд сүүлчийн хувьсагчийн утгыг тооцоолно. Дараа нь энэ утгыг орлуулж дараагийн хувьсагчийн утгыг олно. Дараа нь эдгээр хоёр утгыг орлуулж дараагийн хувьсагчдыг олно...

Урвуу үе шат зөвхөн $O(nm)$ авах ба энэ нь шууд үе шатнаас хамаагүй хурдан. Шууд үе шатанд бид үйлдлийн тоог хоёр дахин багасгаснаар хэрэгжүүлэлтийн ажиллах хугацааг багасгана.

Модуль SLAE-г бодох¶

Ямар нэг модулиар SLAE-г бодоход бид тайлбарласан алгоритмыг ашигласаар болно. Гэвч модуль нь 2-той тэнцүү тохиолдолд бид битийн үйлдэл ба 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
}

Бид бит шахалт ашигладаг тул хэрэгжүүлэлт нь зөвхөн богино төдийгүй 32 дахин хурдан юм.

Тулгуур мөр сонгох өөр өөр эвристикийн тухай бага зэрэг тэмдэглэл¶

Ямар эвристик ашиглах талаар ерөнхий дүрэм байхгүй.

Өмнөх хэрэгжүүлэлтэд ашигласан эвристик практикт нэлээд сайн ажилладаг. Мөн энэ нь "бүрэн тулгах" (тулгуур мөрийг бүхэл дэд матрицын бүх элементийн дундаас (одоогийн мөр ба одоогийн баганаас) хайдаг)-тай бараг ижил хариу өгдөг нь тогтоогддог.

Гэвч хоёр эвристик хоёулаа анхны тэгшитгэлүүд хэр их масштаблагдсанаас хамаардгийг анхаарах хэрэгтэй. Жишээ нь хэрэв тэгшитгэлийн нэг нь $10^6$-аар үржүүлэгдсэн бол энэ тэгшитгэл эхний алхамд тулгуур болгон сонгогдох нь бараг тодорхой. Энэ нь нэлээд хачирхалтай санагдах тул далд тулгах хэмээх илүү төвөгтэй эвристик рүү шилжих нь логиктой мэт.

Далд тулгах нь хамгийн их элемент нэгж болохоор хоёр мөрийг хоёуланг нь нормчилсон мэт элементүүдийг харьцуулна. Энэ аргыг хэрэгжүүлэхийн тулд мөр бүрд хамгийн ихийг хадгалах хэрэгтэй (эсвэл мөр бүрийг хамгийн их нь нэгж болгон хадгалах ч энэ нь хуримтлагдсан алдааг нэмэгдүүлж болзошгүй).

Шийдийг сайжруулах¶

Янз бүрийн эвристикийг үл харгалзан Гаусс-Жорданы алгоритм $50 - 100$ хэмжээтэй тусгай матрицад ч том алдаа гаргаж болно.

Тиймээс гарсан Гаусс-Жорданы шийдийг заримдаа энгийн тоон арга — жишээ нь энгийн итерацийн арга хэрэглэн сайжруулах ёстой.

Ингэснээр шийд хоёр алхамтай болно: Эхлээд Гаусс-Жорданы алгоритмыг хэрэглэх ба дараа нь эхний алхмын шийдийг анхны шийд болгон авах тоон аргыг хэрэглэнэ.