En esta nota revisamos métodos por bloques para resolver sistemas de ecuaciones lineales:
En general son de la forma: $$\begin{array}{ccc} a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n &= & b_1 \\ a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n &= & b_2 \\ \vdots & & \\ a_{m1}x_1 + a_{m2}x_2 + \cdots + a_{mn}x_n &=& b_m \end{array}$$
donde: las $x_i$'s son las incógnitas y las $a_i$'s y $b_i$'s son constantes conocidas.
Las entradas $a_{ij}$'s son llamadas coeficientes del sistema y el conjunto de $b_i$'s se le llama lado derecho del sistema. Si todas las $b_i$'s son iguales a $0$ el sistema se le nombra homogéneo.
3 posibilidades para solución del sistema anterior:
Una única solución: sólo existe uno y sólo un conjunto de valores de $x_i$'s que satisfacen todas las ecuaciones simultáneamente y el sistema se nombra consistente o no singular.
Ninguna solución: no existe ningún conjunto de valores de $x_i$'s que satisfacen todas las ecuaciones simultáneamente (el conjunto solución es vacío) y el sistema se nombra inconsistente o singular.
Infinitas soluciones: hay una infinidad de conjuntos (distintos) de valores de las $x_i$'s que satisfacen todas las ecuaciones simultáneamente. obs: si un sistema tiene más de una solución entonces tiene una infinidad de soluciones y el sistema se nombra consistente o no singular.
Resolver un sistema de ecuaciones lineales equivale a encontrar la intersección entre rectas, planos o hiperplanos (2,3 o n dimensiones respectivamente). Por ejemplo para un caso de dos dimensiones se tiene:
El inciso a) representa un sistema de ecuaciones lineales sin solución, el inciso b) infinitas soluciones (en el dibujo ligeramente se desplazó hacia abajo una de las rectas para mostrar ambas) y el inciso c) una única solución.
Nombres dados a los sistemas dependiendo de sus dimensiones:
Si $m > n$ el sistema $Ax=b$ se nombra overdetermined: hay más ecuaciones que incógnitas. En general no existe $x \in \mathbb{R}^n$ que satisfaga el sistema por lo que se busca resolver el problema de encontrar $x \in \mathbb{R}^n$ tal que $||Ax-b||$ sea mínima. Si $||\cdot||$ es la norma Euclidiana, el problema anterior se le llama mínimos cuadrados.
Si $m < n$ el sistema $Ax=b$ se nombra underdetermined: hay menos ecuaciones que incógnitas. El sistema tiene una infinidad de soluciones ó ninguna solución. Lo anterior depende del $rank(A)$. Si $rank(A)=m$ entonces $Ax=b$ tiene infinitas soluciones*. Típicamente se busca resolver el sistema para $x$ que tiene norma mínima (lo cual conduce a plantear un problema de optimización numérica).
*Para ver esto considérese el resultado del álgebra lineal Rank-nullity theorem que en este contexto se escribe como $n = rank(A) + nullity(A) = m + nullity(A)$ $\therefore nullity(A) = n-m$. Si $x_0 \in \mathbb{R}^n$ es una solución de $Ax=b$ entonces $x_0 + N(A)$ es el conjunto de soluciones de $Ax=b$, donde $x_0 + N(A) = \{x_0 + z | z \in N(A)\}$.
Existen una gran cantidad de métodos para resolver SEL. Típicamente se elige el método de acuerdo a los coeficientes de la matriz del sistema que determinan propiedades de la misma y sus dimensiones.
Son SEL cuya matriz del sistema es triangular inferior o superior. Un SEL triangular inferior se resuelve con el método de sustitución hacia delante que implica realizar $\mathcal{O}(n^2)$ flops. Si el SEL es triangular superior se tiene el método de sustitución hacia atrás con la misma cantidad de operaciones realizadas que el de hacia delante.
Para SEL más generales (no tienen estructura identificable) se tienen los métodos iterativos y directos o basados en factorizaciones matriciales (FM). Entre los directos o basados en FM se encuentran:
y como ejemplo de los iterativos están:
Ambos métodos: iterativos y directos o basados en FM encuentran sistemas de ecuaciones equivalentes a partir de operaciones básicas del álgebra lineal. Dos sistemas de ecuaciones lineales son equivalentes si tienen el mismo conjunto solución.
Los métodos o algoritmos iterativos anteriores no realizan factorizaciones matriciales. Hay características de la matriz (por ejemplo simetría) que no explotan en su implementación pero se benefician de otras (por ejemplo si $A$ es sparse o rala: muchos ceros). Otras características de tales métodos son: uso de menos memoria que los directos o basados en FM, convergen bajo condiciones simples de verificar y son sencillos de implementar.
Los métodos directos o basados en FM en general consumen más memoria que los iterativos, tienen el potencial de utilizar operaciones de nivel BLAS $3$, siempre encuentran solución si la matriz del sistema es no singular y en general requieren mejor dominio de la programación para su implementación que los iterativos. Una vez calculada la FM de la matriz $A$ se utilizan métodos para resolver SEL que resultan de la FM. Un ejemplo de tales métodos son los de sustitución hacia delante y hacia atrás para resolver SEL triangulares.
Nombremos a los métodos o algoritmos numéricos anteriores nuestros Métodos o Algoritmos Básicos del Álgebra Lineal para resolver SEL o FM (MBAL/ABAL para SEL o FM). Cada uno puede utilizarse para resolver SEL, FM o para otros propósitos. Por ejemplo la factorización QR también se utiliza en el cálculo de eigenvalores y eigenvectores en el algoritmo QR, ver QR algorithm.
Los MBAL/ABAL para SEL o FM sirven de apoyo para resolver SEL por bloques.
Aquí se presenta una tabla resúmen que muestra el número de flops si se desea resolver $Ax=b$ con $A \in \mathbb{R}^{n \times n}$:
Como se revisó en la nota 3.1.Vectorizacion_BLAS_y_el_uso_del_cache_eficientemente los blocking algorithms son ricos en operaciones de nivel 3 de BLAS, simplifican notación matemática y su escritura ayuda a pensar en algoritmos paralelizables. Una desventaja es la cantidad de memoria RAM que utilizan vs métodos que no trabajan por bloques.
Este método consiste en eliminar un subconjunto de variables y resolver un sistema más pequeño. Si el sistema más pequeño asociado a la submatriz puede resolverse con alguno de los MBAL/ABAL para SEL o FM anteriores entonces este método puede tener más eficiencia que uno que no trabaja por bloques.
Considérese un sistema de la forma: $Ax=b$ escrito como:
con $x_1 \in \mathbb{R}^{n_1}, x_2 \in \mathbb{R}^{n_2}, A_{11} \in \mathbb{R}^{n_1 \times n_1}, A_{22} \in \mathbb{R}^{n_2 \times n_2}, b_1 \in \mathbb{R}^{n_1}, b_2 \in \mathbb{R}^{n_2}$.
El método de eliminación por bloques consiste en que si $A_{11}$ es invertible, entonces se puede eliminar $x_1$ de las ecuaciones como sigue:
De la primer ecuación por bloques: $A_{11}x_1 + A_{12}x_2 = b_1$ se obtiene: $x_1 = A_{11}^{-1}(b_1-A_{12}x_2)$.
Sustituyendo esta relación en la segunda ecuación por bloques obtenemos la ecuación reducida: $(A_{22}-A_{21}A_{11}^{-1}A_{12})x_2 = b_2 - A_{21}A_{11}^{-1}b_1$.
Obs: estas dos ecuaciones anteriores son equivalentes al sistema original $Ax=b$. La matriz $A_{22}-A_{21}A_{11}^{-1}A_{12}$ se le nombra complemento de Schur del bloque $A_{11}$ en $A$ (ver Schur complement) esto es: $S = A_{22}-A_{21}A_{11}^{-1}A_{12}$. Se tiene la propiedad siguiente: $S$ es no singular si y sólo si $A$ es no singular.
Por lo anterior el método de eliminación por bloque consiste:
Algoritmo
Sean $A$ y $A_{11}$ no singulares.
1) Calcular $A_{11}^{-1}A_{12}$ y $A_{11}^{-1}b_1$ teniendo cuidado en no calcular la inversa sino un sistema de ecuaciones lineales:
2) Calcular el complemento de Schur del bloque $A_{11}$ en $A$: $S = A_{22}-A_{21}A_{11}^{-1}A_{12}$. Calcular $ \hat{b} = b_2-A_{21}A_{11}^{-1}b_1$.
3) Resolver $Sx_2 = \hat{b}$.
4) Resolver $A_{11}x_1 = b_1-A_{12}x_2$.
Comentario: En los pasos que se requieren resolver sistemas de ecuaciones lineales se utilizan los MBAL/ABAL para SEL o FM.
El método de factor solve consiste en expresar a la matriz del sistema $Ax = b$ como un producto de la forma: $A = A_1A_2 \cdots A_k$. Si $A$ es no singular entonces: $x = A^{-1}b = A_k^{-1}A_{k-1}^{-1} \cdots A_1^{-1}b$.
De esta forma, $x$ puede calcularse de "derecha a izquierda":
Algoritmo
Comentarios:
El paso que consiste en factorizar a la matriz $A$ es el paso nombrado factor y el usado en resolver el sistema $A_iz_i = z_{i-1}$ es nombrado solve.
Normalmente $A_i$ tiene una estructura aprovechable (por ejemplo triangular superior/inferior) para resolver el sistema $A_iz_i = z_{i-1}$. Por ejemplo en el caso de la factorización $LU$ de $A$ se cumple $PA = LU$ con $P$ matriz de permutación, $L$ triangular inferior con $1$'s en la diagonal y $U$ triangular superior. Entonces:
El método de block elimination puede escribirse como un factor-solve method:
Paso $1$:
$$\begin{array}{l} \left[ \begin{array}{cc} A_{11} & A_{12}\\ A_{21} & A_{22} \end{array} \right]= \left[\begin{array}{cc} A_{11} & 0 \\ A_{21} & S \end{array} \right] \left[\begin{array}{cc} I & A_{11}^{-1}A_{12} \\ 0 & I \end{array} \right] \end{array} $$Paso $2$. Resolver por el método de sustitución hacia delante por bloques:
$$\begin{array}{l} \left[ \begin{array}{cc} A_{11} & 0\\ A_{21} & S \end{array} \right] \cdot \left[\begin{array}{c} z_1 \\ z_2 \end{array} \right] = \left[\begin{array}{c} b_1\\ b_2 \end{array} \right] \end{array} $$Paso $3$. Resolver por el método de sustitución hacia atrás por bloques:
$$\begin{array}{l} \left[ \begin{array}{cc} I & A_{11}^{-1}A_{12}\\ 0 & I \end{array} \right] \cdot \left[\begin{array}{c} x_1 \\ x_2 \end{array} \right] = \left[\begin{array}{c} z_1\\ z_2 \end{array} \right] \end{array} $$Comentario: En los pasos que se requieren resolver sistemas de ecuaciones lineales se utilizan los MBAL/ABAL para SEL o FM.
Referencias:
G. H. Golub, C. F. Van Loan, Matrix Computations, John Hopkins University Press, 2013.
Ver 1_ecuaciones_lineales para una introducción a sistemas de ecuaciones lineales y factorizaciones matriciales con Python3.