¿Hay alguna regla general simple para decir si vale la pena hacer SOR en lugar de Gauss-Seidel? (y posible forma de estimar el parámetro de realxation )$\omega$

Me refiero solo a mirar la matriz , o al conocimiento de un problema particular que representa la matriz.

Estaba leyendo la respuesta a estas preguntas: ¿Hay alguna heurística para optimizar el método de relajación excesiva sucesiva (SOR)? pero es un poco demasiado sofisticado No veo heurísticas simples sobre cómo estimar el radio espectral simplemente mirando la matriz (o el problema que representa).

Me gustaría algo mucho más simple: solo algunos ejemplos de matrices (problemas) para los cuales SOR convergen más rápido.

---

Estaba experimentando con SOR para la matriz de este rey: donde es la matriz de identidad, y s son números aleatorios de distribución unifrom tal que . Estaba pensando que habría cierta dependencia de óptima en los parámetros .I C i j = c ∀ i , j R i j | R i j | < r ω c , $A = I + C + R$$I$$C_{ij}=c$ $\forall i,j$$R_{ij}$$|R_{ij}|<r$$\omega$$c,r$

EDITAR: utilicé muy pequeño para asegurarme de que sea ​​fuertemente diagonalmente dominante. ( , para la matriz de dimensión 5-10). También debería decir que estas eran reales y simétricas.A | c | < 0.1 r < 2 | c | $c,r$$A$$|c|<0.1$$r<2|c|$$A$

Sin embargo, descubrí que Gauss-Seidel ( ) es casi siempre el mejor (?)$\omega=1$ . ¿Esto significa que tiene que haber alguna correlación más entre s para aprovechar SOR? ¿O hice algo mal? $A_{ij}$

---

Sé que SOR no es el solucionador más eficiente (en comparación con CG, GMRES ...) pero es fácil de implementar y paraelizar, y modificar para un problema particular. Seguro bueno para la creación de prototipos.

La convergencia de los solucionadores iterativos clásicos para sistemas lineales está determinada por el radio espectral de la matriz de iteración, . Para un sistema lineal general, es difícil determinar un parámetro SOR óptimo (o incluso bueno) debido a la dificultad para determinar el radio espectral de la matriz de iteración. A continuación, he incluido muchos detalles adicionales, incluido un ejemplo de un problema real en el que se conoce el peso óptimo de SOR.$\rho(\mathbf{G})$

Radio espectral y convergencia

El radio espectral se define como el valor absoluto del valor propio de mayor magnitud. Un método convergerá si y un radio espectral más pequeño significa una convergencia más rápida. SOR funciona alterando la división de la matriz utilizada para derivar la matriz de iteración en función de la elección de un parámetro de ponderación ω , con suerte disminuyendo el radio espectral de la matriz de iteración resultante.$\rho<1$$\omega$

División de la matriz

Para la discusión a continuación, supondré que el sistema a resolver está dado por

$$\mathbf{A}x=b,$$

con una iteración de la forma

$$x^{(k+1)}=v+\mathbf{G}x^{(k)},$$

donde es un vector, y el número de iteración k se denota x ( k ) .$v$$k$$x^{(k)}$

SOR toma un promedio ponderado de la iteración anterior y una iteración de Gauss-Seidel. El método Gauss-Seidel se basa en una división matricial de la forma

$$\mathbf{A} = \mathbf{D} + \mathbf{L} + \mathbf{U}$$

donde es la diagonal de A , L es una matriz triangular inferior que contiene todos los elementos de A estrictamente por debajo de la diagonal y R es una matriz triangular superior que contiene todos los elementos de A estrictamente por encima de la diagonal. La iteración de Gauss-Seidel viene dada por$\mathbf{D}$$\mathbf{A}$$\mathbf{L}$$\mathbf{A}$$\mathbf{R}$$\mathbf{A}$

$$x^{(k+1)}=(\mathbf{D}+\mathbf{L})^{-1}b+\mathbf{G}_{\rm G-S}x^{(k)}$$

y la matriz de iteración es

$$\mathbf{G}_{\rm G-S}=-(\mathbf{D}+\mathbf{L})^{-1}\mathbf{U}.$$

SOR puede entonces escribirse como

$$x^{(k+1)}=\omega(\mathbf{D}+\omega\mathbf{L})^{-1}b+\mathbf{G}_{\rm SOR}x^{(k)}$$

dónde

$$\mathbf{G}_{\rm SOR} = (\mathbf{D}+\omega\mathbf{L})^{-1}((1-\omega)\mathbf{D}-\omega\mathbf{U}).$$

$\omega$

SOR óptimo

Un ejemplo realista donde se conoce el coeficiente de ponderación óptimo surge en el contexto de la resolución de una ecuación de Poisson:

$$\nabla^2 u = f~{\rm in}~\Omega\\ u = g~{\rm on}~\partial\Omega$$

Discretizar este sistema en un dominio cuadrado en 2D usando diferencias finitas de segundo orden con un espaciado de cuadrícula uniforme da como resultado una matriz simétrica con bandas con 4 en la diagonal, -1 inmediatamente arriba y debajo de la diagonal, y dos bandas más de -1 a cierta distancia del diagonal. Hay algunas diferencias debido a las condiciones de contorno, pero esa es la estructura básica. Dada esta matriz, la opción demostrablemente óptima para el coeficiente SOR viene dada por

$$\omega = \dfrac{2}{1+\sin(\pi \Delta x/L)}$$

$\Delta x$$L$

Como puede ver, SOR alcanza la precisión de la máquina en aproximadamente 100 iteraciones, momento en el que Gauss-Seidel es aproximadamente 25 órdenes de magnitud peor. Si quieres jugar con este ejemplo, he incluido el código MATLAB que utilicé a continuación.

clear all
close all

%number of iterations:
niter = 150;

%number of grid points in each direction
N = 16;
% [x y] = ndgrid(linspace(0,1,N),linspace(0,1,N));
[x y] = ndgrid(linspace(-pi,pi,N),linspace(-pi,pi,N));
dx = x(2,1)-x(1,1);
L = x(N,1)-x(1,1);

%desired solution:
U = sin(x/2).*cos(y);

% Right hand side for the Poisson equation (computed from U to produce the
% desired known solution)
Ix = 2:N-1;
Iy = 2:N-1;
f = zeros(size(U));
f(Ix,Iy) = (-4*U(Ix,Iy)+U(Ix-1,Iy)+U(Ix+1,Iy)+U(Ix,Iy-1)+U(Ix,Iy+1));

figure(1)
clf
contourf(x,y,U,50,'linestyle','none')
title('True solution')

%initial guess (must match boundary conditions)
U0 = U;
U0(Ix,Iy) = rand(N-2);

%Gauss-Seidel iteration:
UGS = U0; EGS = zeros(1,niter);
for iter=1:niter
    for iy=2:N-1
        for ix=2:N-1
            UGS(ix,iy) = -1/4*(f(ix,iy)-UGS(ix-1,iy)-UGS(ix+1,iy)-UGS(ix,iy-1)-UGS(ix,iy+1));
        end
    end

    %error:
    EGS(iter) = sum(sum((U-UGS).^2))/sum(sum(U.^2));
end

figure(2)
clf
contourf(x,y,UGS,50,'linestyle','none')
title(sprintf('Gauss-Seidel approximate solution, iteration %d', iter))
drawnow

%SOR iteration:
USOR = U0; ESOR = zeros(1,niter);
w = 2/(1+sin(pi*dx/L));
for iter=1:niter
    for iy=2:N-1
        for ix=2:N-1
            USOR(ix,iy) = (1-w)*USOR(ix,iy)-w/4*(f(ix,iy)-USOR(ix-1,iy)-USOR(ix+1,iy)-USOR(ix,iy-1)-USOR(ix,iy+1));
        end
    end

    %error:
    ESOR(iter) = sum(sum((U-USOR).^2))/sum(sum(U.^2));
end

figure(4)
clf
contourf(x,y,USOR,50,'linestyle','none')
title(sprintf('Gauss-Seidel approximate solution, iteration %d', iter))
drawnow


figure(5)
clf
semilogy(EGS,'b')
hold on
semilogy(ESOR,'r')
title('L2 relative error')
xlabel('Iteration number')
legend('Gauss-Seidel','SOR','location','southwest')

Este lado de las cosas no es realmente mi especialidad, pero no creo que sea una prueba súper justa para muchas aplicaciones realistas.

No estoy seguro de qué valores estaba usando para c y r , pero sospecho que estaba trabajando con matrices extremadamente mal condicionadas. (A continuación se muestra un código de Python que muestra que estas podrían no ser las matrices más invertibles).

>>> import numpy
>>> for n in [100, 1000]:
...     for c in [.1, 1, 10]:
...         for r in [.1, 1, 10]:
...             print numpy.linalg.cond(
...                 numpy.eye(n) + 
...                 c * numpy.ones((n, n)) + 
...                 r * numpy.random.random((n, n))
...             )
... 
25.491634739
2034.47889101
2016.33059429
168.220149133
27340.0090644
5532.81258852
1617.33518781
42490.4410689
5326.3865534
6212.01580004
91910.8386417
50543.9269739
24737.6648458
271579.469212
208913.592289
275153.967337
17021788.5576
117365.924601

Si realmente necesitara invertir matrices tan mal acondicionadas, a) usaría un método especializado yb) probablemente debería buscar un nuevo campo 😉

Para matrices bien acondicionadas de cualquier tamaño, es probable que SOR sea más rápido. Para problemas reales en los que la velocidad es importante, sería raro usar SOR: en el lado sofisticado, hay mucho mejor en estos días; en el lado lento pero confiable, SOR no es lo mejor que puede hacer.

OK, entonces para matrices simétricas de este rey:

1 t 0 0 0 0 t t 0 0 
t 1 0 0 0 0 0 0 0 0 
0 0 1 0 0 0 0 t 0 t 
0 0 0 1 0 0 0 0 0 t 
0 0 0 0 1 t 0 0 0 0 
0 0 0 0 t 1 0 t 0 0 
t 0 0 0 0 0 1 0 0 0 
t 0 t 0 0 t 0 1 0 0 
0 0 0 0 0 0 0 0 1 t 
0 0 t t 0 0 0 0 t 1 

$t$$t$$t$

$t_i = c + random(-r,r)$

$t$$c=0,r=0.1$$t$

(Esto es solo observación emperical, nada riguroso)