LU décomposition en R et Python

Question

J'ai la matrice suivante A dans R:

  # [,1]  [,2] [,3]  [,4] 
# [1,] -1.1527778 0.4444444 0.375 0.3333333 
# [2,] 0.5555556 -1.4888889 0.600 0.3333333 
# [3,] 0.6250000 0.4000000 -1.825 0.8000000 
# [4,] 0.6666667 0.6666667 0.200 -1.5333333 

A <- structure(c(-1.15277777777778, 0.555555555555556, 0.625, 0.666666666666667, 
0.444444444444444, -1.48888888888889, 0.4, 0.666666666666667, 
0.375, 0.6, -1.825, 0.2, 0.333333333333333, 0.333333333333333, 
0.8, -1.53333333333333), .Dim = c(4L, 4L), .Dimnames = list(NULL, 
    NULL)) 

Je Calculons son LU-décomposition comme suit:

library(Matrix) 
ex <- expand(lu(t(A))) 
L <- ex$L 
P <- ex$P 
C <- U <- ex$U 
C[lower.tri(U)] <- L[lower.tri(L)] 

print(C) 

# 4 x 4 Matrix of class "dgeMatrix" 
      # [,1]  [,2]  [,3]   [,4] 
# [1,] -1.1527778 0.5555556 0.6250000 6.666667e-01 
# [2,] -0.3855422 -1.2746988 0.6409639 9.236948e-01 
# [3,] -0.3253012 -0.6124764 -1.2291115 9.826087e-01 
# [4,] -0.2891566 -0.3875236 -1.0000000 -2.220446e-16 

D'autre part, c'est le Python code pour le même travail:

lu, piv = scipy.linalg.lu_factor(A.T, check_finite=False) 

print(lu) 

# [[ -1.15277778e+00 5.55555556e-01 6.25000000e-01 6.66666667e-01] 
# [ -3.85542169e-01 -1.27469880e+00 6.40963855e-01 9.23694779e-01] 
# [ -2.89156627e-01 -3.87523629e-01 1.22911153e+00 -9.82608696e-01] 
# [ -3.25301205e-01 -6.12476371e-01 -1.00000000e+00 7.69432827e-16]] 

Je me demande pourquoi les deux C et lu matrices dans R et Python (respectivement) ne sont pas les mêmes. Le point est que je dois obtenir le même résultat que la version Python (c'est-à-dire la matrice lu). Avez-vous une idée de ce que je fais mal?

Answer 1

Vous n'avez pas de matrice complète.

det(A) 
# [1] 3.145929e-15 

La factorisation LU n'est pas prévue pour une matrice carrée à défilement de rang. Généralement à travers le logiciel, la façon dont une matrice déficiente est traitée est différente, alors ne soyez pas si surpris de la différence. Pourquoi ne pas essayer une matrice complète de rang, comme ce qui suit:

A <- structure(c(0.923065107548609, 0.922819485189393, 0.277002309216186, 
0.532856695353985, 0.481061384081841, 0.0952619954477996, 
0.261916425777599, 0.433514681644738, 0.677919807843864, 
0.771985625848174, 0.705952850636095, 0.873727774480358, 
0.28782021952793, 0.863347264472395, 0.627262107795104, 
0.187472499441355), .Dim = c(4L, 4L)) 

Si vous voyez toujours la différence entre R et Python, lisez la suite.

En fait, R et Python donnent la factorisation correcte. Cependant, la factorisation LU n'est pas unique (même en cas de pivotement de ligne), car les pivots peuvent différer en signe. Par exemple, le 3ème pivot est -1.23 dans R mais 1.23 dans Python. Une fois le pivot différent, la matrice d'élimination, c'est-à-dire la matrice L, sera très différente. Vous n'avez pas besoin d'être si pointilleux que vous voulez quelque chose de pareil. Des logiciels différents peuvent utiliser des algorithmes différents.