Quelle est la fonction LAPACK correspondante derrière Matlab [Q, R, E] = qr (A)?

Je tente actuellement de calculer à peu de frais une bonne estimation de rang pour une matrice . Par conséquent, je calcule une décompostion QR pivotante en utilisant$A$

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

dans Matlab. J'évalue le rang de utilisant$A$

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

Cela fonctionne bien et un tracé sur toutes les entrées diagonales de R ressemble à:

Si vous portez l'algorithme entier sur C / Fortran, je remplace [Q, R, E] = qr (A) en utilisant DGEQP3 de LAPACK, qui calcule également une décomposition QR pivotante de colonne. Mais si j'utilise la même estimation pour le rang, je me trompe surtout. Le même tracé pour le produit à partir de DGEQP3 ressemble $R$

La matrice d'entrée est exactement la même pour les deux expériences.

Ma question est maintenant sur quelle fonction LAPACK la décomposition QR pivotante de colonne de Matlab s'appuie-t-elle?

Merci pour toute aide, Grisu

Edit: DGEQPF donne le même mauvais résultat.

Edit2:

• La matrice d'entrée est dense et construite commeE + s i g n ( E , F $A$$E+sign(E,F)$
• $A$ est disponible ici: http://www-e.uni-magdeburg.de/makoehle/A.mtx.gz (Format MatrixMarket)
• Le mauvais : http://www-e.uni-magdeburg.de/makoehle/R_wrong.mtx.gz$R$
• J'ai utilisé LAPACK 3.4.0 avec OpenBlas / GotoBLAS (64 bits)
• Matlab 7, 2007b, 2010b Linux 32 bits

Edit3: - En utilisant GDB, j'ai découvert que Matlab 2010b appelle DGEQP3: # 3 0xaa46ce2f dans dgeqp3_ () à partir de /usr/ubuntu10.04/matlabr2010b/bin/glnx86/../../bin/glnx86/../. ./bin/glnx86/mllapack.so Pourquoi est-ce que j'obtiens le mauvais résultat en utilisant LAPACK (3.4.0 inclut les correctifs mentionnés dans la note de travail 176)?

Il y a deux problèmes à résoudre ici:

• Est dense ou clairsemé?$A$
• Avez-vous la même pile logicielle que les bibliothèques internes de MATLAB?

Dense ou clairsemé?

MATLAB ne mentionne plus explicitement les routines LAPACK qu'il appelle pour obtenir une factorisation QR si est dense. Si les informations de la documentation de MATLAB R2008b s'appliquent également aux versions ultérieures, MATLAB appelle depuis LAPACK lorsque vous appelez . Si est rare, MATLAB appelle SuiteSparseQR , hors du groupe de Tim Davis, qui est fourni avec UMFPACK dans la bibliothèque SuiteSparse.$A$DGEQP3[Q,R,E] = qr(A)$A$

Avez-vous la même pile logicielle que les bibliothèques internes de MATLAB?

Probablement pas, ce qui peut être l'une des raisons pour lesquelles vous obtenez des résultats différents.

J'ai rencontré ce problème lors du test d'unité d'une bibliothèque que j'écrivais et qui utilisait des factorisations QR. J'ai utilisé MATLAB pour prototyper mon travail et j'ai obtenu des résultats différents de ceux utilisant LAPACK ou NumPy. Pour autant que je sache, comme MathWorks ne facilite pas la recherche de ces informations, MATLAB utilise un verson de LAPACK au plus tôt en version 3.1.1 et la bibliothèque MKL BLAS d'Intel (pour Windows, Intel Mac et Linux) version 9.1 ou supérieur (voir ici ). Je n'ai rien trouvé sur la version de SuiteSparse MATLAB. En fouillant en ligne ou en consultant les fichiers de bibliothèque de votre système, vous pourrez glaner des informations supplémentaires. Vous pouvez essayer de changer les bibliothèques auxquelles MATLAB est lié afin de pouvoir comparer avec les mêmes bibliothèques à travers les progiciels; Eric Chu fournit une belle rédactionqui montre au moins comment vous pouvez remplacer la bibliothèque BLAS de MATLAB par la vôtre (bien sûr, vous le faites à vos risques et périls). Il suggère que vous pouvez également faire la même chose avec LAPACK. Il peut même être possible de remplacer la version de SuiteSparse que MATLAB utilise par votre propre version.

Les versions de LAPACK changent, ainsi que les versions de BLAS; ils peuvent utiliser des algorithmes différents d'une version à l'autre ou des conventions de classement différentes, même si avec orthogonal et est triangulaire supérieur, quelle que soit la version. Ces changements rendent la reproductibilité un défi. Même la tolérance que vous utilisez pour votre détermination de rang numérique est un appel au jugement; vous semblez utiliser une tolérance standard.Q $A = QR$$Q$$R$

J'ai fini par utiliser NumPy pour prototyper mes résultats pour la factorisation QR, car il utilise les bibliothèques système BLAS et LAPACK. NumPy et SciPy ne remplacent pas MATLAB, car les deux bibliothèques combinées manquent de certaines fonctionnalités de MATLAB, mais pour cette tâche d'algèbre linéaire particulière, Python + NumPy + SciPy + Matplotlib devrait bien fonctionner.

Voir la page de Leslie Foster sur les logiciels révélateurs de rang . Voir également cette note de travail LAPACK analysant les échecs des QR révélateurs de rangxGEQP3 .

Vous devriez être en mesure de savoir quelles routines MATLAB utilise en définissant des points d'arrêt dans un débogueur et en examinant la pile. La dernière fois que j'ai regardé, il y a plusieurs années, MATLAB a utilisé des bibliothèques partagées, auquel cas les noms des symboles ne peuvent pas être supprimés, vous verrez donc les noms des fonctions sur la pile des appels (mais pas les arguments car ils ne conservent définitivement pas les informations de débogage).