Método de Gauss para resolver sistemas de ecuaciones lineales

Se da un sistema de $n$ ecuaciones algebraicas lineales (SLAE) con $m$ incógnitas. Se pide resolver el sistema: determinar si no tiene solución, si tiene exactamente una solución o si tiene un número infinito de soluciones. Y en caso de que tenga al menos una solución, hallar cualquiera de ellas.

Formalmente, el problema se formula así: resolver el sistema:

$\begin{aligned} a_{11} x_{1} + a_{12} x_{2} + & \dots + a_{1 m} x_{m} = b_{1} \\ a_{21} x_{1} + a_{22} x_{2} + & \dots + a_{2 m} x_{m} = b_{2} \\ ⋮ \\ a_{n 1} x_{1} + a_{n 2} x_{2} + & \dots + a_{n m} x_{m} = b_{n} \end{aligned}$

donde los coeficientes $a_{ij}$ (para $i$ de 1 a $n$ , $j$ de 1 a $m$ ) y $b_i$ ( $i$ de 1 a $n$ ) son conocidos y las variables $x_i$ ( $i$ de 1 a $m$ ) son incógnitas.

Este problema también tiene una representación matricial sencilla:

$Ax = b,$

donde $A$ es una matriz de tamaño $n \times m$ de coeficientes $a_{ij}$ y $b$ es el vector columna de tamaño $n$ .

Cabe señalar que el método presentado en este artículo también se puede usar para resolver la ecuación módulo cualquier número $p$ , es decir:

$\begin{aligned} a_{11} x_{1} + a_{12} x_{2} + & \dots + a_{1 m} x_{m} \equiv b_{1} (m o d p) \\ a_{21} x_{1} + a_{22} x_{2} + & \dots + a_{2 m} x_{m} \equiv b_{2} (m o d p) \\ ⋮ \\ a_{n 1} x_{1} + a_{n 2} x_{2} + & \dots + a_{n m} x_{m} \equiv b_{n} (m o d p) \end{aligned}$

Gauss

Estrictamente hablando, el método que se describe a continuación debería llamarse “Gauss-Jordan”, o eliminación de Gauss-Jordan, porque es una variante del método de Gauss, descrita por Jordan en 1887.

Visión general

El algoritmo es una eliminación secuencial de las variables en cada ecuación, hasta que cada ecuación tenga solo una variable restante. Si $n = m$ , se puede pensar como transformar la matriz $A$ en la matriz identidad, y resolver la ecuación en este caso obvio, donde la solución es única y es igual al coeficiente $b_i$ .

La eliminación gaussiana se basa en dos transformaciones sencillas:

Es posible intercambiar dos ecuaciones
Cualquier ecuación se puede reemplazar por una combinación lineal de esa fila (con coeficiente no nulo) y de algunas otras filas (con coeficientes arbitrarios).

En el primer paso, el algoritmo de Gauss-Jordan divide la primera fila por $a_{11}$ . Luego, el algoritmo suma la primera fila a las filas restantes de modo que los coeficientes de la primera columna se vuelvan todos ceros. Para lograrlo, en la $i$ -ésima fila debemos sumar la primera fila multiplicada por $- a_{i1}$ . Nótese que esta operación también debe realizarse sobre el vector $b$ . En cierto sentido, se comporta como si el vector $b$ fuera la $(m+1)$ -ésima columna de la matriz $A$ .

Como resultado, después del primer paso, la primera columna de la matriz $A$ consistirá en $1$ en la primera fila, y $0$ en las demás filas.

De manera similar, realizamos el segundo paso del algoritmo, donde consideramos la segunda columna de la segunda fila. Primero, la fila se divide por $a_{22}$ , luego se resta de las otras filas de modo que toda la segunda columna se vuelva $0$ (excepto la segunda fila).

Continuamos este proceso para todas las columnas de la matriz $A$ . Si $n = m$ , entonces $A$ se convertirá en la matriz identidad.

Búsqueda del elemento pivote

El esquema descrito dejó fuera muchos detalles. En el $i$ -ésimo paso, si $a_{ii}$ es cero, no podemos aplicar directamente el método descrito. En su lugar, primero debemos seleccionar una fila pivote: hallar una fila de la matriz donde la $i$ -ésima columna no sea cero, y luego intercambiar las dos filas.

Nótese que aquí intercambiamos filas pero no columnas. Esto se debe a que si se intercambian columnas, entonces al hallar una solución hay que recordar intercambiar de vuelta a los lugares correctos. Así, intercambiar filas es mucho más fácil.

En muchas implementaciones, cuando $a_{ii} \neq 0$ , se ve que igual se intercambia la $i$ -ésima fila con alguna fila pivote, usando alguna heurística como elegir la fila pivote con máximo valor absoluto de $a_{ji}$ . Esta heurística se usa para reducir el rango de valores de la matriz en pasos posteriores. Sin esta heurística, incluso para matrices de tamaño alrededor de $20$ , el error será demasiado grande y puede causar desbordamiento en los tipos de punto flotante de C++.

Casos degenerados

En el caso en que $m = n$ y el sistema es no degenerado (es decir, tiene determinante no nulo, y tiene solución única), el algoritmo descrito arriba transformará $A$ en la matriz identidad.

Ahora consideramos el caso general, donde $n$ y $m$ no son necesariamente iguales, y el sistema puede ser degenerado. En estos casos, el elemento pivote en el $i$ -ésimo paso puede no hallarse. Esto significa que en la $i$ -ésima columna, a partir de la línea actual, todas contienen ceros. En este caso, o bien no hay ningún valor posible de la variable $x_i$ (lo que significa que el SLAE no tiene solución), o $x_i$ es una variable independiente y puede tomar un valor arbitrario. Al implementar Gauss-Jordan, hay que continuar el trabajo para las variables posteriores y simplemente saltar la $i$ -ésima columna (esto es equivalente a eliminar la $i$ -ésima columna de la matriz).

Así, algunas de las variables en el proceso pueden resultar independientes. Cuando el número de variables, $m$ , es mayor que el número de ecuaciones, $n$ , entonces se hallarán al menos $m - n$ variables independientes.

En general, si se halla al menos una variable independiente, puede tomar cualquier valor arbitrario, mientras que las otras (dependientes) se expresan a través de ella. Esto significa que cuando trabajamos en el cuerpo de los números reales, el sistema potencialmente tiene infinitas soluciones. Pero hay que recordar que cuando hay variables independientes, el SLAE puede no tener solución en absoluto. Esto ocurre cuando las ecuaciones restantes sin tratar tienen al menos un término constante no nulo. Se puede comprobar esto asignando ceros a todas las variables independientes, calculando las otras variables, y luego sustituyendo en el SLAE original para verificar si lo satisfacen.

Implementación

A continuación hay una implementación de Gauss-Jordan. La elección de la fila pivote se hace con una heurística: elegir el valor máximo en la columna actual.

La entrada de la función gauss es la matriz del sistema $a$ . La última columna de esta matriz es el vector $b$ .

La función devuelve el número de soluciones del sistema $(0, 1,\textrm{or } \infty)$ . Si existe al menos una solución, entonces se devuelve en el 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;
}

Notas de implementación:

La función usa dos punteros: la columna actual $col$ y la fila actual $row$ .
Para cada variable $x_i$ , el valor $where(i)$ es la línea donde esta columna no es cero. Este vector es necesario porque algunas variables pueden ser independientes.
En esta implementación, la $i$ -ésima línea actual no se divide por $a_{ii}$ como se describió arriba, así que al final la matriz no es la matriz identidad (aunque aparentemente dividir la $i$ -ésima línea puede ayudar a reducir errores).
Después de hallar una solución, se inserta de vuelta en la matriz, para comprobar si el sistema tiene al menos una solución o no. Si la prueba de la solución es exitosa, entonces la función devuelve 1 o $\inf$ , según haya al menos una variable independiente o no.

Complejidad

Ahora debemos estimar la complejidad de este algoritmo. El algoritmo consiste en $m$ fases; en cada fase:

Buscar y reordenar la fila pivote. Esto toma $O(n + m)$ cuando se usa la heurística mencionada arriba.
Si se halla el elemento pivote en la columna actual, entonces debemos sumar esta ecuación a todas las demás ecuaciones, lo que toma tiempo $O(nm)$ .

Así, la complejidad final del algoritmo es $O(\min (n, m) . nm)$ . En el caso $n = m$ , la complejidad es simplemente $O(n^3)$ .

Nótese que cuando el SLAE no está sobre números reales, sino módulo dos, entonces el sistema se puede resolver mucho más rápido, lo que se describe a continuación.

Aceleración del algoritmo

La implementación anterior se puede acelerar dos veces, dividiendo el algoritmo en dos fases: directa e inversa:

Fase directa: similar a la implementación anterior, pero la fila actual solo se suma a las filas posteriores. Como resultado, obtenemos una matriz triangular en lugar de diagonal.
Fase inversa: cuando la matriz es triangular, primero calculamos el valor de la última variable. Luego sustituimos este valor para hallar el valor de la siguiente variable. Luego sustituimos estos dos valores para hallar las siguientes variables…

La fase inversa solo toma $O(nm)$ , que es mucho más rápido que la fase directa. En la fase directa, reducimos el número de operaciones a la mitad, reduciendo así el tiempo de ejecución de la implementación.

Resolver SLAE modular

Para resolver un SLAE en algún módulo, todavía podemos usar el algoritmo descrito. Sin embargo, en el caso de que el módulo sea igual a 2, podemos realizar la eliminación de Gauss-Jordan de forma mucho más efectiva usando operaciones bit a bit y los tipos de datos bitset de 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;
	}
        // The rest of implementation is the same as above
}

Como usamos compresión de bits, la implementación no solo es más corta, sino también 32 veces más rápida.

Una nota breve sobre distintas heurísticas de elección de la fila pivote

No hay una regla general sobre qué heurísticas usar.

La heurística usada en la implementación anterior funciona bastante bien en la práctica. También resulta dar casi las mismas respuestas que el “pivoteo completo” (full pivoting), donde la fila pivote se busca entre todos los elementos de toda la submatriz (desde la fila actual y la columna actual).

Aunque hay que notar que ambas heurísticas dependen de cuánto se hayan escalado las ecuaciones originales. Por ejemplo, si una de las ecuaciones se multiplicó por $10^6$ , entonces esta ecuación es casi segura de ser elegida como pivote en el primer paso. Esto parece bastante extraño, así que parece lógico cambiar a una heurística más complicada, llamada pivoteo implícito (implicit pivoting).

El pivoteo implícito compara elementos como si ambas líneas estuvieran normalizadas, de modo que el elemento máximo sería la unidad. Para implementar esta técnica, hay que mantener el máximo en cada fila (o mantener cada línea de modo que el máximo sea la unidad, pero esto puede llevar a un aumento del error acumulado).

Mejorar la solución

A pesar de varias heurísticas, el algoritmo de Gauss-Jordan todavía puede producir errores grandes en matrices especiales incluso de tamaño $50 - 100$ .

Por lo tanto, la solución resultante de Gauss-Jordan a veces debe mejorarse aplicando un método numérico sencillo, por ejemplo, el método de iteración simple.

Así, la solución se convierte en dos pasos: primero se aplica el algoritmo de Gauss-Jordan, y luego un método numérico que toma como solución inicial la solución del primer paso.