Como reordenar variáveis ​​para produzir uma matriz em faixas de largura de banda mínima?

15

Estou tentando resolver uma equação de Poisson 2D por diferenças finitas. No processo, obtenho uma matriz esparsa com apenas variáveis ​​em cada equação. Por exemplo, se as variáveis ​​fossem U , a discretização produziria:5U

Ui1,j+Ui+1,j4Ui,j+Ui,j1+Ui,j+1=fi,j

Eu sei que posso resolver esse sistema por um método iterativo, mas ocorreu-me que, se eu ordenasse as variáveis ​​adequadamente, seria possível obter uma matriz em faixas que poderia ser resolvida por um método direto (por exemplo, eliminação gaussiana w / o pivotante). Isso é possível? Existem estratégias para fazer isso em outros sistemas esparsos, talvez menos estruturados?

Paulo
fonte
2
Algo como Cuthill-McKee, então?
JM
Interessante ... nunca ouvi falar do algoritmo Cuthill-McKee antes! :)
Paul
11
Há também um Cuthill-McKee Reverso também.
precisa saber é o seguinte
11
Espero que fique claro nas respostas, mas você não deseja usar um solucionador de faixas para esse problema, nem escolher uma ordem que minimize a largura de banda. Talvez a pergunta ou a resposta escolhida possa ser editada para deixar isso claro, caso contrário, temo que esse mito seja perpetuado. Fiz uma comparação visual e comparei o preenchimento em scicomp.stackexchange.com/a/880/119 .
perfil completo de Jed Brown
@JedBrown: Na verdade, não estou trabalhando muito com um problema de poisson, por si só ... Meu problema tem uma estrutura semelhante ao problema de poisson ... As indicações das variáveis ​​(i e j) são exatamente as mesmas, e a matriz é diagonalmente dominante com as entradas fora da diagonal (dentro da mesma linha) adicionadas exatamente à soma da entrada na diagonal.
Paul

Respostas:

13

Esse é um problema bem estudado no campo dos solucionadores diretos esparsos. Eu recomendo a leitura da visão geral de Joseph Liu do método multifrontal para ter uma idéia melhor de como os novos pedidos e supernós afetam o tempo de preenchimento e solução.

A dissecção aninhada é uma maneira extremamente comum de gerar a reordenação e consiste essencialmente em particionamento de gráfico recursivo. O MeTiS é o padrão de fato para particionamento de gráficos, e você pode ler sobre algumas das idéias por trás dele aqui . Outro pacote comumente usado é o SCOTCH , e o Chaco também é importante, pois seus autores introduziram o particionamento de gráficos em vários níveis , que também é a idéia fundamental por trás do MeTiS .

George e Liu mostrou em seu clássico livro que as soluções esparsas-direta 2d requerem apenas trabalho e O ( n log n ) de memória, enquanto 3d esparsa-direta requer O ( n 2 ) trabalho e O ( n 4 / 3 ) memória.O(n3/2)O(nlogn)O(n2)O(n4/3)

Jack Poulson
fonte
Você tem uma citação para a referência de George e Liu?
Paul
Adicionado; Eu estava prestes a sair do carro quando o enviei pela primeira vez. Sei que existe uma versão disponível gratuitamente do livro on-line em algum lugar (Jed sabe onde está), mas não consegui encontrá-lo.
Jack Poulson
Atualizei o link para apontar para o PDF do livro em vez da resenha.
Jed Brown
@JedBrown Essa foi uma ótima referência! Muito obrigado! :)
Paul
11
@Alexander Todo mundo atribui o 3D vinculado a George e Liu, embora eu não saiba se eles apontam explicitamente no livro. É óbvio da teoria no entanto. O separador de vértice mínimo de grade é n 2 / 3 = m x m . A matriz densa associado com que supernó tem ( n 2 / 3 ) 2 = N 4 / 3 entradas e requer ( n 2 / 3 ) 3 = n 2n=m×m×mn2/3=m×m(n2/3)2=n4/3(n2/3)3=n2operações a fatorar. O termo logarítmico no caso 2D é mais sutil e é tratado no Capítulo 8 sobre Dissecação aninhada, que atinge o limite inferior.
Jed Brown
5

Cuthill-McKee é o padrão de fato para o que você deseja fazer. Se você quiser jogar com esse método, há uma implementação fácil de usar do algoritmo (e seu reverso) na Boost Graph Library (BGL), e a documentação contém exemplos de como usá-lo.

Wolfgang Bangerth
fonte
na verdade, inverta Cuhill-McKee; geralmente dá menos preenchimento. Mas uma ordem de dissecção aninhada é muito superior a uma ordem de baixa largura de banda.
Arnold Neumaier 29/03/2012
4

Falando em métodos multifrontais, Tim Davis , que trabalha em métodos multifrontais para fatoração de LU ( UMFPACK ), possui várias rotinas que reorganizam matrizes para minimizar o preenchimento. Você pode encontrá-los aqui como parte do SuiteSparse . O SuiteSparse usa o MeTiS.

Outra coisa a ser observada: em alguns problemas, você pode ser inteligente ao ordenar variáveis, para obter padrões unidos ou quase unidos, o que pode poupar o problema (e o tempo de CPU) de chamar esses algoritmos. No entanto, essa reordenação inteligente requer insights de sua parte e não é nem de longe tão geral quanto os algoritmos de reordenação baseados na teoria dos grafos que as pessoas mencionaram em suas respostas aqui.

Geoff Oxberry
fonte
De nada, Paul. Se você gosta, vote.
precisa saber é o seguinte
3

Existe um algoritmo chamado ADI (Alternating Direction Implicit) nos círculos matemáticos aplicados e o operador Split nos círculos da física que faz basicamente o que você descreve. É um método iterativo e segue este procedimento básico:

  1. yx

  2. xy

  3. Repita 1 e 2 até que o erro seja tão pequeno quanto você deseja.

Não conheço a complexidade formal desse algoritmo, mas descobri que ele converge em menos iterações do que coisas como Jacobi e Gauss-Seidel toda vez que o uso.

Dan
fonte
2
Se você decidir seguir a rota de divisão do operador, algo que você precisará ter cuidado é que a divisão do operador é conhecida por resultar em erros em soluções de estado estacionário em alguns casos. (Um dos meus colegas de laboratório desenvolveu uma maneira de superar essa dificuldade, mas eu não acredito que ele a tenha publicado ainda.) Além disso, a divisão do operador é conhecida por causar erros numéricos. Existem maneiras bem estabelecidas de estimar esses erros a posteriori ; Don Estep fez um excelente trabalho nessa área.
precisa saber é o seguinte
@ GeoffOxberry Parece que você está se referindo a uma divisão diferente. Você pode usar o ADI em um esquema totalmente implícito que não possui erro de divisão, porque na verdade resolve o sistema. Existem também métodos IMEX que controlam rigorosamente os erros de divisão.
perfil completo de Jed Brown
xy
Nunca ouvi falar de Godunov e Strang se separarem. Costumo dividir meu operador com a fórmula de Baker-Campbell-Hausdorf. Aquilo é a mesma coisa?
Dan