Y a-t-il une règle empirique simple pour dire s'il vaut la peine de faire le SOR au lieu de Gauss-Seidel? (et manière possible d'estimer le paramètre de réallocation )$\omega$

Je veux dire juste en regardant la matrice , ou la connaissance d' un problème particulier que la matrice représente?

Je lisais la réponse à ces questions: Existe-t-il des heuristiques pour optimiser la méthode de sur-relaxation successive (SOR)? mais c'est un peu trop sophistiqué. Je ne vois pas d'heuristique simple comment estimer le rayon spectral en regardant simplement la matrice (ou le problème qu'elle représente).

Je voudrais quelque chose de beaucoup plus simple - juste quelques exemples de matrices (problèmes) pour lesquelles le SOR converge plus rapidement.

---

J'expérimentais avec SOR pour la matrice de ce roi: où I est la matrice d'identité, C i j = c ∀ i , j et R i j s sont des nombres aléatoires de distribution unifrom tels que | R i j | < r . Je pensais qu'il y aurait une certaine dépendance du ω optimal aux paramètres c , r .$A = I + C + R$$I$$C_{ij}=c$ $\forall i,j$$R_{ij}$$|R_{ij}|<r$$\omega$$c,r$

EDIT: J'ai utilisé de très petits pour m'assurer que A est fortement dominant en diagonale. ( | c | < 0,1 , r < 2 | c | pour la matrice de dimension 5-10). Je dois également dire que ces A étaient réels et symétriques.$c,r$$A$$|c|<0.1$$r<2|c|$$A$

Cependant, j'ai trouvé que Gauss-Seidel ( ) est presque toujours le meilleur (?)$\omega=1$ . Cela signifie-t-il qu'il doit y avoir une corrélation plus importante entre s pour bénéficier du SOR? Ou j'ai fait quelque chose de mal? $A_{ij}$

---

Je sais que SOR n'est pas le solveur le plus efficace (par rapport à CG, GMRES ...) mais il est simple à implémenter et à paramétrer et à modifier pour un problème particulier. Bien sûr pour le prototypage.

La convergence des solveurs itératifs classiques pour les systèmes linéaires est déterminée par le rayon spectral de la matrice d'itération, .$\rho(\mathbf{G})$Pour un système linéaire général, il est difficile de déterminer un paramètre SOR optimal (ou même bon) en raison de la difficulté à déterminer le rayon spectral de la matrice d'itération. Ci-dessous, j'ai inclus de nombreux détails supplémentaires, y compris un exemple d'un problème réel où le poids SOR optimal est connu.

Rayon spectral et convergence

Le rayon spectral est défini comme la valeur absolue de la valeur propre de plus grande amplitude. Une méthode convergera si et un rayon spectral plus petit signifie une convergence plus rapide. SOR fonctionne en modifiant le fractionnement de la matrice utilisé pour dériver la matrice d'itération en fonction du choix d'un paramètre de pondération $\rho<1$$\omega$ , ce qui, espérons-le, diminue le rayon spectral de la matrice d'itération résultante.

Fractionnement de matrice

Pour la discussion ci-dessous, je suppose que le système à résoudre est donné par

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

avec une itération du formulaire

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

où est un vecteur et le numéro d'itération k est noté x ( k ) .$v$$k$$x^{(k)}$

SOR prend une moyenne pondérée de l'ancienne itération et une itération de Gauss-Seidel. La méthode de Gauss-Seidel repose sur un découpage matriciel de la forme

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

où est la diagonale de A , L est une matrice triangulaire inférieure contenant tous les éléments de A strictement en dessous de la diagonale et R est une matrice triangulaire supérieure contenant tous les éléments de A strictement au-dessus de la diagonale. L'itération de Gauss-Seidel est alors donnée par$\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)}$$

et la matrice d'itération est

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

Le SOR peut alors s'écrire

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

où

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

$\omega$

SOR optimal

Un exemple réaliste où le coefficient de pondération optimal est connu se présente dans le contexte de la résolution d'une équation de Poisson:

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

La discrétisation de ce système sur un domaine carré en 2D en utilisant des différences finies de second ordre avec un espacement de grille uniforme donne une matrice à bandes symétriques avec 4 sur la diagonale, -1 immédiatement au-dessus et en dessous de la diagonale, et deux autres bandes de -1 à une certaine distance du diagonale. Il y a quelques différences en raison des conditions aux limites, mais c'est la structure de base. Étant donné cette matrice, le choix prouvablement optimal pour le coefficient SOR est donné par

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

où $\Delta x$$L$

Comme vous pouvez le voir, le SOR atteint la précision de la machine en environ 100 itérations, moment auquel Gauss-Seidel est environ 25 ordres de grandeur pire. Si vous voulez jouer avec cet exemple, j'ai inclus le code MATLAB que j'ai utilisé ci-dessous.

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

Ce côté des choses n'est pas vraiment ma spécialité, mais je ne pense pas que ce soit un test super juste pour de nombreuses applications réalistes.

Je ne sais pas pour quelles valeurs vous utilisiez c et r , mais je soupçonne que vous travailliez avec des matrices extrêmement mal conditionnées. (Ci-dessous, un code Python montrant que ce ne sont peut-être pas les matrices les plus inversibles.)

>>> 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 vous aviez réellement besoin d'inverser des matrices aussi mal conditionnées, vous auriez a) utiliser une méthode spécialisée, et b) devrait probablement simplement trouver un nouveau champ 😉

Pour les matrices bien conditionnées de toute taille, le SOR est susceptible d'être plus rapide. Pour de vrais problèmes où la vitesse est importante, il serait rare d'utiliser SOR - du côté sophistiqué, il y a beaucoup mieux de nos jours; Du côté lent mais fiable, le SOR n'est pas le meilleur que vous puissiez faire.

OK, donc pour les matrices symétriques de ce roi:

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$s sont similaires. J'utilisais$t$s généré comme ceci:

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

Si $t$s varient beaucoup et sont centrés autour de 0 ( $c=0,r=0.1$) que Gauss-Seidel est plus rapide. Gauss-Seidel est également plus rapide si chaque ligne est remplie à plus de la moitié par$t$s. Cela signifie également que le SOR est meilleur pour les matrices très grandes et très clairsemées.

(Ce n'est qu'une observation emperique, rien de rigoureux)