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 : Sélectionner tout - Visualiser dans une fenêtre à part
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.