Dado un sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\), un método iterativo para resolver dicho sistema consiste en hallar una sucesión de soluciones aproximadas \((\mathbf{x}^{(k)})_{k=0}^\infty\) tal que \(\displaystyle\lim_{k\to\infty}\mathbf{x}^{(k)}=\mathbf{x}\) donde \(\mathbf{x}\) es la solución del sistema.

Cuando escribimos que \(\displaystyle\lim_{k\to\infty}\mathbf{x}^{(k)}=\mathbf{x}\), queremos decir que dada una norma matricial \(\|\cdot\|\), \(\displaystyle\lim_{k\to\infty}\|\mathbf{x}^{(k)}-\mathbf{x}\|=0\).

Los algoritmos numéricos basado en métodos iterativos tienen las propiedades siguientes:

• obtienen una solución aproximada, al contrario de lo que pasaba con los métodos directos que hallaban una solución exacta del sistema lineal. Es decir, suponiendo que no hay errores en los datos ni en las operaciones aritméticas, la solución obtenida por un método iterativo es aproximada.

• el número de pasos que hay que realizar para hallar la solución de un sistema lineal no está predeterminado de entrada.

Tenemos entonces que establecer un criterio de parada de cara a obtener una solución aproximada del sistema lineal dado un umbral de error \(\epsilon\).

Concretamente, sea \((\mathbf{x}^{(k)})_{k=0}^\infty\) la sucesión de soluciones aproximadas cuyo límite es la solución exacta \(\mathbf{x}\) del sistema lineal, \(\|\cdot\|\) una norma vectorial y \(\epsilon\) un umbral de error.

Existen dos criterios de parada:

• hallar el entero \(k\) tal que \(\|\mathbf{x}^{(k)}-\mathbf{x}^{(k-1)}\|<\epsilon\), (criterio del error absoluto)
• hallar el entero \(k\) tal que \(\frac{\|\mathbf{x}^{(k)}-\mathbf{x}^{(k-1)}\|}{\|\mathbf{x}^{(k)}\|}<\epsilon\), (criterio del error relativo)

y adoptamos \(\mathbf{x}^{(k)}\) como solución aproximada del sistema lineal.

Los métodos iterativos que vamos a ver construyen la sucesión de aproximaciones a la solución de forma iterativa usando una expresión de la forma siguiente: \[ \mathbf{x}^{(k)}=\mathbf{T}\mathbf{x}^{(k-1)}+\mathbf{c}, \] donde \(\mathbf{T}\) es la llamada matriz de iteración y \(\mathbf{c}\) el vector independiente del método iterativo en cuestión.

El método iterativo necesita una aproximación inicial \(\mathbf{x}^{(0)}\) que podemos suponer que vale \(\mathbf{0}\) o cualquier valor que queramos.

Veremos posteriomente que la convergencia del método iterativo depende sólo de la matriz de iteración \(\mathbf{T}\).

¿Cómo se construye la expresión anterior que nos define la sucesión de aproximaciones \(\mathbf{x}^{(k)}\)?

La idea clave es transformar el sistema a resolver \(\mathbf{A}\mathbf{x}=\mathbf{b}\) en un sistema equivalente de la forma \(\mathbf{x}=\mathbf{T}\mathbf{x}+\mathbf{c}\) ya que si la sucesión de aproximaciones \(\mathbf{x}^{(k)}\) converge a la solución del sistema \(\mathbf{x}\), dicha solución debe verificar \(\mathbf{x}=\mathbf{T}\mathbf{x}+\mathbf{c}\).

Supongamos que los coeficientes \(a_{ii}\) de la matriz \(\mathbf{A}=(a_{ij})_{i,j=1,\ldots,n}\) del sistema a resolver \(\mathbf{A}\mathbf{x}=\mathbf{b}\) son todos diferentes de cero.

Fijaos que siempre existe una permutación de filas para que dicha condición se cumpla ya que, en caso contrario, la matriz del sistema sería singular y el sistema no tendría solución única.

Entonces para transformar el sistema \(\mathbf{A}\mathbf{x}=\mathbf{b}\) en un sistema equivalente de la forma \(\mathbf{x}=\mathbf{T}\mathbf{x}+\mathbf{c}\), podemos aislar las incógnitas \(x_i\) de la ecuación \(i\)-ésima de la forma siguiente: \[ \begin{align*} E_i:& a_{i1}x_1+a_{i2}x_2+\cdots +a_{ii}x_i+\cdots +a_{in}x_n=b_i,\ \Rightarrow\\ & x_i=\frac{1}{a_{ii}}(b_i-(a_{i1}x_1+\cdots+a_{i,i-1}x_{i-1}+a_{i,i+1}x_{i+1}+\cdots+a_{in}x_n)) \end{align*} \]

Entonces, la matriz de iteración \(\mathbf{T}_J\) y el vector independiente \(\mathbf{c}_J\) serán los siguientes:

\[ \mathbf{T}_J=\begin{bmatrix}0&-\frac{a_{12}}{a_{11}}&\ldots&-\frac{a_{1n}}{a_{11}} \\-\frac{a_{21}}{a_{22}}&0&\ldots&-\frac{a_{2n}}{a_{22}} \\\vdots&\vdots&\ddots&\vdots \\-\frac{a_{n1}}{a_{nn}}&-\frac{a_{n2}}{a_{nn}}&\ldots&0 \\\end{bmatrix},\quad \mathbf{c}_J=\begin{bmatrix}\frac{b_{1}}{a_{11}} \\\frac{b_2}{a_{22}} \\\vdots \\\frac{b_n}{a_{nn}} \\\end{bmatrix}. \]

La sucesión \(\mathbf{x}^{(k)}\) se generará de la forma siguiente: \[ \begin{bmatrix}x_{1}^{(k)} \\x_2^{(k)} \\\vdots \\x_n^{(k)} \\\end{bmatrix}=\begin{bmatrix}0&-\frac{a_{12}}{a_{11}}&\ldots&-\frac{a_{1n}}{a_{11}} \\-\frac{a_{21}}{a_{22}}&0&\ldots&-\frac{a_{2n}}{a_{22}} \\\vdots&\vdots&\ddots&\vdots \\-\frac{a_{n1}}{a_{nn}}&-\frac{a_{n2}}{a_{nn}}&\ldots&0 \\\end{bmatrix}\begin{bmatrix}x_{1}^{(k-1)} \\x_2^{(k-1)} \\\vdots \\x_n^{(k-1)} \\\end{bmatrix}+\begin{bmatrix}\frac{b_{1}}{a_{11}} \\\frac{b_2}{a_{22}} \\\vdots \\\frac{b_n}{a_{nn}} \\\end{bmatrix} \]

Vamos a escribir la matriz de iteración del método de Jacobi de otra forma.

Descomponemos la matriz del sistema \(\mathbf{A}=(a_{ij})_{i,j=1,\ldots,n}\) de la forma siguiente: \[ \mathbf{A}=\mathbf{D}+\mathbf{L}+\mathbf{U}, \] donde las matrices \(\mathbf{D}\), \(\mathbf{L}\) y \(\mathbf{U}\) contienen los elementos diagonales, los elementos por debajo de la diagonal y por encima de la diagonal respectivamente:

\[ \mathbf{D}=\begin{bmatrix}a_{11}&0&\ldots&0 \\0&a_{22}&\ldots&0 \\\vdots&\vdots&\ddots&\vdots \\0&0&\ldots&a_{nn} \\\end{bmatrix},\ \mathbf{L}=\begin{bmatrix}0&0&\ldots&0 \\a_{21}&0&\ldots&0 \\\vdots&\vdots&\ddots&\vdots \\a_{n1}&a_{n2}&\ldots&0 \\\end{bmatrix},\ \mathbf{U}=\begin{bmatrix}0&a_{12}&\ldots&a_{1n} \\0&0&\ldots&a_{2n} \\\vdots&\vdots&\ddots&\vdots \\0&0&\ldots&0 \\\end{bmatrix}. \]

El sistema \(\mathbf{A}\mathbf{x}=\mathbf{b}\) puede escribirse como: \[ (\mathbf{D}+\mathbf{L}+\mathbf{U})\mathbf{x}=\mathbf{b}. \] Vamos a deducir el método de Jacobi a partir de las matrices anteriores.

• Dejamos los elementos diagonales y “pasamos” los demás elementos al término de la derecha: \[ \mathbf{D}\mathbf{x}=-(\mathbf{L}+\mathbf{U})\mathbf{x}+\mathbf{b}. \]
• Despejamos el vector \(\mathbf{x}\): \[ \mathbf{x}=-\mathbf{D}^{-1}(\mathbf{L}+\mathbf{U})\mathbf{x}+\mathbf{D}^{-1}\mathbf{b}. \]

La sucesión de aproximaciones será: \[ \mathbf{x}^{(k)}=-\mathbf{D}^{-1}(\mathbf{L}+\mathbf{U})\mathbf{x}^{(k-1)}+\mathbf{D}^{-1}\mathbf{b}. \] Entonces, la matriz de iteración y el vector independiente del método de Jacobi serán: \[ \mathbf{T}_J=-\mathbf{D}^{-1}(\mathbf{L}+\mathbf{U}),\quad \mathbf{c}_J=\mathbf{D}^{-1}\mathbf{b}. \]

Vamos a dar el pseudocódigo del método de Jacobi usando como criterio de parada el criterio del error absoluto.

Para tener el pseudocódigo usando como criterio de parada el criterio del error relativo basta cambiar la línea de código

\(\|\mathbf{x}-\mathbf{X0}\| <\) TOL then

por

\(\frac{\|\mathbf{x}-\mathbf{X0}\|}{\|\mathbf{x}\|} <\) TOL then.

• INPUT matriz del sistema \(\mathbf{A}=(a_{ij})_{i=1,\ldots,n,j=1,\ldots,n}\), vector de términos independientes \(\mathbf{b}=(b_i)_{i=1,\ldots,n}\), valor inicial \(\mathbf{X0}=\mathbf{x}^{(0)}\), error absoluto o tolerancia TOL, número máximo de iteraciones Nmax.
• Set k=1.

• While k <= Nmax
  • For i=1,...,n
    • Set \(x_i=\frac{1}{a_{ii}}\left(-\sum_{j=1,j\neq i}^n a_{ij}X0_j+b_i\right)\).
  • If \(\|\mathbf{x}-\mathbf{X0}\| <\) TOL then
    • Print \(x_1,\ldots,x_n\)
    • STOP.
  • Set k=k+1.
  • For i=1,...,n
    • Set \(X0_i =x_i\).
• Print número máximo de iteraciones alcanzado. El método no converge.
• STOP.

El método de Gauss-Seidel consiste en aplicar el método de Jacobi con una modificación: las componentes del vector \(\mathbf{x}^{(k)}\) de la iteración actual, \(x_i^{(k)}\), que se van calculando, se usan para el cálculo de las componentes \(x_{i+1}^{(k)},\ldots,x_n^{(k)}\) de dicho vector.

Es decir, la sucesión de aproximaciones se calcula de la forma siguiente: \[ \begin{align*} x_1^{(k)} & =\frac{1}{a_{11}}(b_1-(a_{12}x_2^{(k-1)}+\cdots +a_{1n}x_n^{(k-1)})),\\ x_2^{(k)} & =\frac{1}{a_{22}}(b_2-(a_{21}x_1^{(k)}+a_{23}x_3^{(k-1)}+\cdots +a_{2n}x_n^{(k-1)})),\\ \vdots & \\ x_i^{(k)}& =\frac{1}{a_{ii}}(b_i-(a_{i1}x_1^{(k)}+\cdots +a_{i,i-1}x_{i-1}^{(k)}+a_{i,i+1}x_{i+1}^{(k-1)}+\cdots +a_{in}x_n^{(k-1)})),\\ \vdots & \\ x_n^{(k)} & =\frac{1}{a_{nn}}(b_n-(a_{n1}x_1^{(k)}+\cdots +a_{n,n-1}x_{n-1}^{(k)})). \end{align*} \]

Matricialmente, la sucesión de aproximaciones se calcularía de la forma siguiente:

\[ \begin{bmatrix}x_{1}^{(k)} \\x_2^{(k)} \\\vdots \\x_n^{(k)} \\\end{bmatrix}=\begin{bmatrix}0&0&\ldots&0 \\-\frac{a_{21}}{a_{22}}&0&\ldots&0 \\\vdots&\vdots&\ddots&\vdots \\-\frac{a_{n1}}{a_{nn}}&-\frac{a_{n2}}{a_{nn}}&\ldots&0 \\\end{bmatrix}\begin{bmatrix}x_{1}^{(k)} \\x_2^{(k)} \\\vdots \\x_n^{(k)} \\\end{bmatrix}+\begin{bmatrix}0&-\frac{a_{12}}{a_{11}}&\ldots&-\frac{a_{1n}}{a_{11}} \\0&0&\ldots&-\frac{a_{2n}}{a_{22}} \\\vdots&\vdots&\ddots&\vdots \\0&0&\ldots&0 \\\end{bmatrix}\begin{bmatrix}x_{1}^{(k-1)} \\x_2^{(k-1)} \\\vdots \\x_n^{(k-1)} \\\end{bmatrix}+\begin{bmatrix}\frac{b_{1}}{a_{11}} \\\frac{b_2}{a_{22}} \\\vdots \\\frac{b_n}{a_{nn}} \\\end{bmatrix}. \]

Llamando \(\mathbf{T}_1\) y \(\mathbf{T}_2\) a las matrices siguientes: \[ \mathbf{T}_1=\begin{bmatrix}0&0&\ldots&0 \\-\frac{a_{21}}{a_{22}}&0&\ldots&0 \\\vdots&\vdots&\ddots&\vdots \\-\frac{a_{n1}}{a_{nn}}&-\frac{a_{n2}}{a_{nn}}&\ldots&0 \\\end{bmatrix},\ \mathbf{T}_2=\begin{bmatrix}0&-\frac{a_{12}}{a_{11}}&\ldots&-\frac{a_{1n}}{a_{11}} \\0&0&\ldots&-\frac{a_{2n}}{a_{22}} \\\vdots&\vdots&\ddots&\vdots \\0&0&\ldots&0 \\\end{bmatrix}, \] el método de Gauss-Seidel puede escribirse de la forma siguiente: \[ \mathbf{x}^{(k)}=\mathbf{T}_1\mathbf{x}^{(k)}+\mathbf{T}_2\mathbf{x}^{(k-1)}+\mathbf{c}_J, \]

o, si se quiere: \[ \begin{align*} (\mathbf{Id}-\mathbf{T}_1)\mathbf{x}^{(k)} & =\mathbf{T}_2\mathbf{x}^{(k-1)}+\mathbf{c}_J,\ \Rightarrow\\ \mathbf{x}^{(k)} & =(\mathbf{Id}-\mathbf{T}_1)^{-1}\mathbf{T}_2\mathbf{x}^{(k-1)}+(\mathbf{Id}-\mathbf{T}_1)^{-1}\mathbf{c}_J. \end{align*} \] La matriz de iteración del método de Gauss Seidel será: \[ \mathbf{T}_{GS}=(\mathbf{Id}-\mathbf{T}_1)^{-1}\mathbf{T}_2, \] y el vector independiente será: \[ \mathbf{c}_{GS}=(\mathbf{Id}-\mathbf{T}_1)^{-1}\mathbf{c}_J. \]

Observación.

La matriz \(\mathbf{Id}-\mathbf{T}_1\) es invertible ya que: \[ \mathbf{Id}-\mathbf{T}_1=\begin{bmatrix}1&0&\ldots&0 \\\frac{a_{21}}{a_{22}}&1&\ldots&0 \\\vdots&\vdots&\ddots&\vdots \\\frac{a_{n1}}{a_{nn}}&\frac{a_{n2}}{a_{nn}}&\ldots&1 \\\end{bmatrix}, \] es triangular inferior con unos en la diagonal. Su determinante será, por tanto, \(1\).

Deduzcamos la matriz de iteración y el vector independiente del método de Gauss-Seidel a partir de la descomposición de la matriz del sistema \(\mathbf{A}\) que vimos anteriormente: \[ \mathbf{A}=\mathbf{D}+\mathbf{L}+\mathbf{U}. \] El sistema que tenemos que resolver és \(\mathbf{A}\mathbf{x}=\mathbf{b}\) o, si se quiere: \((\mathbf{D}+\mathbf{L}+\mathbf{U})\mathbf{x}=\mathbf{b}\):

• Dejamos los elementos diagonales y “pasamos” los demás elementos al término de la derecha: \[ \mathbf{D}\mathbf{x}=-(\mathbf{L}+\mathbf{U})\mathbf{x}+\mathbf{b}. \]
• “Separamos” las componentes calculadas de las no calculadas del vector \(\mathbf{x}\): \[ \mathbf{D}\mathbf{x}=-\mathbf{L}\mathbf{x}-\mathbf{U}\mathbf{x}+\mathbf{b}. \]

• Si escribimos la expresión anterior en función de la sucesión de aproximaciones \(\mathbf{x}^{(k)}\) obtenemos: \[ \mathbf{D}\mathbf{x}^{(k)}=-\mathbf{L}\mathbf{x}^{(k)}-\mathbf{U}\mathbf{x}^{(k-1)}+\mathbf{b}. \]
• “Pasamos” al término de la izquierda todo lo que depende de \(\mathbf{x}^{(k)}\): \[ (\mathbf{D}+\mathbf{L})\mathbf{x}^{(k)}=-\mathbf{U}\mathbf{x}^{(k-1)}+\mathbf{b}. \]
• Despejamos \(\mathbf{x}^{(k)}\): \[ \mathbf{x}^{(k)}=-(\mathbf{D}+\mathbf{L})^{-1}\mathbf{U}\mathbf{x}^{(k-1)}+(\mathbf{D}+\mathbf{L})^{-1}\mathbf{b}. \]

Entonces, la matriz de iteración y el vector independiente del método de Gauss-Seidel serán: \[ \mathbf{T}_{GS}=-(\mathbf{D}+\mathbf{L})^{-1}\mathbf{U},\quad \mathbf{c}_{GS}=(\mathbf{D}+\mathbf{L})^{-1}\mathbf{b}. \]

Vamos a dar el pseudocódigo del método de Gauss-Seidel usando como criterio de parada el criterio del error absoluto.

• INPUT matriz del sistema \(\mathbf{A}=(a_{ij})_{i=1,\ldots,n,j=1,\ldots,n}\), vector de términos independientes \(\mathbf{b}=(b_i)_{i=1,\ldots,n}\), valor inicial \(\mathbf{X0}=\mathbf{x}^{(0)}\), error absoluto o tolerancia TOL, número máximo de iteraciones Nmax.
• Set k=1.

• While k <= Nmax
  • For i=1,...,n
    • Set \(x_i=\frac{1}{a_{ii}}\left(-\sum_{j=1}^{i-1} a_{ij}x_j-\sum_{j=i+1}^n a_{ij}X0_j+b_i\right)\).
  • If \(\|\mathbf{x}-\mathbf{X0}\| <\) TOL then
    • Print \(x_1,\ldots,x_n\)
    • STOP.
  • Set k=k+1.
  • For i=1,...,n
    • Set \(X0_i =x_i\).
• Print número máximo de iteraciones alcanzado. El método no converge.
• STOP.

En esta sección vamos a estudiar un método iterativo general para resolver un sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\) donde la sucesión de aproximaciones \(\mathbf{x}^{(k)}\) viene dada en forma recurrente por la expresión siguiente: \[ \mathbf{x}^{(k)}=\mathbf{T}\mathbf{x}^{(k-1)}+\mathbf{c}, \] donde \(\mathbf{T}\) es la matriz de iteración y \(\mathbf{c}\), el vector independiente, usando un vector \(\mathbf{x}^{(0)}\) inicial.

Los métodos de Jacobi y Gauss-Seidel serían casos particulares de métodos iterativos definidos usando la expresión anterior.

Tal como comentamos en su momento, para hallar la recurrencia que gobierna la sucesión de aproximaciones \(\mathbf{x}^{(k)}\), tranformamos el sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\) en el sistema lineal equivalente \(\mathbf{x}=\mathbf{T}\mathbf{x}+\mathbf{c}\) y, a partir de la última expresión, definimos la recurrencia de las aproximaciones \(\mathbf{x}^{(k)}\).

Entonces, si la sucesión de aproximaciones \(\mathbf{x}^{(k)}\) converge a la solución \(\mathbf{x}\), ésta verificará: \[ \mathbf{x}=\mathbf{T}\mathbf{x}+\mathbf{c},\ \Rightarrow (\mathbf{Id}-\mathbf{T})\mathbf{x}=\mathbf{c},\ \Rightarrow \mathbf{x}=(\mathbf{Id}-\mathbf{T})^{-1}\mathbf{c}. \]

Veamos a continuación si somos capaces de escribir la sucesión de aproximaciones \(\mathbf{x}^{(k)}\) en función de la aproximación inicial \(\mathbf{x}^{(0)}\): \[ \begin{align*} \mathbf{x}^{(1)} & = \mathbf{T}\mathbf{x}^{(0)}+\mathbf{c},\\ \mathbf{x}^{(2)} & = \mathbf{T}\mathbf{x}^{(1)}+\mathbf{c}=\mathbf{T}(\mathbf{T}\mathbf{x}^{(0)}+\mathbf{c})+\mathbf{c}=\mathbf{T}^2\mathbf{x}^{(0)}+\mathbf{T}\mathbf{c}+\mathbf{c},\\ \mathbf{x}^{(3)} & = \mathbf{T}\mathbf{x}^{(2)}+\mathbf{c}=\mathbf{T}(\mathbf{T}^2\mathbf{x}^{(0)}+\mathbf{T}\mathbf{c}+\mathbf{c})+\mathbf{c}=\mathbf{T}^3\mathbf{x}^{(0)}+\mathbf{T}^2\mathbf{c}+\mathbf{T}{c}+\mathbf{c}. \end{align*} \] En general, puede demostrarse fácilmente por inducción que: \[ \mathbf{x}^{(k)}=\mathbf{T}^k\mathbf{x}^{(0)}+\mathbf{T}^{k-1}\mathbf{c}+\cdots +\mathbf{T}\mathbf{c}+\mathbf{c}=\mathbf{T}^k\mathbf{x}^{(0)}+\left(\sum_{i=0}^{k-1}\mathbf{T}^i\right)\mathbf{c}. \]

Teniendo en cuenta que el límite \(\mathbf{x}\) de la sucesión de aproximaciones, en caso de existir, verifica que \(\mathbf{x}=(\mathbf{Id}-\mathbf{T})^{-1}\mathbf{c}\), intuitivamente y usando la expresión anterior, \[ \mathbf{x}^{(k)}=\mathbf{T}^k\mathbf{x}^{(0)}+\left(\sum_{i=0}^{k-1}\mathbf{T}^i\right)\mathbf{c}, \] parece que para que la sucesión de aproximaciones sea convergente, la sucesión de matrices \(\mathbf{T}^k\) debe tender a la matriz \(\mathbf{0}\) y que el sumatorio anterior de matrices \(\displaystyle\left(\sum_{i=0}^{k-1}\mathbf{T}^i\right)\) debe tender a \((\mathbf{Id}-\mathbf{T})^{-1}\).

El siguiente lema es un primer intento de responder a las cuestiones planteadas:

Lema. Convergencia de una serie de matrices.

Sea \(\mathbf{T}\) una matriz cuadrada \(n\times n\) cuyo radio espectral es menor que \(1\), \(\rho(\mathbf{T})<1\). Entonces, \[ (\mathbf{Id}-\mathbf{T})^{-1}=\mathbf{Id}+\mathbf{T}+\mathbf{T}^2+\cdots =\sum_{j=0}^\infty \mathbf{T}^j. \] Es decir, la serie o la suma infinita de matrices \(\displaystyle\sum_{j=0}^\infty \mathbf{T}^j\) es convergente si el radio espectral de la matriz \(\mathbf{T}\) de la que hallamos las potencias es menor que \(1\).

Ya tenemos los ingredientes necesarios para demostrar el teorema que dice cuando un método iterativo es convergente:

Teorema. Condición necesaria y suficiente de convergencia.

Sea \(\mathbf{x}^{(k)}\) una sucesión de aproximaciones para hallar la solución del sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\) a partir de la recurrencia \(\mathbf{x}^{(k)}=\mathbf{T}\mathbf{x}^{(k-1)}+\mathbf{c}\). Entonces la sucesión anterior converge si, y sólo si, el radio espectral de la matriz de iteración \(\mathbf{T}\) es menor que \(1\), \(\rho(\mathbf{T})<1\).

El teorema anterior es importante porque nos da una condición necesaria y suficiente para que la sucesión de aproximaciones \(\mathbf{x}^{(k)}\) converja.

Sin embargo, no es aplicable en la práctica ya que hallar el radio espectral de una matriz es un problema mucho más complejo desde el punto de vista del análisis numèrico que resolver un sistema lineal.

El corolario siguiente intenta resolver este problema, es decir, nos da una condición suficiente de convergencia en función de la norma de la matriz del sistema:

Corolario.

Sea \(\mathbf{x}^{(k)}\) una sucesión de aproximaciones para hallar la solución del sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\) a partir de la recurrencia \(\mathbf{x}^{(k)}=\mathbf{T}\mathbf{x}^{(k-1)}+\mathbf{c}\).

Si existe una norma matricial para la cual, la norma de la matriz de iteración \(\mathbf{T}\) es menor que \(1\), \(\|\mathbf{T}\|<1\), entonces la sucesión de aproximaciones \(\mathbf{x}^{(k)}\) converge a la solución del sistema anterior \(\mathbf{x}\) para cualquier aproximación inicial \(\mathbf{x}^{(0)}\) y además los errores cometidos vienen acotados por las expresiones siguientes:

• \(\|\mathbf{x}^{(k)}-\mathbf{x}\|\leq \|\mathbf{T}\|^k\|\mathbf{x}^{(0)}-\mathbf{x}\|\),
• \(\|\mathbf{x}^{(k)}-\mathbf{x}\|\leq \frac{\|\mathbf{T}\|^k}{1-\|\mathbf{T}\|}\|\mathbf{x}^{(1)}-\mathbf{x}^{(0)}\|\).

Observación.

La primera desigualdad del teorema no tiene aplicación práctica al no conocer \(\mathbf{x}\) pero nos indica que cuanto menor sea la norma de la matriz de iteración, más “rápidamente” tenderá a cero el error cometido \(\|\mathbf{x}^{(k)}-\mathbf{x}\|\).

Sin embargo, usando la segunda desigualdad podemos estimar el número de iteraciones máximo a realizar para que el error \(\|\mathbf{x}^{(k)}-\mathbf{x}\|\) sea menor que una cierta tolerancia \(\epsilon\) ya que: \[ \|\mathbf{x}^{(k)}-\mathbf{x}\|\leq \frac{\|\mathbf{T}\|^k}{1-\|\mathbf{T}\|}\|\mathbf{x}^{(1)}-\mathbf{x}^{(0)}\|\leq \epsilon,\ \Rightarrow \|\mathbf{T}\|^k \leq \frac{\epsilon(1-\|\mathbf{T}\|)}{\|\mathbf{x}^{(1)}-\mathbf{x}^{(0)}\|}. \]

Observación (continuación).

Tomando logaritmos en la expresión anterior, tenemos que el número máximo de iteraciones a realizar para obtener una tolerancia menor que \(\epsilon\) vale: \[ k\leq \frac{\ln\left(\frac{\epsilon(1-\|\mathbf{T}\|)}{\|\mathbf{x}^{(1)}-\mathbf{x}^{(0)}\|}\right)}{\ln\|\mathbf{T}\|} = \frac{\ln\epsilon+\ln(1-\|\mathbf{T}\|)-\ln(\|\mathbf{x}^{(1)}-\mathbf{x}^{(0)}\|)}{\ln\|\mathbf{T}\|}. \]

Los métodos de Jacobi y Gauss-Seidel definían la sucesión de aproximaciones \(\mathbf{x}^{(k)}\) de la forma siguiente:

• Método de Jacobi: \[ \mathbf{x}^{(k)}=\mathbf{T}_J\mathbf{x}^{(k-1)}+\mathbf{c}_J. \]
• Método de Gauss-Seidel: \[ \mathbf{x}^{(k)}=\mathbf{T}_{GS}\mathbf{x}^{(k-1)}+\mathbf{c}_{GS}. \]

Las matrices de iteración eran las siguientes: \[ \mathbf{T}_J =-\mathbf{D}^{-1}(\mathbf{L}+\mathbf{U}),\quad\mathbf{T}_{GS} =-(\mathbf{D}+\mathbf{L})^{-1}\mathbf{U}, \] donde las matrices \(\mathbf{D}\), \(\mathbf{L}\) y \(\mathbf{U}\) eran la matriz diagonal, la matriz triangular inferior y la matriz triangular superior de la matriz del sistema a resolver \(\mathbf{A}\): \(\mathbf{A}=\mathbf{D}+\mathbf{L}+\mathbf{U}\).

Aplicando lo visto anteriormente, si el radio espectral de la matriz de iteración de Jacobi \(\rho(\mathbf{T}_J)<1\) es menor que \(1\), el método de Jacobi será convergente y lo mismo para el método de Gauss Seidel, es decir, si \(\rho(\mathbf{T}_{GS})<1\) es menor que \(1\), el método de Gauss-Seidel será convergente

De la misma manera, si existe una norma matricial para la que \(\|\mathbf{T}_J\|<1\), el método de Jacobi también será convergente y lo mismo para el método de Gauss-Seidel, es decir, si existe una norma matricial para la que \(\|\mathbf{T}_{GS}\|<1\), el método de Gauss-Seidel también será convergente y además en estos casos tenemos una manera de estimar el número de iteraciones para que el error entre \(\mathbf{x}^{(k)}\) y la solución exacta \(\mathbf{x}\) sea menor que una cierta tolerancia \(\epsilon\).

Veamos que si la matriz del sistema \(\mathbf{A}=(a_{ij})_{i,j=1,\ldots,n}\) es estrictamente diagonal dominante, donde recordemos que significa que \[ |a_{ii}| > \sum_{j=1,j\neq i}^n |a_{ij}|,\ i=1,\ldots,n, \] los métodos de Jacobi y Gauss-Seidel convergen:

Proposición. Convergencia de los métodos de Jacobi y Gauss-Seidel cuando la matriz del sistema es estrictamente diagonal dominante.

Sea \(\mathbf{x}^{(k)}_J\), \(\mathbf{x}^{(k)}_{GS}\) las sucesiones de aproximaciones generadas por los métodos de Jacobi y Gauss-Seidel, respectivamente, para hallar la solución del sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\). Si la matriz del sistema \(\mathbf{A}\) es estrictamente diagonal dominante, las sucesiones anteriores convergen a la solución de sistema \(\mathbf{x}\).

Para demostrar la proposición anterior veremos que \(\|\mathbf{T}_J\|_\infty <1\) y \(\|\mathbf{T}_{GS}\|_\infty <1\). Entonces, usando el corolario anterior, podemos afirmar que las sucesiones de aproximaciones generadas por los métodos de Jacobi y Gauss-Seidel convergen.

Hemos visto en la sección anterior que la clave para ver si la sucesión de aproximaciones \(\mathbf{x}^{(k)}\) es el radio espectral de la matriz de iteración \(\mathbf{T}\) de tal forma que la sucesión de aproximaciones se define como: \[ \mathbf{x}^{(k)}=\mathbf{T}\mathbf{x}^{(k-1)}+\mathbf{c}. \] Concretamente, si el radio espectral de la matriz de iteración es menor que \(1\), \(\rho(\mathbf{T})<1\), la sucesión de aproximaciones \(\mathbf{x}^{(k)}\) convergerá a la solución exacta del sistema lineal \(\mathbf{x}\).

La idea sería diseñar un método dependiendo de un parámetro y estudiar para qué valores del parámetro el radio espectral de la matriz de iteración correspondiente tiene un valor lo más pequeño posible.

Para ello, vamos a hallar un método iterativo dependiendo de un parámetro \(w\) hallando un sistema equivalente al sistema a resolver \(\mathbf{A}\mathbf{x}=\mathbf{b}\).

Recordemos que podemos descomponer la matriz del sistema \(\mathbf{A}\) en la parte diagonal \(\mathbf{D}\), la parte triangular inferior \(\mathbf{L}\) y la parte triangular superior \(\mathbf{U}\): \(\mathbf{A}=\mathbf{D}+\mathbf{L}+\mathbf{U}\).

Escribamos el sistema \(\mathbf{A}\mathbf{x}=\mathbf{b}\) en un sistema equivalente introduciendo un parámetro \(w\): \[ \begin{align*} \mathbf{A}\mathbf{x} & = \mathbf{b},\ \Rightarrow (\mathbf{D}+\mathbf{L}+\mathbf{U})\mathbf{x} = \mathbf{b},\\ w(\mathbf{D}+\mathbf{L}+\mathbf{U})\mathbf{x} & = w\mathbf{b},\\ (w\mathbf{D}+w\mathbf{L}+w\mathbf{U})\mathbf{x} & = w\mathbf{b},\\ (\mathbf{D}+(w-1)\mathbf{D}+w\mathbf{L}+w\mathbf{U})\mathbf{x} & = w\mathbf{b},\\ (\mathbf{D}+w\mathbf{L})\mathbf{x} & = ((1-w)\mathbf{D}-w\mathbf{U})\mathbf{x}+w\mathbf{b}. \end{align*} \]

La sucesión de aproximaciones \(\mathbf{x}^{(k)}_w\) definida por los métodos anteriores parametrizados por el parámetro \(w\) viene dada por: \[ (\mathbf{D}+w\mathbf{L})\mathbf{x}^{(k)}_w = ((1-w)\mathbf{D}-w\mathbf{U})\mathbf{x}^{(k-1)}_w +w\mathbf{b} \]

La matriz de iteración del método anterior parametrizado por el parámetro \(w\) será: \(\mathbf{T}_w =(\mathbf{D}+w\mathbf{L})^{-1}((1-w)\mathbf{D}-w\mathbf{U})\) con vector independiente \(\mathbf{c}_w =w(\mathbf{D}+w\mathbf{L})^{-1}\mathbf{b}\).

Fijarse que para \(w=1\), redescubrimos el conocido método de Gauss-Seidel. Por tanto, los métodos anteriores pueden considerarse una generalización del método de Gauss-Seidel.

Los métodos anteriores se denominan métodos de relajación. Si \(0<w<1\), el método se denomina método de infrarelajación y si \(w>1\), método de sobrerelajación. Dichos métodos se abrevian con métodos SOR. (successive over-relaxation)

Vamos a dar el pseudocódigo del método SOR usando como criterio de parada el criterio del error absoluto.

• INPUT matriz del sistema \(\mathbf{A}=(a_{ij})_{i=1,\ldots,n,j=1,\ldots,n}\), vector de términos independientes \(\mathbf{b}=(b_i)_{i=1,\ldots,n}\), valor inicial \(\mathbf{X0}=\mathbf{x}^{(0)}\), error absoluto o tolerancia TOL, número máximo de iteraciones Nmax, valor del parámetro \(w\).
• Set k=1.

• While k <= Nmax
  • For i=1,...,n
    • Set \(x_i=(1-w)X0_i+\frac{w}{a_{ii}}\left(-\sum_{j=1}^{i-1} a_{ij}x_j-\sum_{j=i+1}^n a_{ij}X0_j+b_i\right)\).
  • If \(\|\mathbf{x}-\mathbf{X0}\| <\) TOL then
    • Print \(x_1,\ldots,x_n\)
    • STOP.
  • Set k=k+1.
  • For i=1,...,n
    • Set \(X0_i =x_i\).
• Print número máximo de iteraciones alcanzado. El método no converge.
• STOP.

El siguiente resultado nos dice que para que un método iterativo SOR converja, el parámetro tiene que estar necesariamente entre \(0\) y \(2\):

Teorema de Kahan.

Sea \(\mathbf{x}^{(k)}_w\) una sucesión de aproximaciones para hallar la solución del sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\) a partir de la recurrencia \(\mathbf{x}^{(k)}_w=\mathbf{T}_w\mathbf{x}^{(k-1)}_w+\mathbf{c}_w\) dada por un método SOR.

Supongamos que \(a_{ii}\neq 0\), para todo \(i=1,2,\ldots,n\). Entonces el radio espectral de la matriz de iteración \(\mathbf{T}_w\) cumple: \(\rho(\mathbf{T}_w)\geq |w-1|.\)

Corolario.

Sea \(\mathbf{x}^{(k)}_w\) una sucesión de aproximaciones para hallar la solución del sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\) a partir de la recurrencia \(\mathbf{x}^{(k)}_w=\mathbf{T}_w\mathbf{x}^{(k-1)}_w+\mathbf{c}_w\) dada por un método SOR.

Si \(a_{ii}\neq 0\), \(i=1,\ldots,n\) y si la sucesión de aproximaciones anterior converge, el valor del parámetro \(w\) está entre \(0\) y \(2\): \(0<w<2\).

Tal como vimos en el ejemplo anterior, la convergencia de un método SOR no está asegurada para cualquier valor \(w\in (0,2)\).

Sin embargo, existe el resultado siguiente:

Teorema de Ostrowski-Reich.

Si la matriz del sistema \(\mathbf{A}\) es simétrica y definida positiva y \(0<w<2\), el método iterativo SOR converge para cualquier aproximación inicial \(\mathbf{x}^{(0)}\).

El resultado siguiente nos dice cuál es el valor de \(w\) óptimo, es decir, el valor que minimiza el radio espectral de la matriz de iteración \(\mathbf{T}_w\) en el caso en que la matriz del sistema \(\mathbf{A}\) cumpla las condiciones siguientes:

Teorema.

Si la matriz \(\mathbf{A}\) del sistema es simétrica, definida positiva y tridiagonal, el radio espectral de la matriz de iteración del método de Gauss-Seidel es el cuadrado del radio espectral de la matriz de iteración del método de Jacobi, \(\rho(\mathbf{T}_{GS})=\rho(\mathbf{T}_J)^2 <1\) y además dichos métodos son convergentes.

El valor del parámetro \(w\) óptimo es el siguiente: \(\tilde{w}=\frac{2}{1+\sqrt{1-\rho(\mathbf{T}_J)^2}}=\frac{2}{\rho(\mathbf{T}_J)^2}\left(1-\sqrt{1-\rho(\mathbf{T}_J)^2}\right),\) con valor de radio espectral: \(\rho(\mathbf{T}_{\tilde{w}})=\tilde{w}-1=\frac{2}{\rho(\mathbf{T}_J)^2}\left(2-\rho(\mathbf{T}_J)^2-2\sqrt{1-\rho(\mathbf{T}_J)^2}\right).\)

• Cuando la dimensión de la matriz del sistema lineal a resolver \(\mathbf{A}\mathbf{x}=\mathbf{b}\) no es demasiado grande, es mejor usar los métodos directos al ser éstos más eficientes en este caso y proporcionar soluciones exactas.

• En cambio, si la dimensión de la matriz del sistema lineal es muy grande y es una matriz sparse que significa que tiene muchos ceros sin ninguna estructura predefinida, es mejor usar métodos iterativos ya que ahorran en almacenamiento y tienen un coste computacional más bajo que los métodos directos. Los sistemas de este segundo tipo aparecen en problemas de frontera en ecuaciones en derivadas parciales así como en problemas de redes neuronales en el campo de la inteligencia artificial.

El método del gradiente conjugado, como veremos, puede usarse como método directo donde resuelve el sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\) en \(n\) pasos o como método iterativo hallando soluciones aproximadas del sistema anterior.

Como método directo, es computacionalmente más costoso que el método de Gauss con pivotaje.

Sin embargo, como método iterativo, es muy útil para resolver sistemas sparse, es decir, sistemas lineales donde la matriz del sistema \(\mathbf{A}\) tiene dimensiones \(n\times n\) muy grandes y con muchos ceros. Este tipo de sistemas aparecen en el análisis de redes neuronales.

Recordemos que dentro del ámbito del álgebra lineal, suponemos que los vectores vienen representados en columnas. En este contexto, \(\mathbf{u}^\top \mathbf{v}\) simboliza el producto escalar entre los vectores \(\mathbf{u}\) y \(\mathbf{v}\).

El método del gradiente conjugado está basado en el resultado siguiente:

Teorema.

Consideremos el sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\), donde la matriz del sistema \(\mathbf{A}\) es simétrica y definida positiva.

Entonces, el vector \(\tilde{\mathbf{x}}\) es solución del sistema anterior \(\mathbf{A}\mathbf{x}=\mathbf{b}\) si, y sólo si, \(\tilde{\mathbf{x}}\) minimiza la función siguiente \(g(\mathbf{x})=\mathbf{x}^\top (\mathbf{A}\mathbf{x})-2(\mathbf{x}^\top \mathbf{b})\).

Para desarrollar el algoritmo del gradiente conjugado, supongamos que tenemos \(\mathbf{x}\) una solución aproximada del sistema \(\mathbf{A}\mathbf{x}=\mathbf{b}\) y consideramos un vector \(\mathbf{v}\neq\mathbf{0}\).

La idea es buscar una dirección para avanzar a partir de la solución aproximada anterior \(\mathbf{x}\) que nos dé una aproximación mejor de tal forma que minimice la función \(g(\mathbf{x})\) definida anteriormente ya que el teorema anterior dice que cuanto mejor minimicemos la función \(g(\mathbf{x})\), mejor solución hallaremos del sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\).

Sea \(\mathbf{r}=\mathbf{b}-\mathbf{A}\mathbf{x}\) el residuo asociado a la solución aproximada \(\mathbf{x}\) y sea \(t=\frac{\mathbf{v}^\top \mathbf{r}}{\mathbf{v}^\top \mathbf{A}\mathbf{v}},\) el valor que hace que la función \(g(\mathbf{x}+t\mathbf{v})\) sea mínimo en la dirección \(\mathbf{v}\), tal como vimos en la demostración del teorema.

Definimos \(\hat{\mathbf{x}}=\mathbf{x}+t\mathbf{v}\) la nueva aproximación a partir de la dirección \(\mathbf{v}\) que sabemos que es la mejor ya que \(\hat{\mathbf{x}}\) minimiza \(g(\mathbf{x})\) en la dirección \(\mathbf{v}\), es decir, minimiza la función \(h(t)=g(\mathbf{x}+t\mathbf{v})\).

A partir de las consideraciones anteriores, ya tenemos el algoritmo preliminar del gradiente conjugado:

• Sea \(\mathbf{x}^{(0)}\) una aproximación inicial y sea \(\mathbf{v}^{(1)}\neq\mathbf{0}\) una dirección de búsqueda inicial. Sea \(\mathbf{r}^{(0)}=\mathbf{b}-\mathbf{A}\mathbf{x}^{(0)}\).
• Calculamos \(\mathbf{x}^{(1)}\) de la siguiente forma: \[ t_1=\frac{(\mathbf{v}^{(1)})^\top (\mathbf{r}^{(0)})}{(\mathbf{v}^{(1)})^\top (\mathbf{A}\mathbf{v}^{(1)})},\quad \mathbf{x}^{(1)}=\mathbf{x}^{(0)}+t_1\mathbf{v}^{(1)}. \]
• En general, supongamos que tenemos \(\mathbf{x}^{(0)},\ldots,\mathbf{x}^{(k-1)}\) y direcciones de búsqueda \(\mathbf{v}^{(1)},\ldots,\mathbf{v}^{(k)}\). Sea \(\mathbf{r}^{(k-1)}=\mathbf{b}-\mathbf{A}\mathbf{x}^{(k-1)}\). Definimos: \[ t_k=\frac{(\mathbf{v}^{(k)})^\top (\mathbf{r}^{(k-1)})}{(\mathbf{v}^{(k)})^\top (\mathbf{A}\mathbf{v}^{(k)})},\quad \mathbf{x}^{(k)}=\mathbf{x}^{(k-1)}+t_k\mathbf{v}^{(k)}. \]

En la descripción del algoritmo anterior, no se indica cómo encontrar las direcciones de búsqueda. Éste es el motivo que sea una versión preliminar.

Dada una aproximación inicial, la dirección de búsqueda tiene que ser aquélla que tenga una máxima disminución en la función \(g(\mathbf{x})\).

Dada una función de \(n\) variables \(g(\mathbf{x})\) y una aproximación inicial \(\mathbf{x}\) sabemos a partir del Análisis multivariante que la dirección de disminución máxima de la función \(g(\mathbf{x})\) viene dada por el gradiente de \(g\): \[ \nabla g(\mathbf{x})=\left(\frac{\partial g}{\partial x_1}(\mathbf{x}),\frac{\partial g}{\partial x_2}(\mathbf{x}),\ldots,\frac{\partial g}{\partial x_n}(\mathbf{x})\right). \]

Hallemos pues el gradiente de la función \(g(\mathbf{x})\).

Recordemos que la función \(g(\mathbf{x})=g(x_1,\ldots,x_n)\) era la siguiente: \[ g(\mathbf{x})=\mathbf{x}^\top (\mathbf{A}\mathbf{x})-2(\mathbf{x}^\top \mathbf{b})=\sum_{i,j=1}^n a_{ij}x_i x_j-2\sum_{i=1}^n x_i b_i. \] El valor de \(\frac{\partial g}{\partial x_k}(\mathbf{x})\) será: \(\displaystyle\frac{\partial g}{\partial x_k}(\mathbf{x})=2\sum_{i=1}^n a_{ki}x_i-2 b_k.\)

Matricialmente: \[ \nabla g(\mathbf{x})=2\mathbf{A}\mathbf{x}-2\mathbf{b}=2(\mathbf{A}\mathbf{x}-\mathbf{b})=-2\mathbf{r}, \] donde recordemos que \(\mathbf{r}\) es el residuo: \(\mathbf{r}=\mathbf{b}-\mathbf{A}\mathbf{x}\).

Entonces en el algoritmo anterior, la dirección de búsqueda en el paso \(k\)-ésimo será: \(\mathbf{v}^{(k)}=\mathbf{r}^{(k-1)}=\mathbf{b}-\mathbf{A}\mathbf{x}^{(k-1)}\).

El método que usa la dirección de búsqueda anterior se denomina método del máximo descenso. Dicho método es raramente usado en la resolución de sistemas lineales debido a la lentitud de la convergencia del mismo. Sin embargo, sí se usa en problemas de resolución de sistemas no lineales y problemas de optimización.

Una manera de acelerar la convergencia es usar un conjunto de direcciones de búsqueda \(\mathbf{v}^{(1)},\ldots, \mathbf{v}^{(k)}\) que sean \(\mathbf{A}\)-ortogonales, es decir, \((\mathbf{v}^{(i)})^\top (\mathbf{A}\mathbf{v}^{(j)})=0\), si \(i\neq j\).

Antes de describir cómo hallar las direcciones de búsqueda anteriores, enunciemos el teorema que nos asegura que si las direcciones de búsqueda son \(\mathbf{A}\)-ortogonales, resolvemos el sistema lineal en \(n\) pasos:

Teorema.

Sean \(\mathbf{v}^{(1)},\ldots,\mathbf{v}^{(n)}\), \(n\) vectores no nulos \(\mathbf{A}\)-ortogonales. Sea \(\mathbf{x}^{(0)}\) un vector inicial arbitrario. Definimos, \[ t_k =\frac{(\mathbf{v}^{(k)})^\top (b-\mathbf{A}\mathbf{x}^{(k-1)})}{(\mathbf{v}^{(k)})^\top (\mathbf{A}\mathbf{v}^{(k)})},\quad \mathbf{x}^{(k)}=\mathbf{x}^{(k-1)}+t_k\mathbf{v}^{(k)}, \] para \(k=1,\ldots,n\). Entonces, suponiendo que no hay errores de redondeo ni errores en las operaciones, el vector \(\mathbf{x}^{(n)}\) es la solución del sistema lineal \(\mathbf{A}\mathbf{x}=\mathbf{b}\), es decir, \(\mathbf{A}\mathbf{x}^{(n)}=\mathbf{b}\).

Antes de demostrar el teorema, demostremos el lema siguiente:

Lema.

Sean \(\mathbf{v}^{(1)},\ldots,\mathbf{v}^{(n)}\), \(n\) vectores no nulos \(\mathbf{A}\)-ortogonales, donde la matriz \(\mathbf{A}\) es simétrica, definida positiva y no singular.. Entonces los vectores anteriores \(\mathbf{v}^{(1)},\ldots,\mathbf{v}^{(n)}\) son linealmente independientes.

Veamos a continuación la demostración del teorema.

Antes de enunciar el algoritmo veamos que si los vectores \(\mathbf{v}^{(j)}\) son \(\mathbf{A}\)-ortogonales, entonces los residuos \(\mathbf{r}^{(k)}=\mathbf{b}-\mathbf{A}\mathbf{x}^{(k)}\) y los vectores \(\mathbf{v}^{(j)}\) son ortogonales:

Proposición.

Sean \(\mathbf{v}^{(1)},\ldots,\mathbf{v}^{(n)}\), \(n\) vectores no nulos \(\mathbf{A}\)-ortogonales. Sea \(\mathbf{x}^{(0)}\) un vector inicial arbitrario y sea \(\mathbf{r}^{(0)}=\mathbf{b}-\mathbf{A}\mathbf{x}^{(0)}\). Definimos: \[ t_k =\frac{(\mathbf{v}^{(k)})^\top (b-\mathbf{A}\mathbf{x}^{(k-1)})}{(\mathbf{v}^{(k)})^\top (\mathbf{A}\mathbf{v}^{(k)})},\quad \mathbf{x}^{(k)}=\mathbf{x}^{(k-1)}+t_k\mathbf{v}^{(k)},\quad \mathbf{r}^{(k)}=\mathbf{b}-\mathbf{A}\mathbf{x}^{(k)}, \] para \(k=1,\ldots,n\). Entonces los vectores \(\mathbf{v}^{(j)}\), \(j=1,\ldots,k\) son ortogonales al vector \(\mathbf{r}^{(k)}\): \((\mathbf{r}^{(k)})^\top \mathbf{v}^{(j)}=0\), \(j=1,\ldots,k\).

El algoritmo del gradiente conjugado se basa en definir la sucesión de soluciones aproximadas \(\mathbf{x}^{(k)}\) a partir de una solución inicial \(\mathbf{x}^{(0)}\) usando la recurrencia vista anteriormente: \[ t_k =\frac{(\mathbf{v}^{(k)})^\top (b-\mathbf{A}\mathbf{x}^{(k-1)})}{(\mathbf{v}^{(k)})^\top (\mathbf{A}\mathbf{v}^{(k)})}=\frac{(\mathbf{v}^{(k)})^\top \mathbf{r}^{(k-1)}}{(\mathbf{v}^{(k)})^\top (\mathbf{A}\mathbf{v}^{(k)})},\quad \mathbf{x}^{(k)}=\mathbf{x}^{(k-1)}+t_k\mathbf{v}^{(k)}, \] eligiendo los vectores \(\mathbf{v}^{(k)}\), \(\mathbf{A}\)-ortogonales y además se verificará que los residuos \(\mathbf{r}^{(k)}=\mathbf{b}-\mathbf{A}\mathbf{x}^{(k)}\) serán ortogonales entre sí: \({\mathbf{r}^{(i)}}^\top \mathbf{r}^{(j)}=0\), para \(i\neq j\).

Veamos cómo funciona.

Sea \(\mathbf{x}^{(0)}\) una aproximación inicial con residuo inicial \(\mathbf{r}^{(0)}=\mathbf{b}-\mathbf{A}\mathbf{x}^{(0)}\).

• Elegimos \(\mathbf{v}^{(1)}=\mathbf{r}^{(0)}\). La aproximación \(\mathbf{x}^{(1)}\) es de la forma \(\mathbf{x}^{(1)}=\mathbf{x}^{(0)}+t_1\mathbf{v}^{(1)}\), donde \(t_1 =\frac{(\mathbf{v}^{(1)})^\top (b-\mathbf{A}\mathbf{x}^{(0)})}{(\mathbf{v}^{(1)})^\top (\mathbf{A}\mathbf{v}^{(1)})}=\frac{(\mathbf{v}^{(1)})^\top \mathbf{r}^{(0)}}{(\mathbf{v}^{(1)})^\top (\mathbf{A}\mathbf{v}^{(1)})}\).

• Seguidamente tenemos que hallar \(\mathbf{x}^{(2)}=\mathbf{x}^{(1)}+t_2\mathbf{v}^{(2)}\).

  • Escribimos \(\mathbf{v}^{(2)}\) de la forma: \(\mathbf{v}^{(2)}=\mathbf{r}^{(1)}+s_1\mathbf{v}^{(1)}\).
  • Como \(\mathbf{v}^{(1)}\) y \(\mathbf{v}^{(2)}\) son \(\mathbf{A}\)-ortogonales, \[ {\mathbf{v}^{(1)}}^\top \mathbf{A}\mathbf{v}^{(2)}=0,\ \Rightarrow {\mathbf{v}^{(1)}}^\top \mathbf{A}\mathbf{r}^{(1)}+s_1{\mathbf{v}^{(1)}}^\top \mathbf{A}\mathbf{v}^{(1)}=0, \] de donde deducimos que \(s_1=-\frac{{\mathbf{v}^{(1)}}^\top \mathbf{A}\mathbf{r}^{(1)}}{{\mathbf{v}^{(1)}}^\top \mathbf{A}\mathbf{v}^{(1)}}.\)

• • El siguiente paso es hallar \(t_2\): \[ \begin{align*} t_2 & =\frac{(\mathbf{v}^{(2)})^\top \mathbf{r}^{(1)}}{(\mathbf{v}^{(2)})^\top (\mathbf{A}\mathbf{v}^{(2)})}=\frac{(\mathbf{r}^{(1)}+s_1\mathbf{v}^{(1)})^\top \mathbf{r}^{(1)}}{(\mathbf{v}^{(2)})^\top (\mathbf{A}\mathbf{v}^{(2)})} ,\\ & = \frac{(\mathbf{r}^{(1)})^\top \mathbf{r}^{(1)}}{(\mathbf{v}^{(2)})^\top (\mathbf{A}\mathbf{v}^{(2)})}+s_1 \frac{(\mathbf{v}^{(1)})^\top \mathbf{r}^{(1)}}{(\mathbf{v}^{(2)})^\top (\mathbf{A}\mathbf{v}^{(2)})} \end{align*} \]
  • Por último tenemos la aproximación \(\mathbf{x}^{(2)}=\mathbf{x}^{(1)}+t_2\mathbf{v}^{(2)}\).

• En general, supongamos que hemos hallado las aproximaciones \(\mathbf{x}^{(0)},\mathbf{x}^{(1)},\ldots,\mathbf{x}^{(k-1)}\) con direcciones conjugadas \(\mathbf{v}^{(1)},\ldots,\mathbf{v}^{(k-1)}\) que suponemos \(\mathbf{A}\)-ortogonales, es decir \((\mathbf{v}^{(i)})^\top \mathbf{A}\mathbf{v}^{(j)}=0\), para \(i,j=1,\ldots,k-1\) con \(i\neq j\).
  • La nueva aproximación \(\mathbf{x}^{(k)}\) se escribe de la forma: \(\mathbf{x}^{(k)}=\mathbf{x}^{(k-1)}+t_k\mathbf{v}^{(k)}\).

• • Escribimos la nueva dirección conjugada \(\mathbf{v}^{(k)}\) como \(\mathbf{v}^{(k)}=\mathbf{r}^{(k-1)}+s_{k-1}\mathbf{v}^{(k-1)}\). Como \(\mathbf{v}^{(k)}\) y \(\mathbf{v}^{(k-1)}\) son \(\mathbf{A}\)-ortogonales, tenemos que: \[ {\mathbf{v}^{(k-1)}}^\top \mathbf{A}\mathbf{v}^{(k)}=0,\ \Rightarrow {\mathbf{v}^{(k-1)}}^\top \mathbf{A}\mathbf{r}^{(k-1)}+s_{k-1}{\mathbf{v}^{(k-1)}}^\top \mathbf{A}\mathbf{v}^{(k-1)}=0, \] de donde deducimos que \(s_{k-1}=-\frac{{\mathbf{v}^{(k-1)}}^\top \mathbf{A}\mathbf{r}^{(k-1)}}{{\mathbf{v}^{(k-1)}}^\top \mathbf{A}\mathbf{v}^{(k-1)}}.\)

Se puede demostrar usando técnicas avanzadas de álgebra lineal que \((\mathbf{v}^{(i)})^\top \mathbf{A}\mathbf{v}^{(k)}=0\), para \(i=1,\ldots,k-2\). Por tanto, los vectores \(\mathbf{v}^{(1)},\ldots,\mathbf{v}^{(k)}\) serán \(A\)-ortogonales.

• • El siguiente paso es hallar \(t_k\): \[ \begin{align*} t_k & =\frac{(\mathbf{v}^{(k)})^\top \mathbf{r}^{(k-1)}}{(\mathbf{v}^{(k)})^\top (\mathbf{A}\mathbf{v}^{(k)})}=\frac{(\mathbf{r}^{(k-1)}+s_{k-1}\mathbf{v}^{(k-1)})^\top \mathbf{r}^{(k-1)}}{(\mathbf{v}^{(k)})^\top (\mathbf{A}\mathbf{v}^{(k)})} ,\\ & = \frac{(\mathbf{r}^{(k-1)})^\top \mathbf{r}^{(k-1)}}{(\mathbf{v}^{(k)})^\top (\mathbf{A}\mathbf{v}^{(k)})}+s_{k-1} \frac{(\mathbf{v}^{(k-1)})^\top \mathbf{r}^{(k-1)}}{(\mathbf{v}^{(k)})^\top (\mathbf{A}\mathbf{v}^{(k)})}=\frac{(\mathbf{r}^{(k-1)})^\top \mathbf{r}^{(k-1)}}{(\mathbf{v}^{(k)})^\top (\mathbf{A}\mathbf{v}^{(k)})}. \end{align*} \] En la última igualdad hemos usado la proposición anterior que asegura que \((\mathbf{r}^{(k-1)})^\top \mathbf{v}^{(k-1)}=0\).
  • Por último tenemos la aproximación \(\mathbf{x}^{(k)}=\mathbf{x}^{(k-1)}+t_k\mathbf{v}^{(k)}\).

Observación.

Los residuos \(\mathbf{r}^{(k)}\) cumplen la expresión siguiente: \[ \begin{align*} \mathbf{r}^{(k)} & =\mathbf{b}-\mathbf{A}\mathbf{x}^{(k)}=\mathbf{b}-\mathbf{A}(\mathbf{x}^{(k-1)}+t_k\mathbf{v}^{(k)})=\mathbf{b}-\mathbf{A}\mathbf{x}^{(k-1)}-t_k\mathbf{A}\mathbf{v}^{(k)}\\ & =\mathbf{r}^{(k-1)}-t_k\mathbf{A}\mathbf{v}^{(k)} \end{align*} \] A partir de la expresión anterior, tenemos que: \[ (\mathbf{r}^{(k)})^\top \mathbf{r}^{(k)}=(\mathbf{r}^{(k)})^\top \mathbf{r}^{(k-1)}-t_k(\mathbf{r}^{(k)})^\top (\mathbf{A}\mathbf{v}^{(k)})=-t_k (\mathbf{r}^{(k)})^\top \mathbf{A}\mathbf{v}^{(k)}. \] Además, a partir de la definición de \(t_k\), tenemos que: \[ (\mathbf{v}^{(k)})^\top \mathbf{A}{\mathbf{v}^{(k)}}=\frac{1}{t_{k}}(\mathbf{v}^{(k)})^\top \mathbf{r}^{(k-1)}. \]

Observación (continuación).

Por último, podemos simplificar el valor de \(s_{k}\) de la forma siguiente: \[ \begin{align*} s_{k} & =-\frac{(\mathbf{v}^{(k)})^\top \mathbf{A}\mathbf{r}^{(k)}}{{\mathbf{v}^{(k)}}^\top \mathbf{A}\mathbf{v}^{(k)}}=-\frac{(\mathbf{r}^{(k)})^\top \mathbf{A}\mathbf{v}^{(k)}}{{\mathbf{v}^{(k)}}^\top \mathbf{A}\mathbf{v}^{(k)}}= \frac{\frac{1}{t_{k}}(\mathbf{r}^{(k)})^\top \mathbf{r}^{(k)}}{\frac{1}{t_{k}}(\mathbf{v}^{(k)})^\top \mathbf{r}^{(k-1)}}\\ &= \frac{(\mathbf{r}^{(k)})^\top \mathbf{r}^{(k)}}{(\mathbf{r}^{(k-1)})^\top \mathbf{v}^{(k)}}= \frac{(\mathbf{r}^{(k)})^\top \mathbf{r}^{(k)}}{(\mathbf{r}^{(k-1)})^\top (\mathbf{r}^{(k-1)}+s_{k-1}\mathbf{v}^{(k-1)})}\\ & =\frac{(\mathbf{r}^{(k)})^\top \mathbf{r}^{(k)}}{(\mathbf{r}^{(k-1)})^\top \mathbf{r}^{(k-1)}+s_{k-1}(\mathbf{r}^{(k-1)})^\top \mathbf{v}^{(k-1)}}=\frac{(\mathbf{r}^{(k)})^\top \mathbf{r}^{(k)}}{(\mathbf{r}^{(k-1)})^\top \mathbf{r}^{(k-1)}}. \end{align*} \] En la última igualdad hemos usado la proposición anterior que asegura que \((\mathbf{r}^{(k-1)})^\top \mathbf{v}^{(k-1)}=0\).

En resumen, el algoritmo del gradiente conjugado se compone de los pasos siguientes:

• Paso \(0\): empezamos con una aproximación inicial \(\mathbf{x}^{(0)}\) y un residuo inicial \(\mathbf{r}^{(0)}=\mathbf{b}-\mathbf{A}\mathbf{x}^{(0)}\).
• Paso \(1\): consideramos como dirección inicial \(\mathbf{v}^{(1)}\), \(\mathbf{v}^{(1)}=\mathbf{r}^{(0)}\).
  • Calculamos \(\mathbf{x}^{(1)}=\mathbf{x}^{(0)}+t_1\mathbf{v}^{(1)}\), con \(t_1=\frac{(\mathbf{r}^{(0)})^\top \mathbf{r}^{(0)}}{(\mathbf{v}^{(1)})^\top (\mathbf{A}\mathbf{v}^{(1)})}.\)
• \(\cdots\)

• Paso \(k\): suponemos calculadas las aproximaciones \(\mathbf{x}^{(0)},\ldots,\mathbf{x}^{(k-1)}\), los residuos \(\mathbf{r}^{(0)},\ldots,\mathbf{r}^{(k-1)}\) y las direcciones \(\mathbf{v}^{(1)},\ldots,\mathbf{v}^{(k-1)}\).
  • Calculamos \(\mathbf{v}^{(k)}\): \(\mathbf{v}^{(k)}=\mathbf{r}^{(k-1)}+s_{k-1}\mathbf{v}^{(k-1)}\), con \(s_{k-1}=\frac{(\mathbf{r}^{(k-1)})^\top \mathbf{r}^{(k-1)}}{(\mathbf{r}^{(k-2)})^\top \mathbf{r}^{(k-2)}}.\)
  • Calculamos \(\mathbf{x}^{(k)}=\mathbf{x}^{(k-1)}+t_k\mathbf{v}^{(k)}\), con \(t_k=\frac{(\mathbf{r}^{(k-1)})^\top \mathbf{r}^{(k-1)}}{(\mathbf{v}^{(k)})^\top (\mathbf{A}\mathbf{v}^{(k)})}.\)

hasta llegar al paso \(k=n\) donde \(\mathbf{x}^{(n)}\) es la solución del sistema \(\mathbf{A}\mathbf{x}=\mathbf{b}\).

• INPUT matriz del sistema \(\mathbf{A}=(a_{ij})_{i=1,\ldots,n,j=1,\ldots,n}\), vector de términos independientes \(\mathbf{b}=(b_i)_{i=1,\ldots,n}\), valor inicial \(\mathbf{x}=\mathbf{x}^{(0)}\).

• Set \(\mathbf{r}=\mathbf{b}-\mathbf{A}\mathbf{x}\). (Calculamos \(\mathbf{r}^{(0)}\))

• Set \(\mathbf{v}=\mathbf{r}\). (Calculamos \(\mathbf{v}^{(1)}\))

• Set \(t=\frac{\mathbf{r}^\top \mathbf{r}}{\mathbf{v}^\top \mathbf{A}\mathbf{v}}\). (Calculamos \(t_1\))

• Set \(\mathbf{x}=\mathbf{x}+t\mathbf{v}\). (Calculamos \(\mathbf{x}^{(1)}\))

• Set \(\mathbf{r}_1=\mathbf{b}-\mathbf{A}\mathbf{x}\). (Calculamos \(\mathbf{r}^{(1)}\))

• For k=2,...,n
  • Set \(s=\frac{(\mathbf{r}_{1})^\top \mathbf{r}_{1}}{\mathbf{r}^\top \mathbf{r}}\). (Calculamos \(s_{k-1}\))
  • Set \(\mathbf{v}=\mathbf{r}_1+s\mathbf{v}\). (Calculamos \(\mathbf{v}^{(k)}\))
  • Set \(t=\frac{(\mathbf{r}_{1})^\top \mathbf{r}_{1}}{\mathbf{v}^\top \mathbf{A}\mathbf{v}}\). (Calculamos \(t_k\))
  • Set \(\mathbf{x}=\mathbf{x}+t\mathbf{v}\). (Calculamos \(\mathbf{x}^{(k)}\))
  • Set \(\mathbf{r}=\mathbf{r}_1\). (Actualizamos \(\mathbf{r}^{(k-1)}\) para el cálculo de \(s_{k-1}\) en la próxima iteración)
  • Set \(\mathbf{r}_{1}=\mathbf{b}-\mathbf{A}\mathbf{x}\). (Calculamos \(\mathbf{r}^{(k)}\))
• Print \(\mathbf{x}\) (Solución del sistema \(\mathbf{A}\mathbf{x}=\mathbf{b}\))