Implementações eficientes em memória de decomposições parciais de valores singulares (SVD)

10

Para redução do modelo, quero calcular os vetores singulares esquerdos associados aos - digamos 20 - maiores valores singulares de uma matriz , em que N 10 6 e k 10 3 . Infelizmente, minha matriz A será densa sem nenhuma estrutura.UMARN,kN106k103UMA

Se eu apenas chamar a svdrotina do numpy.linalgmódulo em Python para uma matriz aleatória desse tamanho, ocorrerá um erro de memória. Isto é devido à atribuição de para a decomposição A = V S L .VRN,NUMA=VSvocê

Existem algoritmos por aí que evitam essa armadilha? Por exemplo, configurando apenas os vetores singulares associados a valores singulares diferentes de zero.

Estou pronto para negociar com precisão e tempo de computação.

Jan
fonte
11
Interessante, parece que Numpy não sabe como fazer um SVD fina ...
JM
Obrigado pela dica. De fato, numpy.linalg.svd tem a opção full_matricesque está definida como False, para que apenas as partes 'diferentes de zero' sejam computadas. No entanto, existe uma maneira de reduzir ainda mais a computação?
Jan
3
O numpyback-end usa o código fortran, a LAPACKE_dgesvdrotina para o svd padrão. No entanto, normalmente sua matriz é C_CONTIGOUS(verifique com matrix.flags). Portanto, ele copia os dados para o alinhamento do fortran. Além disso, durante a execução da rotina lapack dgesvd, é necessária outra cópia da sua matriz (ou pelo menos a memória). Você pode se livrar de uma cópia se certificar-se de que o alinhamento da memória seja do tipo fortran desde o início.
Bort

Respostas:

6

Se você deseja apenas alguns valores / vetores singulares, o ARPACK deve fazer o truque. Os documentos SVD não são ótimos e essa distribuição é mais atualizada.

Edição: Se você quiser fazer isso em python, SciPy tem um wrapper . Como sua matriz é densa, você pode tentar o formato BSR ( Block Sparse Row ).

Max Hutchinson
fonte
Vou dar uma olhada, como ARPACK integra com python ...
Jan
11
Parece que scipy tem invólucros. Vou adicioná-los para responder ao corpo.
Max Hutchinson
2

Talvez você possa tentar isso.

https://github.com/jakevdp/pypropack

Este é um wrapper Python para o pacote PROPACK, que implementa decomposições eficientes parciais de valor singular de matrizes esparsas grandes e operadores lineares.

Mass Zhou
fonte
2

O Intel MKL implementa o novo algoritmo Jacobi-SVD. Aqui estão os detalhes da implementação: http://www.netlib.org/lapack/lawnspdf/lawn169.pdf http://www.fernuni-hagen.de/MATHPHYS/veselic/downloads/j02.pdf

E a rotina LAPACK: http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-732F9EE1-BCEC-4D9B-9B93-AF5499B21140.htm#DRMAC08-1

Obviamente, o tamanho do trabalho é ajustável. Você pode chamar funções C do Python facilmente usando Cython, SWIG ou qualquer outro mecanismo de empacotamento.

Tolga Birdal
fonte