Como o SVD de uma matriz é calculado na prática

Como o MATLAB, por exemplo, calcula o SVD de uma determinada matriz? Suponho que a resposta provavelmente envolva calcular os autovetores e autovalores de A*A'. Se for esse o caso, eu também gostaria de saber como isso os computa.

$AA^{T}$

Vejo

GH Golub e C. Reinsch. Soluções de Decomposição de Valor Singular e Mínimos Quadrados. Numerische Mathematik 14: 403-420, 1970.

Este material é discutido em muitos livros sobre álgebra linear numérica.

Além do artigo (agora clássico) de Golub-Reinsch, Brian observa em sua resposta (vinculei à versão Handbook do artigo), bem como do artigo (também agora clássico) predecessor de Golub-Kahan , houve vários importantes desenvolvimentos na computação do SVD desde então. Primeiro, tenho que resumir como o método usual funciona.

A idéia de calcular o SVD de uma matriz é qualitativamente semelhante ao método usado para calcular a composição automática de uma matriz simétrica (e, como observado no PO, há uma relação íntima entre eles). Em particular, procede-se em dois estágios: a transformação em uma matriz bidiagonal e, em seguida, localizando o SVD de uma matriz bidiagonal. Isso é completamente análogo ao procedimento de primeiro reduzir uma matriz simétrica à forma tridiagonal e depois calcular a composição automática do tridiagonal resultante.

Para calcular o SVD de uma matriz bidiagonal, uma descoberta particularmente interessante foi o artigo de Jim Demmel e Velvel Kahan , que demonstrou que é possível computar até minúsculos valores singulares de uma matriz bidiagonal com boa precisão, modificando adequadamente o método proposto inicialmente em Golub-Reinsch. Este foi, em seguida, seguido pela (re?) Descoberta do DQD algoritmo , que é uma descendente da idade algoritmo quociente-diferença de Rutishauser. (Beresford Parlett faz uma boa discussão aqui.) Se a memória servir, agora é o método preferido usado internamente pelo LAPACK. Além disso, sempre foi possível derivar versões SVD de desenvolvimentos na solução de autovalores simétricos; por exemplo, existe uma versão SVD de dividir e conquistar, bem como uma versão SVD do antigo algoritmo Jacobi (que pode ser mais preciso em algumas circunstâncias).

Quanto à bidiagonalização, um método aprimorado foi descrito no artigo de Barlow , que requer um pouco mais de trabalho do que o procedimento original de Golub e Reincsh, mas produz matrizes bidiagonais mais precisas.