Matrice d'inversion Méthodes

Question

Quand on a un problème d'une multiplication de matrice inverse avec un vecteur, comme tel:

peut prendre une décomposition Cholesky de A et backsubstitute b pour trouver le vecteur résultant x. Cependant, une matrice inverse est parfois nécessaire lorsque le problème n'est pas formulé comme ci-dessus. Ma question est de savoir quelle est la meilleure façon de gérer une telle situation. Ci-dessous, j'ai comparé différentes façons (en utilisant numpy) pour inverser une matrice définie positive:

Tout d'abord, générer la matrice:

>>> A = np.random.rand(5,5) 
>>> A 
array([[ 0.13516074, 0.2532381 , 0.61169708, 0.99678563, 0.32895589], 
     [ 0.35303998, 0.8549499 , 0.39071336, 0.32792806, 0.74723177], 
     [ 0.4016188 , 0.93897663, 0.92574706, 0.93468798, 0.90682809], 
     [ 0.03181169, 0.35059435, 0.10857948, 0.36422977, 0.54525 ], 
     [ 0.64871162, 0.37809219, 0.35742865, 0.7154568 , 0.56028468]]) 
>>> A = np.dot(A.transpose(), A) 
>>> A 
array([[ 0.72604206, 0.96959581, 0.82773451, 1.10159817, 1.05327233], 
     [ 0.96959581, 1.94261607, 1.53140854, 1.80864185, 1.9766411 ], 
     [ 0.82773451, 1.53140854, 1.52338262, 1.89841402, 1.59213299], 
     [ 1.10159817, 1.80864185, 1.89841402, 2.61930178, 2.01999385], 
     [ 1.05327233, 1.9766411 , 1.59213299, 2.01999385, 2.10012097]]) 

Les résultats de la méthode d'inversion directe sont les suivantes:

>>> np.linalg.inv(A) 
array([[ 5.49746838, -1.92540877, 2.24730018, -2.20242449, 
     -0.53025806], 
     [ -1.92540877, 95.34219156, -67.93144606, 50.16450952, 
     -85.52146331], 
     [ 2.24730018, -67.93144606, 57.0739859 , -40.56297863, 
     58.55694127], 
     [ -2.20242449, 50.16450952, -40.56297863, 30.6441555 , 
     -44.83400183], 
     [ -0.53025806, -85.52146331, 58.55694127, -44.83400183, 
     79.96573405]]) 

Lorsque vous utilisez le pseudo-inverse, les résultats sont les suivants (vous remarquerez peut-être que la précision affichée, les résultats sont les mêmes que l'inversion directe):

>>> np.linalg.pinv(A) 
array([[ 5.49746838, -1.92540877, 2.24730018, -2.20242449, 
     -0.53025806], 
     [ -1.92540877, 95.34219156, -67.93144606, 50.16450952, 
     -85.52146331], 
     [ 2.24730018, -67.93144606, 57.0739859 , -40.56297863, 
     58.55694127], 
     [ -2.20242449, 50.16450952, -40.56297863, 30.6441555 , 
     -44.83400183], 
     [ -0.53025806, -85.52146331, 58.55694127, -44.83400183, 
     79.96573405]]) 

Enfin, lors de la résolution de la matrice d'identité:

>>> np.linalg.solve(A, np.eye(5)) 
array([[ 5.49746838, -1.92540877, 2.24730018, -2.20242449, 
     -0.53025806], 
     [ -1.92540877, 95.34219156, -67.93144606, 50.16450952, 
     -85.52146331], 
     [ 2.24730018, -67.93144606, 57.0739859 , -40.56297863, 
     58.55694127], 
     [ -2.20242449, 50.16450952, -40.56297863, 30.6441555 , 
     -44.83400183], 
     [ -0.53025806, -85.52146331, 58.55694127, -44.83400183, 
     79.96573405]]) 

Encore une fois, vous remarquerez peut-être que sur une inspection rapide, le résultat est le même que les deux méthodes précédentes.

Il est bien connu que l'inversion matricielle est un problème mal posé en raison de l'instabilité numérique qui devrait être évitée autant que possible. Cependant, dans les situations où cela semble inévitable, quelle est l'approche préférable et pourquoi? Pour clarifier, je me réfère à la meilleure approche lors de la mise en œuvre de telles équations dans les logiciels. Un exemple d'un tel problème est fourni avec un autre my questions

Answer 1

La raison pour éviter d'inverser les matrices n'a qu'un rapport avec l'efficacité. Il est plus rapide de résoudre directement les systèmes linéaires. Si vous pensez un peu différemment au problème de votre question liée, vous pouvez appliquer les mêmes principes.

Afin de trouver la matrice inv(K) * Y * T(Y) * inv(K) - D * inv(K) vous pouvez résoudre les systèmes d'équations suivantes:

K * R * K = Y * T(Y) 

Vous pouvez le résoudre en deux parties:

R2 * K = R1 
K * R1 = Y * T(Y) 

Alors vous devez d'abord résoudre pour R1 avec votre méthode habituelle, puis résolvez pour R2 (reconnaissez que vous pouvez résoudre T(K) * T(R2) = T(R1) si vous devez).

Cependant, à ce stade, je ne sais pas si cela sera plus efficace que de calculer explicitement l'inverse sauf si K est symétrique.(Il peut y avoir un moyen d'obtenir efficacement la décomposition de T(K) de K, mais je ne sais pas désinvolture)

Si K est symétrique, vous pouvez calculer votre décomposition sur K une fois et le réutiliser pour les deux arrière-substitution étapes et il pourrait être plus efficace que de calculer l'inverse explicitement.