Entendendo como o Numpy faz SVD

13

Eu tenho usado métodos diferentes para calcular a classificação de uma matriz e a solução de um sistema matricial de equações. Me deparei com a função linalg.svd. Comparando isso com o meu próprio esforço de resolver o sistema com a Eliminação Gaussiana, parece ser mais rápido e preciso. Estou tentando entender como isso é possível.

Até onde eu sei, a função linalg.svd usa um algoritmo QR para calcular os autovalores da minha matriz. Sei como isso funciona matematicamente, mas não sei como o Numpy consegue fazê-lo tão rapidamente e sem perder muita precisão.

Então, minha pergunta: como funciona a função numpy.svd e, mais especificamente, como ela é executada com rapidez e precisão (em comparação com a eliminação gaussiana)?

RobVerheyen
fonte
2
numpy usa a rotina Lapack dgesddpara SVDs com valor real. Portanto, sua pergunta real é provavelmente "como funciona o Lapack dgesdd?", E isso é bastante problemático para o stackoverflow.
talonmies
Se você está realmente curioso, sugiro examinar a fonte do LAPACK.
Obrigado por seus comentários, e minhas desculpas, sou offtopic.
RobVerheyen
Esta postagem é uma publicação cruzada do Stack Overflow . A postagem cruzada geralmente não é recomendada nos sites do Stack Exchange. O protocolo padrão para reposicionar uma pergunta em um site diferente é fechar, excluir ou migrar a postagem original antes de tentar republicar em um site diferente. (Se você migrar a questão, é anunciado automaticamente.)
Geoff Oxberry
Me desculpe, eu não estava ciente do protocolo. Espero ainda conseguir uma resposta.
RobVerheyen

Respostas:

15

Há vários problemas na sua pergunta.

Não use a eliminação gaussiana (fatoração LU) para calcular a classificação numérica de uma matriz. A fatoração da LU não é confiável para esse fim na aritmética de ponto flutuante. Em vez disso, use uma decomposição QR que revele a classificação (como xGEQPXou xGEPQYno LAPACK, onde x é C, D, S ou Z, embora essas rotinas sejam difíceis de rastrear; consulte a resposta de JedBrown em uma pergunta relacionada ) ou use um SVD (decomposição de valor singular, como xGESDDor xGESVD, onde x é novamente C, D, S ou Z). O SVD é um algoritmo mais preciso e confiável para a determinação da classificação numérica, mas requer mais operações de ponto flutuante.

No entanto, para resolver um sistema linear, a fatoração de LU (com rotação parcial, que é a implementação padrão no LAPACK) é extremamente confiável na prática. Existem alguns casos patológicos para os quais a fatoração de LU com pivotamento parcial é instável (ver Lição 22 na Álgebra Linear Numérica).por Trefethen e Bau para detalhes). A fatoração QR é um algoritmo numérico mais estável para resolver sistemas lineares, e é provavelmente por isso que ele fornece resultados tão precisos. No entanto, requer mais operações de ponto flutuante do que a fatoração LU por um fator 2 para matrizes quadradas (acredito; JackPoulson pode me corrigir nisso). Para sistemas retangulares, a fatoração QR é uma escolha melhor, pois produzirá soluções de mínimos quadrados para sistemas lineares sobredeterminados. O SVD também pode ser usado para resolver sistemas lineares, mas será mais caro que a fatoração QR.

janneb está correto que numpy.linalg.svd é um invólucro xGESDDno LAPACK. As decomposições de valores singulares prosseguem em dois estágios. Primeiro, a matriz a ser decomposta é reduzida à forma bidiagonal. O algoritmo usado para reduzir a forma bidiagonal no LAPACK é provavelmente o algoritmo Lawson-Hanson-Chan e usa a fatoração QR em um ponto. A aula 31 de Álgebra Linear Numérica, de Trefethen e Bau, fornece uma visão geral desse processo. Em seguida, xGESDDusa um algoritmo de dividir e conquistar para calcular os valores singulares e os vetores singulares esquerdo e direito da matriz bidiagonal. Para obter informações sobre essa etapa, você precisará consultar Matrix Computations de Golub e Van Loan ou Applied Numerical Linear Algebra de Jim Demmel.

Finalmente, você não deve confundir valores singulares com valores próprios . Esses dois conjuntos de quantidades não são os mesmos. O SVD calcula os valores singulares de uma matriz. A computação numérica de Cleve Moler com o MATLAB fornece uma boa visão geral das diferenças entre valores singulares e valores próprios . Em geral, não há relação óbvia entre os valores singulares de uma matriz e seus valores próprios, exceto no caso de matrizes normais , em que os valores singulares são o valor absoluto dos valores próprios.

Geoff Oxberry
fonte
Eu acho que "não relacionado" é bastante forte para a relação entre autovalores e valores singulares. O relacionamento é bastante obscuro, a menos que você conheça a decomposição completa da sua matriz na Jordânia, mas você pode usar uma para obter estimativas da outra se tiver informações (ou estiver disposto a fazer suposições) sobre a referida decomposição na Jordânia.
Dan
O que você sugeriria?
amigos estão dizendo
Primeiro de tudo, obrigado pela resposta elaborada. Descobri que não posso usar a decomposição da LU para determinar a classificação da matriz da maneira mais difícil. Sua resposta parece sugerir que a fatoração QR seria de fato um método mais rápido de resolver meu problema, correto? Existe uma vantagem distinta no uso de SVD? Eu estava ciente do fato de que valores singulares não são autovalores. Eu estava me referindo ao fato de que valores singulares podem ser calculados como autovalores da matriz multiplicados por sua transposição a partir da esquerda. Me desculpe, isso não estava claro.
RobVerheyen
Devo acrescentar que a matriz que estou resolvendo é realmente singular. De fato, a classificação da matriz é apenas cerca da metade do tamanho da matriz. Talvez isso torne algum método mais preferível?
RobVerheyen
1
@RobVerheyen: o QR será mais lento que o LU, mas consideravelmente mais preciso. O SVD será ainda mais lento que o QR, mas o SVD é considerado o método mais confiável para determinar a classificação numérica (por exemplo, o MATLAB usa o SVD em sua rankfunção). Há também um pouco de discrição ao usar qualquer uma das abordagens; na abordagem SVD, a classificação numérica é o número de valores singulares acima de um ponto de corte especificado (geralmente muito pequeno). (A abordagem QR é semelhante, mas substitui valores singulares com elementos da diagonal da matriz de R.)
Geoff Oxberry
8

Devido à redação da sua pergunta, estou assumindo que sua matriz é quadrada. As rotinas SVD do LAPACK, como zgesvd , prosseguem essencialmente em três estágios para matrizes quadradas:

  1. vocêUMAVUMAUMAB: =vocêUMAHUMAVUMAvocêUMAVUMABO(n3)
  2. {vocêB,VB,Σ}B=vocêBΣVBHO(n2)O(n3)
  3. vocêUMABVUMAH=UMAUMA=(vocêUMAvocêB)Σ(VUMAVB)HvocêUMAVUMAvocêBVBO(n3)
Jack Poulson
fonte
7

numpy.linalg.svd é um wrapper em torno do {Z, D} GESDD do LAPACK. O LAPACK, por sua vez, é cuidadosamente escrito por alguns dos principais especialistas do mundo em álgebra linear numérica. De fato, seria muito surpreendente se alguém não familiarizado com o campo conseguisse vencer o LAPACK (em velocidade ou precisão).

Quanto ao motivo pelo qual o QR é melhor que a eliminação gaussiana, isso provavelmente é mais apropriado para /scicomp//

janneb
fonte
Obrigado pela resposta e pela referência. Vou tentar por lá.
RobVerheyen