Aplicação segura de métodos iterativos em matrizes diagonalmente dominantes

9

Suponha que o sistema linear a seguir seja dado que é o Laplaciano ponderado conhecido como positivo definido, com um espaço nulo unidimensional abrangido por , e a variação de conversão de , ou seja, não altera o valor da função (cuja derivada é ). As únicas entradas positivas de estão na diagonal, que é um somatório dos valores absolutos das entradas negativas fora da diagonal.

(1)Lx=c,
Lsemi1n=(1,,1)RnxRnx+a1n(1)L

Descobri em um trabalho acadêmico altamente citado em seu campo que, embora seja diagonalmente dominante, métodos como o Gradiente Conjugado, Gauss-Seidl, Jacobi ainda poderiam ser usados ​​com segurança para resolver . A lógica é que, devido à invariância da tradução, é seguro fixar um ponto (por exemplo, remova a primeira linha e coluna de e a primeira entrada de ), convertendo em uma matriz diagonalmente dominante. De qualquer forma, o sistema original é resolvido na forma completa de , com .Lnot strictly(1)LcLstrictly(1)LRn×n

Esta suposição está correta e, em caso afirmativo, qual é a lógica alternativa? Estou tentando entender como a convergência dos métodos ainda se mantém.

Se o método de Jacobi é convergente com , o que se poderia afirmar sobre o raio espectral da matriz de iteração , onde é a matriz diagonal com entradas de na diagonal? É , portanto diferente das garantias gerais de convergência para ? Estou perguntando isso, já que os valores próprios dos A matriz laplaciana com as na diagonal deve estar no intervalo .ρ D - 1 ( D - L ) D L ρ ( D - 1 ( D - L ) 1 ρ ( D - 1 ( D - L ) ) < 1 D - 1 L [ 0 , 2 ](1)ρD1(DL)DLρ(D1(DL)1ρ(D1(DL))<1D1L[0,2]

Do trabalho original:

......................................

Em cada iteração, calculamos um novo layout (x (t +1), y (t + 1)) resolvendo o seguinte sistema linear: Sem perda de generalidade, podemos fixar a localização de um dos os sensores (utilizando o grau de liberdade de translação da tensão localizada) e obtêm uma matriz estritamente diagonal dominante. Portanto, podemos usar com segurança a iteração de Jacobi para resolver (8)

(8)L·x(t+1)=L(x(t),y(t))·x(t)L·y(t+1)=L(x(t),y(t))·y(t)

.......................................

Acima, a noção de "iteração" está relacionada ao procedimento de minimização subjacente e não deve ser confundida com a iteração de Jacobi. Portanto, o sistema é resolvido por Jacobi (iterativamente) e, em seguida, a solução é comprada no lado direito de (8), mas agora para outra iteração da minimização subjacente. Espero que isso esclareça o assunto.

Observe que eu achei Quais solucionadores lineares iterativos convergem para matrizes semidefinidas positivas? , mas estou procurando uma resposta mais elaborada.

usero
fonte
Você poderia postar um link ou uma citação no trabalho mais citado?
precisa
Ele pode ser recuperado de: citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.164.1421 Como não se espera que você leia todo o trabalho, dê uma olhada na p.7 (abaixo). Suponho que a escolha de solucionadores iterativos seja justificada, mas sinto que é necessária uma lógica melhor (ou, pelo menos, diferente).
usero
Eu me pergunto se esses caras são da mesma comunidade que os pré-condicionadores combinatórios.
shuhalo

Respostas:

5

A iteração Jacobi pode ser comprovada convergente.

A primeira coisa que você deve ter certeza é que , que é a condição para a existência de solução (presumo L = L T , caso contrário, você precisa c ( K e r L T ) ), porque você disse V 0 : = K e r L = s p a n { 1 n } . Usaremos a convenção de que V 0cT1n=0L=LTc(KerLT)V0:=KerL=span{1n}V0é também a matriz com colunas sendo a base ortonormal dela. No seu caso, .V0:=1n/n

Então, para os erros da iteração Jacobi no sistema original, você tem onde P : = I -

e1=(ID1L)e0=(ID1L)(Pe0+V0a)=(ID1L)Pe0+V0a,
é a projeção ortogonal em V 1 : = V 0 . A partir da iteração acima, sabemos que P e 1 = P ( I - D - 1 L ) P e 0 , da qual temos a matriz de iteração S em V 1 , S : = P ( I - D - 1 L ) P . Isso nãoP:=IV0V0V1:=V0
Pe1=P(ID1L)Pe0,

SV1
S:=P(ID1L)P.
tem o mesmo espectro (exceto zeros) com a seguinte matriz ˜ S : = ( I - D - 1 L ) P P = ( I - D - 1 L ) P = ( I - D - 1 L ) ( I - V 0 V 0 )S Queremos que o raio espectral de S seja menor que um para provar a convergência.
S~:=(ID1L)PP=(ID1L)P=(ID1L)(IV0V0)=ID1LV0V0.
S

A citação a seguir é antiga e mantida apenas para referência. Veja depois para a nova prova.

No seu caso, E você pode verificar seD-1L+V0V0 é estritamente diagonal-dominante usando a suposição de que as entradas deLsão positivas na diagonal e negativas em caso contrário. Para mostrar os autovalores de D-1L+V0V0 são reais, notamos que a matriz é auto-adjunta sob o produto interno<x,y>:=yTDxV0V0=1n1n×n.D1L+V0V0LD1L+V0V0<x,y>:=yTDx.

Se não estiver no seu formulário específico, não encontrei uma resposta para a questão da convergência. Alguém poderia esclarecer isso?V0

Note-se que é o eigen-vector correspondente ao valor próprio uma de eu - D - 1 G . Com base na observação, chamamos o Teorema 2.1 de Autovalores de matrizes atualizadas de primeiro nível com algumas aplicações de Jiu Ding e Ai-Hui Zhou. V01ID1L

Teorema 2,1 Vamos e v ser dois n vectores de coluna -dimensional tal que u é um vector próprio de uma associada com valores próprios λ 1 . Então, os autovalores de A + u v T são { λ 1 + u T v , λ 2 , , λ n } contando a multiplicidade algébrica.uvnuAλ1A+uvT{λ1+uTv,λ2,,λn}

Em seguida, sabemos que os espectros de é o mesmo que eu - D - 1 L , excepto que o valor próprio 1 neste último é deslocado por - uma para o valor próprio zero no ex. Como ρ ( I - D - 1 L ) ( - 1 , 1 ] , temos ρ ( ˜ S ) ( - 1 , 1 ) .S~ID1L11ρ(ID1L)(1,1]ρ(S~)(1,1)

Hui Zhang
fonte
Obrigado por responder. Algo semelhante foi o que eu considerei: a saber, com o Laplaciano ponderado estruturado como o acima, pode ser mostrado que seus autovalores estão dentro de [ 0 , 2 ) , portanto, com o raio espectral dentro ( 0 , 2 ) (um valor próprio é maior que 0 e pelo menos um é 0 ). Portanto, o raio espectral da matriz de iteração I - D - 1 L é menor que 1 , portanto, com Jacobi convergente. Talvez a suposição acima no raio espectral deD1L[0,2)(0,2)00ID1L1 (excluindo 0 ) não é seguro? ID1L0
usero
Eu acho que o espectro de deve estar em [ 0 , 2 ] , que é fechado em 2 . Eu não sei como você pode obter 2 excluídos. Do meu ponto de vista, o (teorema do círculo de Gershgorin) [ en.wikipedia.org/wiki/Gershgorin_circle_theorem] só pode fornecer a estimativa, incluindo 2 . Se for esse o caso, a estimativa do raio espectral de I - D - 1 L é 1 com a igualdade alcançável com os vetores no núcleo de LD1L[0,2]222ID1L1L. Eu acho que a convergência que você quer é que no espaço ortogonal do complemento como observado na 'resposta' acima. V1
Hui Zhang
Você pode dar uma olhada no Lema 1.7 (v) de math.ucsd.edu/~fan/research/cb/ch1.pdf A matriz pode ser considerada como um Laplaciano ponderado em um gráfico completo, portanto, com 2 excluídos . Eu acho que é um argumento suficiente para a prova de convergência? ............ Sua abordagem requer outro pré / pós-processamento das iterações além da centralização c . Estou perguntando porque você introduziu V 0 E em relação aos espectros de I - D - 1 L - V 0 V 0 : dado que o raio espectral ( s r ) de ID1L2cV0ID1LV0V0sr é ( 0 , 1 ] , a adição de - 1ID1L(0,1] , produziria sr<1. Este argumento não é bom o suficiente? 1nsr<1
usero
Olá, obrigado por apontar para um bom livro. Mas descobri que não posso dar uma olhada rápida. Sobre o seu último argumento, ele aparece quase o mesmo que a "resposta" acima. Apenas tenha cuidado, você não está adicionando mas11n, por isso, não é uma simples adição aossrdeeu-D-1G. Geralmente, osrda soma de duas matrizesnãoéuma soma simples dossr's das matrizes individuais. 1n1n×nsrID1Lsrsr
Hui Zhang
Bom que você apontou isso. Sua abordagem requer outro pré / pós-processamento das iterações além da centralização c. Estou perguntando porque você introduziu o e pensei que estava falando em projetar o espaço nulo. Em caso afirmativo, a projeção é o espaço nulo realmente necessário para a convergência? V0
usero
5

Os métodos de Krylov nunca usam explicitamente a dimensionalidade do espaço em que iteram; portanto, você pode executá-los em sistemas singulares, desde que mantenha as iterações no subespaço não nulo. Isso normalmente é feito projetando o espaço nulo em cada iteração. Há duas coisas que podem dar errado: a primeira é muito mais comum que a segunda.

  1. O pré-condicionador é instável quando aplicado ao operador singular. Os solucionadores diretos e a fatoração incompleta podem ter essa propriedade. Por uma questão prática, apenas escolhemos diferentes pré-condicionadores, mas existem maneiras mais básicas de projetar pré-condicionadores para sistemas singulares, por exemplo, Zhang (2010) .
  2. Em alguma iteração, está no subespaço não nulo, mas A x vive inteiramente no espaço nulo. Isso só é possível com matrizes não simétricas. O GMRES não modificado divide-se nesse cenário, mas consulte Reichel e Ye (2005) para obter variantes livres de discriminação.xAx

Para resolver sistemas singulares usando o PETSc, consulte KSPSetNullSpace(). A maioria dos métodos e pré-condicionadores pode resolver sistemas singulares. Na prática, o pequeno espaço nulo para PDEs com condições de contorno de Neumann quase nunca é um problema, desde que você informe o solucionador de Krylov do espaço nulo e escolha um pré-condicionador razoável.

Pelos comentários, parece que você está especificamente interessado em Jacobi. (Por que? Jacobi é útil como um multigrid suave, mas existem métodos muito melhores para usar como solucionadores.) Jacobi aplicado a não converge quando o vetor b tem um componente no espaço nulo de A , no entanto, o parte da solução ortogonal ao espaço nulo converge; portanto, se você projetar o espaço nulo de cada iteração, ele convergirá. Como alternativa, se você escolher um bec consistente e uma estimativa inicial, as iterações (na aritmética exata) não acumulam componentes no espaço nulo.Ax=bbAb

Jed Brown
fonte
QA1=QTAQA1A1Z(IZ)P1Ax=(IZ)P1b
Hmm, parece que eu respondi a um comentário que foi excluído. Deixarei o comentário aqui, caso seja útil.
Jed Brown
diag(A)
Xk+1=D1(b(AD)Xk)Xk+1Xk?
usero
1
In a reasonable basis, Jacobi should be stable. It can also use 1 on the diagonal if the diagonal matrix element is 0, the projection still removes the null space. Are you planning to use a Krylov method like CG or GMRES? If not, why not? If you are, you just need an orthogonal basis for the null space. You only have the constant mode in your null space, so an orthogonal projector into the null space is N=ZZT where Z is the column vector. The orthogonal projector that removes the null space is thus IN. (My first comment had a mistake, if Z is the basis, N=IZZT is the projector.)
Jed Brown