GSL/BLAS: Multiplier une matrice avec une matrice inverse

Question

J'utilise la GNU GSL pour faire des calculs matriciels. J'essaie de multiplier une matrice B avec l'inverse d'une matrice A.

Maintenant, j'ai remarqué que la partie BLAS de GSL a une fonction pour ce faire, mais seulement si A est triangulaire. Y a-t-il une raison spécifique à cela? Aussi, quel serait le moyen le plus rapide de faire ce calcul? Dois-je inverser A en utilisant la décomposition LU, ou y a-t-il un meilleur moyen?

FWIW, A a la forme P 'G P, où P est une matrice normale, P' est son inverse, et G est une matrice diagonale.

Merci un tas :)

Answer 1

en un mot:

Le fait que seules matrices triangulaires sont pris en charge est simplement parce que cette opération est vraiment facile à réaliser pour une matrice triangulaire (elle est appelée back-substitution) et BLAS ne fournit que des routines pour les fonctions de bas niveau. Toute fonction de niveau supérieur se trouve généralement dans LAPACK qui utilise BLAS en interne pour les opérations par blocs.

Dans le cas spécifique où vous avez affaire: A:=P'*G*P, si P est une matrice normale, alors l'inverse de A peut s'écrire inv(A) := P'*inv(G)*P. Vous avez alors deux options:

1. Construire inv(A) et multiplier avec B.
2. séquentiellement multiplier B avec les différents facteurs de inv(A): multiplier B par P (à gauche), puis remettre à l'échelle chaque ligne du résultat avec l'inverse des éléments diagonaux de G, puis de nouveau se multiplier avec P' (nouveau à gauche).

Selon le cas, vous pouvez choisir votre méthode. Quoi qu'il en soit, vous aurez besoin de dgemm (multiplications matrice matricielle, Lvl3 BLAS) et dscal (mise à l'échelle des lignes, Lvl 1 BLAS), en supposant une double précision.

Espérons que cela aide.

A.

Answer 2

Je crois Adrien est juste avec la raison de ne pas avoir BLAS une fonction inverse pour les matrices carrées. Cela dépend de la matrice que vous utilisez pour optimiser le calcul de son inverse.

En général, vous pouvez utiliser la décomposition LU, qui est valable pour n'importe quelle matrice carrée. . I. e, quelque chose comme:

gsl_linalg_LU_decomp(A, p, signum); 
gsl_linalg_LU_invert(A, p, invA); 

où A est la matrice carrée que vous voulez son inverse, p est un objet gsl_permutation (l'objet de permutation où la matrice de permutation est codée), Signum est le signe de la permutation et invA est l'inverse de A.

Puisque vous déclarez que A = P' G P étant P normal et G diagonale, probablement A est une matrice normale. Je ne les ai pas utilisés depuis un moment, mais il doit y avoir un théorème de factorisation pour eux, que vous pouvez trouver dans le Chapter 14 of the GSL reference manual ou mieux, dans un livre d'algèbre linéaire.Par exemple, si vous aviez des matrices positives positives symétriques et que vous vouliez trouver l'inverse, vous pourriez utiliser la factorisation de Cholesky, qui est optimisée pour ce type de matrices. Ensuite, vous pouvez utiliser les fonctions gsl_linalg_cholesky_decomp() et gsl_linalg_cholesky_invert() dans le GSL pour le rendre efficace.

J'espère que ça aide!