Qual é a função LAPACK correspondente por trás do Matlab [Q, R, E] = qr (A)?

12

Eu atualmente tentando calcular de forma barata uma estimativa boa classificação para uma matriz . Portanto, eu calculo uma decomposição QR dinâmica de colunas usandoA

[Q,R,E]=qr(A)

no Matlab. Estimo a classificação de usandoA

tol = size(A,n)*eps*norm(A,'fro'); 
r = sum(abs(diag(R))>tol)

Isso funciona bem e um gráfico sobre todas as entradas diagonais de R se parece com: plot (classificar (abs (diag (R)), 1, 'descendente'), 'r *')

Se for portar todo o algoritmo para C / Fortran, substituo [Q, R, E] = qr (A) usando DGEQP3 do LAPACK, que também calcula uma coluna que gira a decomposição QR. Mas se eu usar a mesma estimativa para a classificação, geralmente entendo algo errado. O mesmo enredo para o produzido a partir do DGEQP3 parece R

A matriz de entrada é exatamente a mesma para os dois experimentos.

Minha pergunta agora é em qual função LAPACK a coluna que gira a decomposição QR do Matlab se baseia?

Obrigado por qualquer ajuda, Grisu

Edit: DGEQPF fornece o mesmo resultado errado.

Edit2:

Edit3: - Usando o GDB descobri, o Matlab 2010b chama DGEQP3: # 3 0xaa46ce2f em dgeqp3_ () em /usr/ubuntu10.04/matlabr2010b/bin/glnx86/../../bin/glnx86/../. ./bin/glnx86/mllapack.so Por que obtenho o resultado errado usando o LAPACK (3.4.0 inclui as correções mencionadas na Nota de trabalho 176)?

MK tcp Grisu
fonte
Você pode provocar o mesmo comportamento com uma matriz menor que possa compartilhar aqui?
GertVdE
Does tem nenhuma estrutura especial? Quero dizer que o MATLAB usa UMFPACK para álgebra linear esparsa, mas não tenho certeza de quais são as bibliotecas de álgebra linear densa subjacentes. A
Aron Ahmadia 27/03/12
é escassa? Você pode definir ">> spparms ('spumoni', 1)" e ver que algo chamado "SuiteSparseQR" seja usado pelo matlab nesse caso. A
dranxo
Grisu - Eu adoraria ver sua matriz. No entanto, o link www-e.uni-magdeburg.de/makoehle/A.mtx.gz não está mais ativo (no momento, de qualquer maneira). Você tem um link atual para a matriz? Obrigado, Les Foster
@LeslieFoster - bem-vindo ao scicomp!
Aron Ahmadia 06/07/12

Respostas:

7

Há duas questões em mãos aqui:

  • é densa ou esparsa?A
  • Você tem a mesma pilha de software que as bibliotecas internas do MATLAB?

Denso ou escasso?

O MATLAB não menciona mais explicitamente as rotinas LAPACK que ele chama para obter uma fatoração QR se for denso. Se as informações na documentação do MATLAB R2008b também se mantiverem para versões posteriores, o MATLAB chamará do LAPACK quando você ligar . Se for escasso, o MATLAB chamará o SuiteSparseQR , fora do grupo de Tim Davis, que é fornecido com o UMFPACK na biblioteca SuiteSparse.AADGEQP3[Q,R,E] = qr(A)A

Você tem a mesma pilha de software que as bibliotecas internas do MATLAB?

Provavelmente não, o que pode ser uma das razões pelas quais você está obtendo resultados diferentes.

Eu me deparei com esse problema ao testar uma unidade que estava escrevendo que usava fatorações QR. Eu usei o MATLAB para prototipar meu trabalho e obtive resultados diferentes do que usar o LAPACK ou o NumPy. Pelo que sei, porque o MathWorks não facilita a localização dessas informações, o MATLAB usa uma versão do LAPACK não anterior à versão 3.1.1 e a biblioteca MKL BLAS da Intel (para Windows, Intel Mac e Linux) versão 9.1 ou superior (veja aqui ). Não consegui encontrar nada sobre a versão do SuiteSparse que o MATLAB usa. Ao pesquisar on-line ou consultar os arquivos da biblioteca do seu sistema, você poderá obter informações adicionais. Você pode tentar alterar as bibliotecas às quais o MATLAB se vincula para poder comparar com as mesmas bibliotecas nos pacotes de software; Eric Chu fornece um ótimo artigoisso mostra pelo menos como você pode substituir a biblioteca BLAS do MATLAB pela sua (é claro, você faz isso por seu próprio risco). Ele sugere que você também possa fazer o mesmo com o LAPACK. Pode até ser possível substituir a versão do SuiteSparse que o MATLAB usa por sua própria versão.

As versões do LAPACK são alteradas, assim como as versões do BLAS; eles podem usar algoritmos diferentes de versão para versão ou convenções de pedidos diferentes, mesmo que com ortogonal e seja triangular superior, independentemente da versão. Essas mudanças tornam a reprodutibilidade um desafio. Mesmo a tolerância que você usa para a determinação da classificação numérica é uma decisão judicial; você parece estar usando uma tolerância padrão.Q RA=QRQR

Acabei usando o NumPy para prototipar meus resultados para a fatoração QR, porque ele usa as bibliotecas do sistema BLAS e LAPACK. O NumPy e o SciPy não substituem o MATLAB, porque as duas bibliotecas combinadas carecem de algumas funcionalidades do MATLAB, mas para essa tarefa de álgebra linear específica, Python + NumPy + SciPy + Matplotlib deve funcionar bem.

Geoff Oxberry
fonte
Conseguir o mesmo empilhamento de software que o Matlab é impossível, eu acho. O uso de outro ambiente para prototipagem também não é desejado. O problema é o código funciona em Matlab corretamente, mas não em C.
MK aka Grisu
@Grisu: Eu acho que seria muito difícil obter a mesma pilha de software, sem tentar vincular suas bibliotecas. O que me deixa confuso é como você sabe que o resultado no MATLAB está correto e o resultado em C está errado. Isso é algum tipo de matriz de teste que possui propriedades conhecidas? Mais precisamente, AronAhmadia está certo; além de replicar a arquitetura e a pilha de software, você não pode esperar obter os mesmos resultados com um algoritmo instável. Disseram-me basicamente a mesma coisa nos fóruns do MATLAB há dois anos.
Geoff Oxberry 28/03/2012
Na minha opinião, o QR não é instável. Não posso verificar diretamente qual decomposição QR está correta, mas, a partir da classificação e da matriz Q, calculo um projetor e lá posso facilmente verificar se um resultado bom ou ruim e o da Matlab são bons. Mas tento vincular-me às bibliotecas do Matlab.
MK aka Grisu
@Grisu: Há uma diferença distinta entre bons resultados e resultados corretos. Recentemente, implementei um cálculo incorreto que fez meus resultados parecerem maravilhosos. No entanto, o cálculo estava errado e o cálculo correto fez com que meus resultados parecessem menos impressionantes (mas, felizmente, ilustra que meus resultados estão corretos). Você está tentando calcular um projetor ortogonal ou um projetor oblíquo? (Eu pergunto porque partes significativas da minha tese estão em projetores oblíquos.)
Geoff Oxberry
1
@GeoffOxberry: fwiw, na minha versão de MATLAB, eu posso ligar internal.matlab.language.versionPlugins.blase internal.matlab.language.versionPlugins.lapackobter versões Blas e LAPACK
Amro
6

Veja a página de Leslie Foster no software revelador de rank . Veja também esta Nota de Trabalho do LAPACK, analisando falhas do QR que revela a classificaçãoxGEQP3 .

Você deve descobrir quais rotinas o MATLAB usa definindo pontos de interrupção em um depurador e examinando a pilha. A última vez que observei, há vários anos, o MATLAB usava bibliotecas compartilhadas; nesse caso, os nomes dos símbolos não podem ser removidos; portanto, você verá os nomes das funções na pilha de chamadas (mas não nos argumentos, porque definitivamente não mantém as informações de depuração).

Jed Brown
fonte
A página do software revelador de classificação não ajudou. O procedimento RRQR descrito foi a primeira coisa que uso na minha ideia, mas fornece resultados ainda piores do que a idéia de giro da coluna.
MK aka Grisu
2
@Grisu - Deveria ter ajudado você. O xGEQP3algoritmo não é completamente seguro para revelar a classificação. Se você deseja garantir que obtém o resultado certo, use o SVD ou um QR mais seguro como xGEQPXou xGEQPY. Você não pode esperar que um algoritmo instável retorne o mesmo resultado em arquiteturas diferentes ou em implementações diferentes (o MATLAB provavelmente está usando um LAPACK mais antigo).
Aron Ahmadia 28/03/12
Eu sei que o GEQP3 não é revelador de classificação, mas fornece resultados mais corretos do que as sub-rotinas RRQR. Usar um SVD é muito caro no meu algoritmo externo. Também vou falar com um dos autores do LAWN-176 e ele acha que esse erro não está coberto pelo bug. Eu também tentei o DGEQPF / DGEQP3 do LAPACK 3.0.0 com os mesmos resultados.
MK aka Grisu