Probleme, bei denen SOR schneller ist als Gauß-Seidel?

9

Gibt es eine einfache Faustregel, um zu sagen, ob es sich lohnt, SOR anstelle von Gauß-Seidel durchzuführen? (und mögliche Methode zur Schätzung des Realxationsparameters )ω

Ich meine, nur indem ich auf die Matrix schaue oder das Wissen über ein bestimmtes Problem, das die Matrix darstellt?

Ich habe die Antwort auf diese Fragen gelesen: Gibt es Heuristiken zur Optimierung der Methode der sukzessiven Überentspannung (SOR)? aber es ist ein bisschen zu raffiniert. Ich sehe keine einfachen Heuristiken, wie man den Spektralradius schätzt, wenn man nur auf die Matrix (oder das Problem, das sie darstellt) schaut.

Ich möchte etwas viel Einfacheres - nur einige Beispiele für Matrizen (Probleme), für die SOR schneller konvergiert.


Ich habe mit SOR für die Matrix dieses Königs experimentiert: wobei I die Identitätsmatrix ist, C i j = c i , j und R i j s Zufallszahlen aus einer einheitlichen Verteilung sind, so dass | R i j | < r . Ich dachte, dass es eine gewisse Abhängigkeit des optimalen ω von den Parametern c , r geben würde .A=I+C+RichC.ichj=c ich,jR.ichj|R.ichj|<rωc,r

EDIT: Ich habe sehr kleine , um sicherzustellen, dass A stark diagonal dominant ist. ( | c | < 0,1 , r < 2 | c | für Matrix der Dimension 5-10). Ich sollte auch sagen, dass diese A real und symmetrisch waren.c,rEIN|c|<0,1r<2|c|EIN

Ich fandω=1 jedoch, dass Gauß-Seidel ( ω = 1 ) fast immer das Beste ist (?) . Bedeutet dies, dass es eine größere Korrelation zwischen s geben muss, um SOR nutzen zu können? Oder habe ich etwas falsch gemacht? EINichj


Ich weiß, dass SOR nicht der effizienteste Löser ist (im Vergleich zu CG, GMRES ...), aber es ist einfach zu implementieren, zu paraelisieren und für bestimmte Probleme zu modifizieren. Sicher gut für das Prototyping.

Prokop Hapala
quelle

Antworten:

5

Die Konvergenz klassischer iterativer Löser für lineare Systeme wird durch den Spektralradius der Iterationsmatrix . Für ein allgemeines lineares System ist es schwierig, einen optimalen (oder sogar guten) SOR-Parameter zu bestimmen, da es schwierig ist, den Spektralradius der Iterationsmatrix zu bestimmen. Im Folgenden habe ich viele zusätzliche Details aufgeführt, einschließlich eines Beispiels für ein echtes Problem, bei dem das optimale SOR-Gewicht bekannt ist.ρ(G)

Spektralradius und Konvergenz

Der Spektralradius ist definiert als der Absolutwert des Eigenwerts der größten Größe. Eine Methode konvergiert, wenn und ein kleinerer Spektralradius eine schnellere Konvergenz bedeutet. SOR ändert die Matrixaufteilung, die zum Ableiten der Iterationsmatrix verwendet wird, basierend auf der Wahl eines Gewichtungsparameters ω , wodurch hoffentlich der Spektralradius der resultierenden Iterationsmatrix verringert wird.ρ<1ω

Matrixaufteilung

Für die folgende Diskussion gehe ich davon aus, dass das zu lösende System gegeben ist durch

EINx=b,

mit einer Iteration des Formulars

x(k+1)=v+Gx(k),

wobei ein Vektor ist und die Iterationszahl k mit x ( k ) bezeichnet wird .vkx(k)

SOR nimmt einen gewichteten Durchschnitt der alten Iteration und eine Gauß-Seidel-Iteration. Die Gauß-Seidel-Methode beruht auf einer Matrixaufteilung der Form

EIN=D.+L.+U.

wobei die Diagonale von A ist , L eine untere Dreiecksmatrix ist, die alle Elemente von A streng unterhalb der Diagonale enthält, und R eine obere Dreiecksmatrix ist, die alle Elemente von A genau oberhalb der Diagonale enthält. Die Gauß-Seidel-Iteration ist dann gegeben durchD.EINL.EINR.EIN

x(k+1)=(D.+L.)- -1b+GG- -S.x(k)

und die Iterationsmatrix ist

GG- -S.=- -(D.+L.)- -1U..

SOR kann dann geschrieben werden als

x(k+1)=ω(D.+ωL.)- -1b+GS.ÖR.x(k)

wo

GS.ÖR.=(D.+ωL.)- -1((1- -ω)D.- -ωU.).

ω

Optimale SOR

Ein realistisches Beispiel, bei dem der optimale Gewichtungskoeffizient bekannt ist, ergibt sich im Zusammenhang mit der Lösung einer Poisson-Gleichung:

2u=f ichn Ωu=G Ön Ω

Die Diskretisierung dieses Systems auf einer quadratischen Domäne in 2D unter Verwendung endlicher Differenzen zweiter Ordnung mit gleichmäßigem Gitterabstand führt zu einer symmetrischen Bandmatrix mit 4 auf der Diagonale, -1 unmittelbar über und unter der Diagonale und zwei weiteren Bändern von -1 in einiger Entfernung von der diagonal. Aufgrund der Randbedingungen gibt es einige Unterschiede, aber das ist die Grundstruktur. In Anbetracht dieser Matrix ist die nachweislich optimale Wahl für den SOR-Koeffizienten gegeben durch

ω=21+Sünde(πΔx/.L.)

ΔxL.

Gauß-Seidel- und SOR-Fehler

Wie Sie sehen können, erreicht SOR die Maschinengenauigkeit in etwa 100 Iterationen, wobei Gauß-Seidel um etwa 25 Größenordnungen schlechter ist. Wenn Sie mit diesem Beispiel herumspielen möchten, habe ich den unten verwendeten MATLAB-Code eingefügt.

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')
Doug Lipinski
quelle
Kennen Sie gute / bekannte Techniken, mit denen der SOR-Parameter im laufenden Betrieb berechnet wird? Ich habe zuvor gehört, dass diese Techniken Schätzungen des Spektralradius verwenden - können Sie erklären, wie sie den Spektralradius verwenden, oder eine gute Referenz liefern?
Nukeguy
Oh, ich sehe, dass dies in der verknüpften Frage scicomp.stackexchange.com/questions/851/… angesprochen wird . Egal, meine Fragen, aber wenn Sie mehr hinzufügen möchten, können Sie dies gerne tun.
Nukeguy
@ Doug Lipinski Ich dachte, dass f mit dx * dy multipliziert werden sollte. Dieser Faktor stammt aus der diskreten zweiten Ableitung (siehe hier zum Beispiel). Übrigens, wenn ich es mache, funktioniert der Algorithmus nicht richtig. Weißt du, warum?
Shamalaia
0

Diese Seite der Dinge ist nicht wirklich meine Spezialität, aber ich denke nicht, dass dies ein super fairer Test für viele realistische Anwendungen ist.

Ich bin nicht sicher, welche Werte Sie für c und r verwendet haben , aber ich vermute, Sie haben mit extrem schlecht konditionierten Matrizen gearbeitet. (Nachfolgend finden Sie Python-Code, der zeigt, dass dies möglicherweise nicht die invertierbarsten Matrizen sind.)

>>> 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

Wenn Sie tatsächlich Matrizen invertieren müssten, die so schlecht konditioniert sind, würden Sie a) eine spezielle Methode verwenden und b) wahrscheinlich nur ein neues Feld suchen 😉

Bei gut konditionierten Matrizen jeder Größe ist SOR wahrscheinlich schneller. Für echte Probleme, bei denen es auf Geschwindigkeit ankommt, ist es selten, SOR zu verwenden - auf der anspruchsvollen Seite gibt es heutzutage viel bessere; Auf der langsamen, aber zuverlässigen Seite ist SOR nicht das Beste, was Sie tun können.

Mike
quelle
Hallo, ich sage nicht, dass mein "Test" fair ist. Ich würde nicht einmal sagen, dass es ein Test ist, es ist nur mein naiver Versuch, eine Vorstellung davon zu bekommen, wie sich SOR und Gauss-Seidel experimentell verhalten. Angenommen, ich bin ein absoluter Neuling auf diesem Gebiet. Meine Parameter lagen im Bereich von 0,01<|c|<0,1r<2|c|
Ich wollte stark diagonal dominant sagen.
Meawoppl
0

OK, also für symetrische Matrizen dieses Königs:

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 

ttt

tich=c+reinndÖm(- -r,r)

tc=0,r=0,1t

(Dies ist nur eine kaiserliche Beobachtung, nichts Strenges)

Prokop Hapala
quelle