Tipo de problemas em que o SOR é mais rápido que o Gauss-Seidel?

9

Existe alguma regra simples para dizer se vale a pena fazer SOR em vez de Gauss-Seidel? (e possível maneira de estimar o parâmetro de realxação )ω

Quero dizer, apenas olhando para a matriz , ou conhecimento de um problema específico que a matriz representa?

Eu estava lendo a resposta para estas perguntas: Existem heurísticas para otimizar o método de relaxamento excessivo (SOR) sucessivo? mas é um pouco sofisticado demais. Não vejo heurísticas simples como estimar o raio espectral apenas olhando para a matriz (ou problema que ela representa).

Gostaria de algo muito mais simples - apenas alguns exemplos de matrizes (problemas) para os quais o SOR converge mais rapidamente.


Eu estava experimentando o SOR para a matriz desse rei: onde eu sou matriz de identidade, C i j = c i , j e R i j s são números aleatórios da distribuição unifrom, tal que | R i j | < r . Eu estava pensando que haveria alguma dependência do ω ótimo nos parâmetros c , r .A=I+C+RICij=c i,jRij|Rij|<rωc,r

EDIT: Usei muito pequeno para garantir que A seja fortemente dominante na diagonal. ( | c | < 0,1 , r < 2 | c | para matriz da dimensão 5-10). Devo dizer também que estes A eram reais e simétricos.c,rA|c|<0.1r<2|c|A

No entanto, descobri que Gauss-Seidel ( ) é quase sempre o melhor (?)ω=1 . Isso significa que deve haver mais correlação entre s para obter vantagem do SOR? Ou eu fiz algo errado? Aij


Eu sei que o SOR não é o solucionador mais eficiente (em comparação com CG, GMRES ...), mas é simples de implementar, paraelizar e modificar para um problema específico. Certamente bom para prototipagem.

Prokop Hapala
fonte

Respostas:

5

A convergência de solucionadores iterativos clássicos para sistemas lineares é determinada pelo raio espectral da matriz de iteração, . Para um sistema linear geral, é difícil determinar um parâmetro SOR ideal (ou até bom) devido à dificuldade em determinar o raio espectral da matriz de iteração. Abaixo, incluí muitos detalhes adicionais, incluindo um exemplo de um problema real em que o peso ideal da SOR é conhecido.ρ(G)

Raio espectral e convergência

O raio espectral é definido como o valor absoluto do maior valor próprio de magnitude. Um método convergirá se e um raio espectral menor significar uma convergência mais rápida. O SOR funciona alterando a divisão da matriz usada para derivar a matriz de iteração com base na escolha de um parâmetro de ponderação ω , diminuindo, esperançosamente, o raio espectral da matriz de iteração resultante.ρ<1ω

Divisão de matriz

Para a discussão abaixo, assumirei que o sistema a ser resolvido é dado por

Ax=b,

com uma iteração do formulário

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

onde é um vetor e o número da iteração k é denotado x ( k ) .vkx(k)

O SOR leva uma média ponderada da iteração antiga e uma iteração de Gauss-Seidel. O método Gauss-Seidel se baseia em uma divisão de matriz da forma

A=D+L+U

onde é a diagonal de A , L é uma matriz triangular inferior contendo todos os elementos de A estritamente abaixo da diagonal e R é uma matriz triangular superior contendo todos os elementos de A estritamente acima da diagonal. A iteração de Gauss-Seidel é então dada porDALARA

x(k+1)=(D+L)1b+GGSx(k)

e a matriz de iteração é

GGS=(D+L)1U.

O SOR pode então ser escrito como

x(k+1)=ω(D+ωL)1b+GSORx(k)

Onde

GSOR=(D+ωL)1((1ω)DωU).

ω

SOR ideal

Um exemplo realista em que o coeficiente de ponderação ideal é conhecido surge no contexto da resolução de uma equação de Poisson:

2u=f in Ωu=g on Ω

A discretização desse sistema em um domínio quadrado em 2D usando diferenças finitas de segunda ordem com espaçamento uniforme de grade resulta em uma matriz simétrica com 4 na diagonal, -1 imediatamente acima e abaixo da diagonal e mais duas bandas de -1 a alguma distância do diagonal. Existem algumas diferenças devido às condições de contorno, mas essa é a estrutura básica. Dada essa matriz, a escolha comprovadamente ótima para o coeficiente SOR é dada por

ω=21+sin(πΔx/L)

ΔxL

Erro de Gauss-Seidel e SOR

Como você pode ver, o SOR atinge a precisão da máquina em cerca de 100 iterações, momento em que Gauss-Seidel é cerca de 25 ordens de magnitude pior. Se você quiser brincar com este exemplo, incluí o código MATLAB que usei abaixo.

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
fonte
Você conhece alguma técnica boa / conhecida que é usada para calcular o parâmetro SOR em tempo real? Ouvi dizer que essas técnicas usam estimativas do raio espectral - você poderia explicar como elas usam o raio espectral ou fornecer uma boa referência?
Nkeguy
Ah, vejo que isso é abordado na pergunta vinculada scicomp.stackexchange.com/questions/851/… . Não importa as minhas perguntas, mas se você tiver mais a acrescentar, sinta-se à vontade para fazê-lo.
Nukeguy 21/07
@Doug Lipinski Eu pensei que f deveria ser multiplicado por dx * dy. Esse fator vem da segunda derivada discreta (veja aqui, por exemplo). Btw, quando eu faço isso, o algoritmo não funciona corretamente. Você sabe por quê?
shamalaia
0

Esse lado das coisas não é realmente minha especialidade, mas não acho que seja um teste super justo para muitas aplicações realistas.

Não tenho certeza de quais valores você estava usando para c e r , mas suspeito que você estava trabalhando com matrizes extremamente mal condicionadas. (Abaixo está um código Python mostrando que essas podem não ser as matrizes mais invertíveis.)

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

Se você realmente precisasse inverter matrizes dessa condição, a) usaria um método especializado eb) provavelmente deveria apenas encontrar um novo campo 😉

Para matrizes bem condicionadas de qualquer tamanho, o SOR provavelmente será mais rápido. Para problemas reais em que a velocidade é importante, seria raro usar o SOR - do lado sofisticado, há muito melhor nos dias de hoje; No lado lento, mas confiável, o SOR não é o melhor que você pode fazer.

Mike
fonte
0.01<|c|<0.1r<2|c|
Eu ia dizer fortemente diagonalmente dominante.
precisa saber é o seguinte
0

OK, então para matrizes simétricas deste rei:

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

ti=c+random(r,r)

tc=0,r=0.1t

(Isso é apenas observação empírica, nada rigoroso)

Prokop Hapala
fonte