Fonction récursive pour implémenter une technique d'inversion de blocs dans R

Question

Je suis très novice en algèbre linéaire et j'essaye d'implémenter une fonction récursive qui inverse n'importe quelle matrice en utilisant la technique d'inversion de blocs sans utiliser la bibliothèque R " résoudre".

Cette question était déjà réponse dans le post suivant: function for matrix

Cependant, cela n'a pas fonctionné pour moi et j'ai essayé de mettre en œuvre ma propre version:

matrixInversion <- function(M){ 
    if(nrow(M) == 2){ 
     a <- M[1,1] 
     b <- M[1,2] 
     c <- M[2,1] 
     d <- M[2,2] 
     deter <- ((a*d)-(b*c)) 
     InverseMatrix <- ((1/deter)*matrix(c(d,-c,-b,a),nrow=2,ncol=2)) 
    } else { 
     x <- (floor(nrow(M)/2)) 
     A <- M[1:x, 1:x, drop=F] 
     B <- M[1:x, -1:-x, drop=F] 
     C <- M[-1:-x, 1:x, drop=F] 
     D <- M[-1:-x, -1:-x, drop=F] 

     Ainv <- matrixInversion(A) 
     common <- matrixInversion(D - C %*% Ainv %*% B) 
     newA <- Ainv+Ainv%*%B%*%common%*%C%*%Ainv 
     newB <- (-Ainv)%*%B%*%common 
     newC <- (-common)%*%C%*%Ainv 
     newD <- (-common) 

     result <- cbind(rbind(newA, newC), rbind(newB, newD)) 
    } 
} 

Cette version fonctionne avec succès pour la matrice avec nombre pair de colonnes, mais pas avec la matrice avec le nombre impair de colonnes. Je ne peux pas comprendre comment l'implémenter correctement. Toute suggestion?

Merci!

Answer 1

En plus des points faits dans les commentaires, vous invalidez incorrectement newD. Votre fonction des changements fonctionne maintenant avec les deux matrices impairs rangs et même rangée:

matrixInversion <- function(M){ 
    if (nrow(M) == 1){ 
    return(1/M) 
    } else if (nrow(M) == 2){ 
    a <- M[1,1] 
    b <- M[1,2] 
    c <- M[2,1] 
    d <- M[2,2] 
    deter <- ((a*d)-(b*c)) 
    return((1/deter)*matrix(c(d,-c,-b,a),nrow=2,ncol=2)) 
    } else { 
    x <- (floor(nrow(M)/2)) 
    A <- M[1:x, 1:x, drop=F] 
    B <- M[1:x, -1:-x, drop=F] 
    C <- M[-1:-x, 1:x, drop=F] 
    D <- M[-1:-x, -1:-x, drop=F] 
    Ainv <- matrixInversion(A) 
    common <- matrixInversion(D - C %*% Ainv %*% B) 
    newA <- Ainv+Ainv%*%B%*%common%*%C%*%Ainv 
    newB <- (-Ainv)%*%B%*%common 
    newC <- (-common)%*%C%*%Ainv 
    newD <- common 
    return(cbind(rbind(newA, newC), rbind(newB, newD))) 
    } 
} 

## random seed not relevant here ... 
m1 <- matrix(runif(150^2), nr=150) 
all.equal(solve(m1), matrixInversion(m1)) 
# [1] TRUE 
any(is.nan(matrixInversion(m1))) # based on your comment 
# [1] FALSE 

m2 <- matrix(runif(151^2), nr=151) 
all.equal(solve(m2), matrixInversion(m2)) 
# [1] TRUE 
any(is.nan(matrixInversion(m2))) 
# [1] FALSE 

Pour mémoire, je trouve cela en vérifiant progressivement des matrices plus grandes. Quelque chose comme:

m1 <- matrix(runif(150^2), nr=150) 
i <- 4; all.equal(solve(m1[1:i,1:i]), matrixInversion(m1[1:i,1:i])) 

En 2, il était égal, et il a commencé à casser au-dessus de cela. J'ai utilisé debugonce(matrixInversion) pour vérifier à l'intérieur des deux premiers blocs conditionnels, mais ils allaient bien. J'ai ensuite comparé vos équations avec le wiki-page initialement cité, et j'ai remarqué votre calcul incorrect. Bien sûr, c'était le dernier que j'ai vérifié :-)