Метод Гаусса розв'язання системи лінійних рівнянь

Дано систему з $n$ лінійних алгебраїчних рівнянь (СЛАР) із $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$ від 1 до $n$, $j$ від 1 до $m$) і $b_i$ ($i$ від 1 до $n$) відомі, а змінні $x_i$ ($i$ від 1 до $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}$$

Коли підходить цей алгоритм?

• Потрібно розв'язати систему лінійних рівнянь $Ax=b$ (над дійсними числами або за модулем) і з'ясувати кількість розв'язків?
• Цікавить лише кількість лінійно незалежних рівнянь / стовпців, а не сам розв'язок? (якщо так → Ранг матриці)
• Потрібен лише визначник квадратної матриці, без розв'язку? (якщо так → Визначник методом Гаусса)

Гаусс

Строго кажучи, описаний нижче метод варто було б назвати «метод Гаусса-Жордана», або виключенням за Гауссом-Жорданом, оскільки це варіація методу Гаусса, описана Жорданом 1887 року.

Огляд

Алгоритм полягає в послідовному виключенні змінних у кожному рівнянні доти, доки в кожному рівнянні не залишиться лише одна змінна. Якщо $n = m$, можна уявляти це як перетворення матриці $A$ на одиничну матрицю та розв'язання рівняння в цьому очевидному випадку, де розв'язок єдиний і дорівнює коефіцієнту $b_i$.

Метод Гаусса ґрунтується на двох простих перетвореннях:

• Можна поміняти місцями два рівняння.
• Будь-яке рівняння можна замінити лінійною комбінацією цього рядка (з ненульовим коефіцієнтом) і деяких інших рядків (з довільними коефіцієнтами).

На першому кроці алгоритм Гаусса-Жордана ділить перший рядок на $a_{11}$. Потім алгоритм додає перший рядок до решти рядків так, щоб коефіцієнти в першому стовпці стали нульовими. Щоб цього досягти, до $i$-го рядка ми маємо додати перший рядок, помножений на $- a_{i1}$. Зауважимо, що цю операцію треба також виконувати над вектором $b$. У певному сенсі це поводиться так, ніби вектор $b$ був $m+1$-м стовпцем матриці $A$.

У результаті після першого кроку перший стовпець матриці $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$ (тобто СЛАР не має розв'язку), або $x_i$ — незалежна змінна й може набувати довільного значення. Реалізуючи метод Гаусса-Жордана, слід продовжити роботу для наступних змінних і просто пропустити $i$-й стовпець (це еквівалентно видаленню $i$-го стовпця матриці).

Отже, у процесі можна виявити, що деякі змінні є незалежними. Коли кількість змінних $m$ більша за кількість рівнянь $n$, то знайдеться щонайменше $m - n$ незалежних змінних.

Загалом, якщо ви знайшли хоча б одну незалежну змінну, вона може набувати будь-якого довільного значення, а решта (залежних) змінних виражаються через неї. Це означає, що, працюючи в полі дійсних чисел, система потенційно має нескінченно багато розв'язків. Але слід пам'ятати, що за наявності незалежних змінних СЛАР може взагалі не мати розв'язку. Це трапляється, коли в решти необроблених рівнянь є хоча б один ненульовий вільний член. Це можна перевірити, присвоївши всім незалежним змінним нулі, обчисливши інші змінні, а потім підставивши їх у початкову СЛАР, щоб перевірити, чи вони її задовольняють.

Реалізація

Нижче наведено реалізацію методу Гаусса-Жордана. Вибір опорного рядка виконується за евристикою: обираємо максимальне значення в поточному стовпці.

Вхідними даними функції gauss є матриця системи $a$. Останній стовпець цієї матриці — це вектор $b$.

Функція повертає кількість розв'язків системи $(0, 1,\textrm{або } \infty)$. Якщо існує хоча б один розв'язок, то він повертається у векторі $ans$.

C++

const double EPS = 1e-9;
const int INF = 2; // насправді це не обов'язково має бути нескінченність чи велике число

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;
}

Python

EPS = 1e-9
INF = 2  # насправді це не обов'язково має бути нескінченність чи велике число


def gauss(a):
    # a — матриця системи; останній стовпець кожного рядка — вектор b.
    # Функція повертає кортеж (кількість розв'язків, вектор ans).
    a = [row[:] for row in a]  # копія, щоб не змінювати вхідні дані
    n = len(a)
    m = len(a[0]) - 1

    where = [-1] * m
    row = 0
    for col in range(m):
        if row >= n:
            break
        sel = row
        for i in range(row, n):
            if abs(a[i][col]) > abs(a[sel][col]):
                sel = i
        if abs(a[sel][col]) < EPS:
            continue
        a[sel], a[row] = a[row], a[sel]
        where[col] = row

        for i in range(n):
            if i != row:
                c = a[i][col] / a[row][col]
                for j in range(col, m + 1):
                    a[i][j] -= a[row][j] * c
        row += 1

    ans = [0.0] * m
    for i in range(m):
        if where[i] != -1:
            ans[i] = a[where[i]][m] / a[where[i]][i]
    for i in range(n):
        s = sum(ans[j] * a[i][j] for j in range(m))
        if abs(s - a[i][m]) > EPS:
            return 0, ans

    for i in range(m):
        if where[i] == -1:
            return INF, ans
    return 1, ans

TypeScript

const EPS = 1e-9;
const INF = 2; // насправді це не обов'язково має бути нескінченність чи велике число

// Повертає [кількість розв'язків, вектор ans].
// a — матриця системи; останній стовпець кожного рядка — вектор b.
function gauss(a: number[][]): [number, number[]] {
    a = a.map((row) => row.slice()); // копія, щоб не змінювати вхідні дані
    const n = a.length;
    const m = a[0].length - 1;

    const where: number[] = new Array(m).fill(-1);
    for (let col = 0, row = 0; col < m && row < n; ++col) {
        let sel = row;
        for (let i = row; i < n; ++i)
            if (Math.abs(a[i][col]) > Math.abs(a[sel][col])) sel = i;
        if (Math.abs(a[sel][col]) < EPS) continue;
        for (let i = col; i <= m; ++i)
            [a[sel][i], a[row][i]] = [a[row][i], a[sel][i]];
        where[col] = row;

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

    const ans: number[] = new Array(m).fill(0);
    for (let i = 0; i < m; ++i)
        if (where[i] !== -1) ans[i] = a[where[i]][m] / a[where[i]][i];
    for (let i = 0; i < n; ++i) {
        let s = 0;
        for (let j = 0; j < m; ++j) s += ans[j] * a[i][j];
        if (Math.abs(s - a[i][m]) > EPS) return [0, ans];
    }

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

Go

const eps = 1e-9
const inf = 2 // насправді це не обов'язково має бути нескінченність чи велике число

// gauss повертає кількість розв'язків і вектор ans.
// a — матриця системи; останній стовпець кожного рядка — вектор b.
func gauss(a [][]float64) (int, []float64) {
	n := len(a)
	m := len(a[0]) - 1

	where := make([]int, m)
	for i := range where {
		where[i] = -1
	}
	for col, row := 0, 0; col < m && row < n; col++ {
		sel := row
		for i := row; i < n; i++ {
			if math.Abs(a[i][col]) > math.Abs(a[sel][col]) {
				sel = i
			}
		}
		if math.Abs(a[sel][col]) < eps {
			continue
		}
		for i := col; i <= m; i++ {
			a[sel][i], a[row][i] = a[row][i], a[sel][i]
		}
		where[col] = row

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

	ans := make([]float64, m)
	for i := 0; i < m; i++ {
		if where[i] != -1 {
			ans[i] = a[where[i]][m] / a[where[i]][i]
		}
	}
	for i := 0; i < n; i++ {
		sum := 0.0
		for j := 0; j < m; j++ {
			sum += ans[j] * a[i][j]
		}
		if math.Abs(sum-a[i][m]) > eps {
			return 0, ans
		}
	}

	for i := 0; i < m; i++ {
		if where[i] == -1 {
			return inf, ans
		}
	}
	return 1, ans
}

Зауваження щодо реалізації:

• Функція використовує два вказівники — поточний стовпець $col$ і поточний рядок $row$.
• Для кожної змінної $x_i$ значення $where(i)$ — це рядок, у якому цей стовпець ненульовий. Цей вектор потрібен, бо деякі змінні можуть бути незалежними.
• У цій реалізації поточний $i$-й рядок не ділиться на $a_{ii}$, як описано вище, тож наприкінці матриця не є одиничною (хоча, вочевидь, ділення $i$-го рядка може допомогти зменшити похибки).
• Знайшовши розв'язок, його підставляють назад у матрицю — щоб перевірити, чи має система хоча б один розв'язок. Якщо пробний розв'язок успішний, то функція повертає 1 або $\inf$ залежно від того, чи є хоча б одна незалежна змінна.

Складність

Тепер оцінимо складність цього алгоритму. Алгоритм складається з $m$ фаз, у кожній фазі:

• Пошук і перестановка опорного рядка. Це займає $O(n + m)$ при використанні згаданої вище евристики.
• Якщо опорний елемент у поточному стовпці знайдено — то ми маємо додати це рівняння до всіх інших рівнянь, що займає час $O(nm)$.

Отже, остаточна складність алгоритму становить $O(\min (n, m) . nm)$. У випадку $n = m$ складність дорівнює просто $O(n^3)$.

Зауважимо, що коли СЛАР задана не над дійсними числами, а за модулем два, то систему можна розв'язати значно швидше, про що йдеться нижче.

Пришвидшення алгоритму

Попередню реалізацію можна пришвидшити вдвічі, поділивши алгоритм на дві фази: пряму та зворотну:

• Пряма фаза: подібна до попередньої реалізації, але поточний рядок додається лише до рядків після нього. У результаті ми отримуємо трикутну матрицю замість діагональної.
• Зворотна фаза: коли матриця трикутна, ми спочатку обчислюємо значення останньої змінної. Потім підставляємо це значення, щоб знайти значення наступної змінної. Потім підставляємо ці два значення, щоб знайти наступні змінні...

Зворотна фаза займає лише $O(nm)$, що набагато швидше за пряму фазу. У прямій фазі ми зменшуємо кількість операцій удвічі, тим самим скорочуючи час роботи реалізації.

Розв'язання модульної СЛАР

Для розв'язання СЛАР за деяким модулем ми все одно можемо використати описаний алгоритм. Однак у випадку, коли модуль дорівнює 2, ми можемо виконати виключення за Гауссом-Жорданом значно ефективніше, використовуючи бітові операції та тип даних bitset у C++:

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;
	}
        // Решта реалізації така сама, як вище
}

Python

def gauss_mod(a, n, m):
    # Рядки подано як цілі-бітові маски: біти 0..m-1 — коефіцієнти, біт m — вектор b.
    # XOR двох масок (a[i] ^= a[row]) замінює побітову операцію над bitset.
    a = a[:]
    where = [-1] * m
    row = 0
    for col in range(m):
        if row >= n:
            break
        for i in range(row, n):
            if (a[i] >> col) & 1:
                a[i], a[row] = a[row], a[i]
                break
        if not ((a[row] >> col) & 1):
            continue
        where[col] = row

        for i in range(n):
            if i != row and ((a[i] >> col) & 1):
                a[i] ^= a[row]
        row += 1
    # Решта реалізації така сама, як вище

TypeScript

// Рядки — бітові маски (bigint): біти 0..m-1 — коефіцієнти, біт m — вектор b.
// XOR двох масок (a[i] ^= a[row]) замінює побітову операцію над bitset.
function gaussMod(a: bigint[], n: number, m: number): [number, bigint] {
    a = a.slice();
    const where: number[] = new Array(m).fill(-1);
    for (let col = 0, row = 0; col < m && row < n; ++col) {
        for (let i = row; i < n; ++i)
            if ((a[i] >> BigInt(col)) & 1n) {
                [a[i], a[row]] = [a[row], a[i]];
                break;
            }
        if (!((a[row] >> BigInt(col)) & 1n)) continue;
        where[col] = row;

        for (let i = 0; i < n; ++i)
            if (i !== row && (a[i] >> BigInt(col)) & 1n) a[i] ^= a[row];
        ++row;
    }
    // Решта реалізації така сама, як вище
}

Go

// Рядки — бітові маски (uint64): біти 0..m-1 — коефіцієнти, біт m — вектор b.
// XOR двох масок (a[i] ^= a[row]) замінює побітову операцію над bitset.
func gaussMod(a []uint64, n, m int) (int, uint64) {
	where := make([]int, m)
	for i := range where {
		where[i] = -1
	}
	for col, row := 0, 0; col < m && row < n; col++ {
		for i := row; i < n; i++ {
			if (a[i]>>uint(col))&1 == 1 {
				a[i], a[row] = a[row], a[i]
				break
			}
		}
		if (a[row]>>uint(col))&1 == 0 {
			continue
		}
		where[col] = row

		for i := 0; i < n; i++ {
			if i != row && (a[i]>>uint(col))&1 == 1 {
				a[i] ^= a[row]
			}
		}
		row++
	}
	// Решта реалізації така сама, як вище
}

Оскільки ми використовуємо бітове стиснення, реалізація не лише коротша, а й у 32 рази швидша.

Невелике зауваження про різні евристики вибору опорного рядка

Загального правила, які евристики використовувати, не існує.

Евристика, використана в попередній реалізації, доволі добре працює на практиці. До того ж вона дає майже ті самі відповіді, що й «повний вибір опорного елемента» (full pivoting) — коли опорний рядок шукається серед усіх елементів усієї підматриці (від поточного рядка та поточного стовпця).

Утім, слід зауважити, що обидві евристики залежать від того, наскільки масштабовані початкові рівняння. Наприклад, якщо одне з рівнянь помножили на $10^6$, то це рівняння майже напевно буде обрано як опорне на першому кроці. Це здається доволі дивним, тож логічно перейти до складнішої евристики, яку називають неявним вибором опорного елемента (implicit pivoting).

Неявний вибір опорного елемента порівнює елементи так, ніби обидва рядки нормалізовано, щоб максимальний елемент дорівнював одиниці. Щоб реалізувати цю техніку, потрібно підтримувати максимум у кожному рядку (або підтримувати кожен рядок так, щоб максимум дорівнював одиниці, але це може призвести до зростання накопиченої похибки).

Покращення розв'язку

Попри різні евристики, алгоритм Гаусса-Жордана все одно може призводити до великих похибок у спеціальних матрицях навіть розміру $50 - 100$.

Тому отриманий розв'язок Гаусса-Жордана іноді потрібно покращити, застосувавши простий чисельний метод — наприклад, метод простої ітерації.

Таким чином, розв'язання стає двокроковим: спочатку застосовується алгоритм Гаусса-Жордана, а потім чисельний метод, що бере початковий розв'язок як розв'язок із першого кроку.

Задачі для практики

• Spoj - Xor Maximization
• Codechef - Knight Moving
• Lightoj - Graph Coloring
• UVA 12910 - Snakes and Ladders
• TIMUS1042 Central Heating
• TIMUS1766 Humpty Dumpty
• TIMUS1266 Kirchhoff's Law
• Codeforces - No game no life

Відеоматеріали