Alternative(s) quand impossible d'inverser une matrice avec "solve"?

Version imprimable

Bonjour à tous,

Je programme actuellement une fonction itérative dans laquelle j'ai besoin de calculer l'inverse d'une matrice. Cependant, il arrive dans certains cas que je ne puisse avoir recours à la fonction "solve". Pour un tel cas j'ai vu plusieurs alternatives : soit utiliser la fonction "solvecov" du package "fpc", soit utiliser une méthode basée sur la décomposition en valeurs singulières (svd) que j'ai vue dans l'ouvrage "Comprendre et réaliser les tests statistiques à l'aide de R", de Gaël Millot.
Le problème est que ces deux méthodes me donnent des résultats différents et je ne sais pas trop comment les départager, donc si quelqu'un connaît des arguments en faveur ou contre une de ces méthodes, je suis preneuse.

Comme un exemple vaut mieux qu'un long discours :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
# création d'un jeu de données test
obs1<-matrix(rep(0,300),nrow=50)
obs1[,1]<-rep(1,50)
for(i in 2:6)
{
	obs1[(10*(i-2)+(1:10)),i]<-1
}
X<-t(obs1)%*%obs1
 
X
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]   50   10   10   10   10   10
[2,]   10   10    0    0    0    0
[3,]   10    0   10    0    0    0
[4,]   10    0    0   10    0    0
[5,]   10    0    0    0   10    0
[6,]   10    0    0    0    0   10
 
solve(X)
Error in solve.default(X) : 
  Lapack routine dgesv: system is exactly singular
# on s'en doutait un peu ^^
 
 
## Alternative 1 : solvecov ##
 
library(fpc)
?solvecov
 
# Description:
#
#     Tries to invert a matrix by ‘solve’. If this fails because of
#     singularity, an eigenvector decomposition is computed, and
#     eigenvalues below ‘1/cmax’ are replaced by ‘1/cmax’, i.e., ‘cmax’
#     will be the corresponding eigenvalue of the inverted matrix.
 
Y<-solvecov(X)$inv
Y
            [,1]        [,2]        [,3]        [,4]        [,5]        [,6]
[1,]  1666666667 -1666666667 -1666666667 -1666666667 -1666666667 -1666666667
[2,] -1666666667  1666666667  1666666667  1666666667  1666666667  1666666667
[3,] -1666666667  1666666667  1666666667  1666666667  1666666667  1666666667
[4,] -1666666667  1666666667  1666666667  1666666667  1666666667  1666666667
[5,] -1666666667  1666666667  1666666667  1666666667  1666666667  1666666667
[6,] -1666666667  1666666667  1666666667  1666666667  1666666667  1666666667
 
 
## Alternative 2 : SVD ##
 
S<-svd(X)
val.crit<-which(S$d<10^-8)
Diag.mod<-diag(1/S$d)
for(i in val.crit)
{
	Diag.mod[i,i]<-0
}
Z<-S$v%*%Diag.mod%*%t(S$u)
 
Z
            [,1]         [,2]         [,3]         [,4]         [,5]         [,6]
[1,] 0.013888889  0.002777778  0.002777778  0.002777778  0.002777778  0.002777778
[2,] 0.002777778  0.080555556 -0.019444444 -0.019444444 -0.019444444 -0.019444444
[3,] 0.002777778 -0.019444444  0.080555556 -0.019444444 -0.019444444 -0.019444444
[4,] 0.002777778 -0.019444444 -0.019444444  0.080555556 -0.019444444 -0.019444444
[5,] 0.002777778 -0.019444444 -0.019444444 -0.019444444  0.080555556 -0.019444444
[6,] 0.002777778 -0.019444444 -0.019444444 -0.019444444 -0.019444444  0.080555556

Merci d'avance ! :)

Cordialement,

A.D.

22/11/2011, 09h44
lilly74

Bonjour A.D.,

je ne connais pas ces méthodes.

Mais selon ce que vous voulez obtenir, avez-vous pensé à:

1/ la décomposition de Choleski? (si matrice symétrique définie positive) -package matrix, il me semble-
2/ la décomposition QR pour déterminer la pseudo-inverse?

Bon courage et bonne journée :)

Lilly.
22/11/2011, 11h17
A. D.

Bonjour Lilly74,

Tout d'abord merci pour votre réponse :)

En fait le truc c'est que je ne sais pas trop à quoi va ressembler ma matrice (cela va dépendre des paramètres d'entrée de ma fonction et comme il s'agit d'une procédure itérative, il y aura des modifications au cours du temps) et je ne sais donc pas trop quelle approche choisir.

A priori ma matrice va être symétrique semi-définie positive (car matrice de variance-covariance), mais je ne suis pas sûre que cette propriété soit impérativement vérifiée à chaque itération.

Savez-vous s'il y a des "conditions" qui feraient qu'il est préférable de choisir une approche plutôt qu'une autre ou bien cela dépend-t-il essentiellement de ce qu'on souhaite obtenir?

Cordialement,

A.D.
22/11/2011, 11h40
lilly74

Bonjour,

je suis désolée, cela fait trop longtemps que je n'ai pas fait d'algèbre et mes pistes reposaient sur de vieux souvenirs :oops:

Je pense que d'autres que moi seront mieux placés pour donner une réponse précise.

Pour ma part, si l'on ne sait pas quelle est la nature de la matrice (définie positive ou pas; cf par exemple la discussion http://forums.futura-sciences.com/ma...reguliere.html), je pencherais pour une décomposition QR qui n'impose pas des hypothèses sur cela...
Mais la pseudo-inverse (http://en.wikipedia.org/wiki/Pseudoi...ldid=112040404) vous conviendrait-elle?

Désolée de ne pouvoir en dire plus sans mes vieux bouquins :roll:
Bonne continuation :)
22/11/2011, 12h35
A. D.

Re-bonjour Lilly74,

Pas de soucis, c'est déjà très gentil d'essayer de me répondre :)

Moi c'est pareil, je pense qu'il va falloir que je me replonge aussi dans mes cours à ce sujet ! Cela dit, si d'autres personnes ont des avis à donner sur le sujet (notamment sur l'utilisation de la fonction "solvecov" du package "fpc"), vous êtes les bienvenus :mrgreen:

Cordialement,

A.D.