Sobrerrelajación sucesiva

En álgebra lineal numérica, el método de sobre-relajación sucesiva (SOR), es una variante del método de Gauss-Seidel para estimar la solución de un sistema lineal de ecuaciones, permitiendo una convergencia más rápida.

Fue propuesto simultáneamente por David M. Young Jr. y Stanley P. Frankel en 1950, con el propósito de resolver sistemas lineales en ordenadores digitales. Anteriormente ya existían métodos de tipo sobre-relajación, como el método de Lewis Fry Richardson, y los métodos desarrollados por R. V. Southwell. Sin embargo, estos últimos estaban diseñados para ser utilizados por calculadoras humanas, requiriendo alguna pericia para asegurar convergencia a la solución, lo que los hacía inaplicables para ordenadores digitales. Se pueden encontrar detalles de estos aspectos en la tesis de David M. Young Jr.^[1]

Formulación matricial

Se considera un sistema lineal cuadrado, de $n$ ecuaciones, con variable desconocida $\mathbf {x}$ :

A\mathbf {x} =\mathbf {b} ,\quad A={\begin{bmatrix}a_{11}&a_{12}&\cdots &a_{1n}\\a_{21}&a_{22}&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &a_{nn}\end{bmatrix}},\qquad \mathbf {x} ={\begin{bmatrix}x_{1}\\x_{2}\\\vdots \\x_{n}\end{bmatrix}},\qquad \mathbf {b} ={\begin{bmatrix}b_{1}\\b_{2}\\\vdots \\b_{n}\end{bmatrix}}.

La matriz A se puede escribir como la suma de: su componente diagonal $D$ , y sus componentes estrictamente triangular inferior y superior, $L$ y $U$ respectivamente: $A=D+L+U;$ dónde:

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

De esta forma, el sistema de ecuaciones lineales puede ser escrito como: $(D+\omega L)\mathbf {x} =\omega \mathbf {b} -[\omega U+(\omega -1)D]\mathbf {x} ;$ para cualquier constante $\omega >1$ , denominada factor de relajación. El método SOR es una técnica iterativa que en cada iteración "despeja" $x$ del lado izquierdo de esta igualdad, utilizando el valor de $x$ del paso anterior en el lado derecho. Analíticamente, esto puede ser escrito como: $\mathbf {x} ^{(k+1)}=(D+\omega L)^{-1}{\big (}-[\omega U+(\omega -1)D]\mathbf {x} ^{(k)}+\omega \mathbf {b} {\big )}=C_{\omega }\mathbf {x} ^{(k)}+\mathbf {c} ;$ dónde $\mathbf {x} ^{(k)}$ es la k-ésima aproximación de $\mathbf {x}$ , y $\mathbf {x} ^{(k+1)}$ es la nueva estimación que se quiere determinar. Como se puede ver, la matriz de iteración del método es: $C_{\omega }=-(D+\omega L)^{-1}{\big (}\omega U+(\omega -1)D{\big )}.$

Formulación por coordenadas

En la práctica se evita hallar la inversa de forma explícita al aplicar SOR. En su lugar se puede resolver el sistema de ecuaciones lineales que se obtiene al multiplicar cada lado de la iteración por ${\big (}D+\omega L{\big )}$ a la izquierda:

${\big (}D+\omega L{\big )}\mathbf {x} ^{(k+1)}=-[\omega U+(\omega -1)D]\mathbf {x} ^{(k)}+\omega \mathbf {b} ;$

Dado que la matriz ${\big (}D+\omega L{\big )}$ es triangular inferior, se puede hallar $x^{(k+1)}$ mediante sustitución hacia adelante. De esta forma se obtiene una expresión para cada coordenada de la estimación $x^{(k+1)}$ :

x_{i}^{(k+1)}=(1-\omega )x_{i}^{(k)}+{\frac {\omega }{a_{ii}}}\left(b_{i}-\sum _{j<i}a_{ij}x_{j}^{(k+1)}-\sum _{j>i}a_{ij}x_{j}^{(k)}\right),\quad i=1,2,\ldots ,n.

Vínculo con el método de Gauss-Seidel

Como se puede ver en la expresión anterior, si se toma $\omega =1$ , se obtiene el método de Gauss-Seidel como caso particular de SOR. Para los demás valores de $\omega$ , la estimación de SOR es una combinación convexa de la estimación SOR del paso anterior, y la estimación que se obtendría aplicando un paso de Gauss-Seidel a la estimación SOR del paso anterior:

$x^{(k+1)}=(1-\omega )x^{(k)}+\omega GS{\big (}x^{(k)}{\big )}.$

Convergencia

Una condición necesaria para que el algoritmo SOR sea convergente, es que se cumpla que $0<\omega <2$ .^[2] En general esta condición no es suficiente para asegurar la convergencia, aunque sí lo es para cierto tipo de matrices. En 1947, Ostrowski probó que si $A$ es simétrica y definida positiva, el método SOR converge si y solo si $0<\omega <2$ .^[2]

Tasa de convergencia

Generalmente es deseable seleccionar $\omega$ de forma que el método no solo sea convergente, sino que logre la convergencia lo más rápido posible. El valor de $\omega$ con el que se logra la mejor tasa de convergencia se denomina factor de relajación óptimo, y se denota $\omega _{\text{opt}}$ . Esta tasa o velocidad de convergencia viene dada por el recíproco del radio espectral $\rho {\big (}C_{\omega }{\big )}$ de la matriz de iteración, por lo que encontrar la mejor tasa de convergencia se reduce a minimizar $\rho {\big (}C_{\omega }{\big )}$ como una función de $\omega$ . Por lo tanto, la elección de $\omega$ no siempre es sencilla, y depende de las propiedades de la matriz $A$ del sistema.

Bajo ciertas hipótesis, es posible obtener una expresión analítica para el radio espectral $\rho {\big (}C_{\omega }{\big )}$ , y hallar en particular el $\omega _{\text{opt}}$ donde éste se minimiza (mayor tasa de convergencia). Supongamos en particular que se cumplen las siguientes hipótesis:^[3]^[4]

el parámetro de relajación cumple la condición necesaria de convergencia: $\omega \in (0,2)$ ,
los valores propios de la matriz de iteración de Jacobi, $C_{\text{Jac}}:=I-D^{-1}A$ son todos reales,
el método de Jacobi es convergente: $\mu :=\rho (C_{\text{Jac}})<1$ , y
la descomposición matricial $A=D+L+U$ satisface la propiedad que $\operatorname {det} (\lambda D+zL+{\tfrac {1}{z}}U)=\operatorname {det} (\lambda D+L+U)$ para todo $z\in \mathbb {C} \setminus \{0\}$ y $\lambda \in \mathbb {C}$ .

Entonces el radio espectral de la matriz de iteración de SOR puede ser expresado como:

\rho (C_{\omega })={\begin{cases}{\frac {1}{4}}\left(\omega \mu +{\sqrt {\omega ^{2}\mu ^{2}-4(\omega -1)}}\right)^{2}\,,&0<\omega \leq \omega _{\text{opt}}\\\omega -1\,,&\omega _{\text{opt}}<\omega <2\end{cases}};

donde el parámetro de relajación óptimo que minimiza el radio espectral es

\omega _{\text{opt}}:=1+\left({\frac {\mu }{1+{\sqrt {1-\mu ^{2}}}}}\right)^{2}\,.

En particular, para $\omega =1$ (Gauss-Seidel) se cumple que $\rho (C_{\omega })=\mu ^{2}=\rho (C_{\text{Jac}})^{2}$ .

La última hipótesis se satisface para matrices tridiagonales, pues $Z(\lambda D+L+U)Z^{-1}=\lambda D+zL+{\tfrac {1}{z}}U$ para la matriz diagonal $Z$ con componentes $Z_{ii}=z^{i-1}$ , y $\operatorname {det} (\lambda D+L+U)=\operatorname {det} (Z(\lambda D+L+U)Z^{-1})$ .

Pseudo-código del algoritmo

Dado que las estimaciones del algoritmo pueden ser sobre-escritas cuando están siendo computadas, la implementación del algoritmo requiere un único vector de almacenamiento en cada iteración. Por lo tanto, en el siguiente pseudo-código se omite la indexación de cada vector.

Entradas:  $A$ ,  $b$ ,  $ω$ 
Salida:  $\phi$

Escoger una estimación inicial  $\phi$  de la solución
while (no converge)
    for i de 1 hasta n hacer
         $\sigma \leftarrow 0$ 
        for j de 1 hasta n hacer
            if  $j$  ≠ i entonces
                $\sigma \leftarrow \sigma +a_{ij}\phi _{j}$ 
            end if
        end (j-bucle)
         $\phi _{i}\leftarrow (1-\omega )\phi _{i}+{\frac {\omega }{a_{ii}}}\left(b_{i}-\sigma \right)$ 
    end (i-bucle)
    verificar posible convergencia
end

Nota: El término: $(1-\omega )\phi _{i}+{\frac {\omega }{a_{ii}}}(b_{i}-\sigma )$ , puede ser escrito como: $\phi _{i}+\omega \left({\frac {b_{i}-\sigma }{a_{ii}}}-\phi _{i}\right)$ . De esta forma se ahorra una multiplicación en cada iteración del bucle for exterior.

Ejemplo

Se considera el siguiente sistema de ecuaciones, de tamaño $n=4$ :

{\begin{aligned}4x_{1}-x_{2}-6x_{3}+0x_{4}=2,\\-5x_{1}-4x_{2}+10x_{3}+8x_{4}=21,\\0x_{1}+9x_{2}+4x_{3}-2x_{4}=-12,\\1x_{1}+0x_{2}-7x_{3}+5x_{4}=-6.\end{aligned}}

Para estimar su solución se aplica SOR con un factor de relajación $\omega =0.5$ y estimación inicial $\phi =(0,0,0,0)$ . En cada iteración se obtienen las estimaciones que se muestran en la siguiente tabla. El método converge a la solución exacta $(3, -2, 2, 1)$ .

Iteración	$x_{1}$	$x_{2}$	$x_{3}$	$x_{4}$
0	0	0	0	0
1	0.25	-2.78125	1.6289062	0.5152344
2	1.2490234	-2.2448974	1.9687712	0.9108547
3	2.070478	-1.6696789	1.5904881	0.76172125
...	...	...	...	...
37	2.9999998	-2.0	2.0	1.0
38	3.0	-2.0	2.0	1.0

Implementación en Python

Una implementación en Python del pseudo-código proporcionado más arriba es:

import numpy as np

def sor_solver(A, b, omega, initial_guess, convergence_criteria):
    """
    This is an implementation of the pseudo-code provided in the Wikipedia article.
    Arguments:
        A: nxn numpy matrix.
        b: n dimensional numpy vector.
        omega: relaxation factor.
        initial_guess: An initial solution guess for the solver to start with.
        convergence_criteria: The maximum discrepancy acceptable to regard the current solution as fitting.
    Returns:
        phi: solution vector of dimension n.
    """
    phi = initial_guess[:]
    residual = np.linalg.norm(np.matmul(A, phi) - b) #Initial residual
    while residual > convergence_criteria:
        for i in range(A.shape[0]):
            sigma = 0
            for j in range(A.shape[1]):
                if j != i:
                    sigma += A[i][j] * phi[j]
            phi[i] = (1 - omega) * phi[i] + (omega / A[i][i]) * (b[i] - sigma)
        residual = np.linalg.norm(np.matmul(A, phi) - b)
        print('Residual: {0:10.6g}'.format(residual))
    return phi

# An example case that mirrors the one in the Wikipedia article
residual_convergence = 1e-8
omega = 0.5 #Relaxation factor

A = np.matrix([[4, -1, -6, 0],
              [-5, -4, 10, 8],
              [0, 9, 4, -2],
              [1, 0, -7, 5]])

b = np.matrix([2, 21, -12, -6])

initial_guess = np.zeros(4)

phi = sor_solver(A, b, omega, initial_guess, residual_convergence)
print(phi)

Véase también

Referencias

↑ Young, David M. (1 de mayo de 1950), Iterative methods for solving partial difference equations of elliptical type, PhD thesis, Harvard University, consultado el 15 de junio de 2009 .
↑ ^a ^b Demmel, J. W. (1997). Applied Numerical Linear Algebra. Society for Industrial and Applied Mathematics (SIAM).
↑ Hackbusch, Wolfgang (2016). «4.6.2». Iterative Solution of Large Sparse Systems of Equations | SpringerLink. Applied Mathematical Sciences (en inglés británico) 95. ISBN 978-3-319-28481-1. doi:10.1007/978-3-319-28483-5.
↑ Greenbaum, Anne (1997). «10.1». Iterative Methods for Solving Linear Systems. Frontiers in Applied Mathematics (en inglés británico) 17. ISBN 978-0-89871-396-1. doi:10.1137/1.9781611970937.

Datos: Q1561270

[1] Young, David M. (1 de mayo de 1950), Iterative methods for solving partial difference equations of elliptical type, PhD thesis, Harvard University, consultado el 15 de junio de 2009 .

[:0-2] Demmel, J. W. (1997). Applied Numerical Linear Algebra. Society for Industrial and Applied Mathematics (SIAM).

[3] Hackbusch, Wolfgang (2016). «4.6.2». Iterative Solution of Large Sparse Systems of Equations | SpringerLink. Applied Mathematical Sciences (en inglés británico) 95. ISBN 978-3-319-28481-1. doi:10.1007/978-3-319-28483-5.

[4] Greenbaum, Anne (1997). «10.1». Iterative Methods for Solving Linear Systems. Frontiers in Applied Mathematics (en inglés británico) 17. ISBN 978-0-89871-396-1. doi:10.1137/1.9781611970937.

[1]

[2]

[3]

[4]