Como o operador de barra invertida do MATLAB resolve

36

Eu estava comparando alguns dos meus códigos com os códigos "estoque" do MATLAB. Estou surpreso com os resultados.

Corri um código de exemplo (Sparse Matrix)

n = 5000;
a = diag(rand(n,1));
b = rand(n,1);
disp('For a\b');
tic;a\b;toc;
disp('For LU');
tic;LULU;toc;
disp('For Conj Grad');
tic;conjgrad(a,b,1e-8);toc;
disp('Inv(A)*B');
tic;inv(a)*b;toc;

Resultados :

    For a\b
    Elapsed time is 0.052838 seconds.

    For LU
    Elapsed time is 7.441331 seconds.

    For Conj Grad
    Elapsed time is 3.819182 seconds.

    Inv(A)*B
    Elapsed time is 38.511110 seconds.

Para matriz densa:

n = 2000;
a = rand(n,n);
b = rand(n,1);
disp('For a\b');
tic;a\b;toc;
disp('For LU');
tic;LULU;toc;
disp('For Conj Grad');
tic;conjgrad(a,b,1e-8);toc;
disp('For INV(A)*B');
tic;inv(a)*b;toc;

Resultados:

For a\b
Elapsed time is 0.575926 seconds.

For LU
Elapsed time is 0.654287 seconds.

For Conj Grad
Elapsed time is 9.875896 seconds.

Inv(A)*B
Elapsed time is 1.648074 seconds.

Como diabos é tão impressionante?

Inquérito
fonte
1
A barra invertida interna do MATLAB, em outras palavras, o solucionador direto de um sistema de equações lineares, usa o método Multifrontal para matriz esparsa, e é por isso que A \ B é tão impressionante.
Shuhao Cao 25/01/12
1
Ele usa o código de Tim Davis disponível em cise.ufl.edu/research/sparse . Também a grandiosidade vai embora quando você tem um problema não-trivial
stali
1
O que é "LULU"? Por que você acha que é uma boa implementação de uma fatoração de LU e subsequente resolução direta?
Jed Brown
3
@Nunoxic Qual implementação? Você mesmo escreveu? Álgebra linear densa de alto desempenho, embora geralmente seja bem compreendida por algoritmos, não é fácil de implementar com eficiência no hardware moderno. As melhores implementações BLAS / Lapack devem se aproximar do pico de uma matriz desse tamanho. Além disso, pelos seus comentários, estou tendo a impressão de que você acha que LU e Eliminação Gaussiana são algoritmos diferentes.
Jed Brown
1
Ele chama um código Fortran escrito usando o Intel MKL.
Inquérito

Respostas:

37

No Matlab, o comando '\' invoca um algoritmo que depende da estrutura da matriz A e inclui verificações (pequenas despesas gerais) nas propriedades de A.

  1. Se A for escasso e com faixas, use um solucionador de faixas.
  2. Se A for uma matriz triangular superior ou inferior, utilize um algoritmo de substituição para trás.
  3. Se A é simétrico e possui elementos diagonais positivos reais, tente uma fatoração de Cholesky. Se A for escasso, use primeiro a reordenação para minimizar o preenchimento.
  4. Se nenhum dos critérios acima for atendido, faça uma fatoração triangular geral usando a eliminação gaussiana com pivô parcial.
  5. Se A for escasso, empregue a biblioteca UMFPACK.
  6. Se A não for quadrado, utilize algoritmos baseados na fatoração QR para sistemas indeterminados.

Para reduzir a sobrecarga, é possível usar o comando linsolve no Matlab e selecionar um solucionador adequado entre essas opções.

Allan P. Engsig-Karup
fonte
Supondo que estou lidando com uma matriz densa não estruturada de 10000x10000 com todos os elementos diferentes de zero (alto nível de densidade), qual seria minha melhor aposta? Eu quero isolar esse algoritmo 1 que funciona para matrizes densas. É LU, QR ou Eliminação Gaussiana?
Inquérito
1
Soa como um passo 4, em que a eliminação gaussiana é invocada, o que corresponde ao caso mais geral em que nenhuma estrutura de A pode ser explorada para aumentar o desempenho. Então, basicamente, essa é uma fatoração da LU e a subsequente subsequente seguida por uma etapa de substituição para trás.
Allan P. Engsig-Karup
Obrigado! Eu acho que isso me dá uma direção para pensar. Atualmente, a eliminação gaussiana é o melhor que temos para resolver problemas não estruturados, está correto?
Inquérito
37

Se você quiser ver o que a\bfaz para sua matriz em particular, pode definir spparms('spumoni',1)e descobrir exatamente qual algoritmo ficou impressionado. Por exemplo:

spparms('spumoni',1);
A = delsq(numgrid('B',256));
b = rand(size(A,2),1);
mldivide(A,b);  % another way to write A\b

irá produzir

sp\: bandwidth = 254+1+254.
sp\: is A diagonal? no.
sp\: is band density (0.01) > bandden (0.50) to try banded solver? no.
sp\: is A triangular? no.
sp\: is A morally triangular? no.
sp\: is A a candidate for Cholesky (symmetric, real positive diagonal)? yes.
sp\: is CHOLMOD's symbolic Cholesky factorization (with automatic reordering) successful? yes.
sp\: is CHOLMOD's numeric Cholesky factorization successful? yes.
sp\: is CHOLMOD's triangular solve successful? yes.

para que eu possa ver que "\" acabou usando "CHOLMOD" nesse caso.

dranxo
fonte
3
+1 para novas configurações do MATLAB de que nunca ouvi falar. Bem jogado, senhor.
Geoff Oxberry
2
Ei, obrigado! Está dentro help mldivide.
dranxo
16

Para matrizes esparsas, o Matlab usa UMFPACK para a \operação " ", que, no seu exemplo, basicamente percorre os valores de a, inverte-os e multiplica-os pelos valores de b. Para este exemplo, no entanto, você deve usar b./diag(a), o que é muito mais rápido.

Para sistemas densos, o operador de barra invertida é um pouco mais complicado. Uma breve descrição do que é feito quando é fornecida aqui . De acordo com essa descrição, no seu exemplo, o Matlab resolveria a\busando a substituição reversa. Para matrizes quadradas gerais, é usada uma decomposição LU.

Pedro
fonte
Regd. Escassez, o inv de uma matriz diag seria simplesmente o recíproco dos elementos diagonais, para que b./diag(a) funcionasse, mas a \ b funciona de maneira impressionante para matrizes esparsas gerais também. Por que o linsolve ou LULU (minha versão otimizada da LU) não é mais rápido que a \ b para matrizes densas nesse caso.
Inquérito
@ Não-tóxico O seu LULU faz alguma coisa para detectar a diagonalidade ou triangularidade da matriz densa? Caso contrário, levará tanto tempo com todas as matrizes, independentemente de seu conteúdo ou estrutura.
Pedro
Um pouco. Mas, com os sinalizadores linsolve OPT, defini tudo o que há para definir sobre a estrutura. No entanto, a \ b é mais rápido.
Inquérito
@ Não tóxico, toda vez que você chama uma função de usuário, você causa sobrecarga. O Matlab faz tudo para a barra invertida internamente, por exemplo, a decomposição e a aplicação subsequente do lado direito, com muito pouca sobrecarga, e, portanto, será mais rápido. Além disso, em seus testes, você deve usar mais de uma chamada para obter intervalos confiáveis, por exemplo tic; for k=1:100, a\b; end; toc.
Pedro
5

Como regra geral, se você tem uma matriz esparsa de complexidade razoável (ou seja, não precisa ser o estêncil de 5 pontos, mas pode de fato ser uma discretização das equações de Stokes para as quais o número de não zeros por linha é muito maior que 5), um solucionador direto esparso, como o UMFPACK, geralmente vence um solucionador iterativo de Krylov se o problema não for maior que cerca de 100.000 incógnitas.

Em outras palavras, para a maioria das matrizes esparsas resultantes de discretizações 2D, um solucionador direto é a alternativa mais rápida. Somente para problemas em 3D nos quais você rapidamente ultrapassa as 100.000 incógnitas, torna-se imperativo usar um solucionador iterativo.

Wolfgang Bangerth
fonte
3
Não está claro para mim como isso responde à pergunta, mas também discuto a premissa. É verdade que os solucionadores diretos tendem a funcionar bem em problemas 2D de tamanho modesto, mas os solucionadores diretos tendem a sofrer em 3D muito antes de 100 mil incógnitas (os separadores de vértices são muito maiores que em 2D). Além disso, afirmo que, na maioria dos casos (por exemplo, operadores elípticos), os solucionadores iterativos podem ser feitos mais rapidamente que os solucionadores diretos, mesmo para problemas 2D de tamanho moderado, mas pode ser necessário um esforço significativo para fazê-lo (por exemplo, usando os ingredientes certos para fazer o desempenho de multigrid) .
Jed Brown
1
Se seus problemas são razoavelmente pequenos e você não quer pensar em projetar solucionadores implícitos, ou se seus problemas são muito difíceis (como o Maxwell de alta frequência) e você não deseja dedicar sua carreira a projetar bons solucionadores, então eu concorda que os solucionadores diretos esparsos são uma ótima opção.
Jed Brown
Na verdade, meu problema é criar um bom solucionador denso. Não tenho uma aplicação específica como tal (portanto, nenhuma estrutura). Eu queria ajustar alguns dos meus códigos e torná-los competitivos com o que é usado atualmente. Eu só estava curioso sobre a \ b e como ele sempre supera a porcaria do meu código.
Inquérito
@JedBrown: Sim, talvez com um esforço significativo seja possível conseguir que os solucionadores iterativos superem o solucionador direto, mesmo em pequenos problemas 2D. Mas por que fazer isso? Para problemas com <100k desconhecidos, os solucionadores diretos esparsos em 2d são bastante rápidos. O que eu queria dizer é: para problemas tão pequenos, não gaste seu tempo mexendo e descobrindo a melhor combinação de parâmetros - basta pegar a caixa preta.
Wolfgang Bangerth
Já existe uma diferença de ordem de magnitude no tempo de execução entre um escasso multigrid de Cholesky e geométrico (baseado em matriz) para problemas 2D "fáceis" com 100k dofs usando um estêncil de 5 pontos (~ 0,5 segundo em comparação com ~ 0,05 segundos). Se o seu estêncil usa segundos vizinhos (por exemplo, discretização de quarta ordem, Newton para algumas opções de reologia não linear, estabilização etc.), o tamanho dos separadores de vértices praticamente dobra, de modo que o custo da solução direta sobe cerca de 8x (o custo é mais dependente do problema para MG). Com muitas etapas de tempo ou loops de otimização / UQ, essas diferenças são significativas.
Jed Brown