Un certain nombre d’op ́erations classiques sur les matrices sont pr ́esentes dans la version de base de R. Avant de les pr ́esenter, d ́efinissons un scalaire λ et deux matrices r ́eelles A et B, et une matrice complexe C dont nous nous servirons par la suite. Le lecteur int ́eress ́e par ce sujet pourra consulter avec profi
Nous allons créer un \(\lambda\) et des matrices plus précisement la matrice réelle, rééllesumétrique, complexe hermitienne et la matrice identité.
lambda <- 2 # Création du scalaire λ.
A <- matrix(c(2,3,5,4),nrow=2,ncol=2) # Matrice réelle.
A ## [,1] [,2]
## [1,] 2 5
## [2,] 3 4B <- matrix(c(1,2,2,7),nrow=2,ncol=2) # Matrice réellesymétrique.
B## [,1] [,2]
## [1,] 1 2
## [2,] 2 7C <- matrix(c(1,1i,-1i,3),ncol=2) # Matrice complexe hermitienne.
C## [,1] [,2]
## [1,] 1+0i 0-1i
## [2,] 0+1i 3+0iI2 <- diag(rep(1,2)) # Matrice identité d’ordre 2.
I2## [,1] [,2]
## [1,] 1 0
## [2,] 0 1Nous allons nous servir de ces objets pour illustrer les diff ́erentes op ́erations possibles sur les matrices.
Pour des fonctions plus élaboréé des de manipulation de matrices, le lecteur pourra utiliser le package Matrix.
Les op ́erations de base sur les matrices en R sont les suivantes : -L’addition d’un scalaire : λ + A
lambda+A## [,1] [,2]
## [1,] 4 7
## [2,] 5 6-L’addition (terme `a terme) : A + B
A+B## [,1] [,2]
## [1,] 3 7
## [2,] 5 11-La soustraction (terme `a terme) : A − B
A-B## [,1] [,2]
## [1,] 1 3
## [2,] 1 -3-La multiplication par un scalaire : λA
lambda*A## [,1] [,2]
## [1,] 4 10
## [2,] 6 8-La transposition : \(A^t\)
t(A)## [,1] [,2]
## [1,] 2 3
## [2,] 5 4La conjugu ́ee :C
Conj(C)## [,1] [,2]
## [1,] 1+0i 0+1i
## [2,] 0-1i 3+0i-La multiplication terme `a terme :
A*B## [,1] [,2]
## [1,] 2 10
## [2,] 6 28-La multiplication matricielle
A%*%B## [,1] [,2]
## [1,] 12 39
## [2,] 11 34-La division terme `a terme :
A/B## [,1] [,2]
## [1,] 2.0 2.5000000
## [2,] 1.5 0.5714286-L’inverse matricielle :\(B^{-1}\)
solve(B)## [,1] [,2]
## [1,] 2.3333333 -0.6666667
## [2,] -0.6666667 0.3333333-La division matricielle \(A^{-1}B\)
solve(A)%*%B## [,1] [,2]
## [1,] 0.8571429 3.857143
## [2,] -0.1428571 -1.142857-Le produit avec la transposé: \(A^tB\)
crossprod(A,B) # ou t(A)%*%B## [,1] [,2]
## [1,] 8 25
## [2,] 13 38Exemple :
M <- matrix(c(1,2,4,1,3,1),nrow = 3,ncol = 2)
N <- matrix(c(3,1,4,4,3,1),nrow = 3,ncol = 2)
O <- matrix(c(3,1,4,3,2,2),nrow = 2,ncol = 3)
P <- matrix(c(3,1,1,4,3,2,2,2,1),nrow = 3,ncol = 3)Nous allons donner les dimesnsions des matrix M,N,O et P
dim(M)## [1] 3 2dim(N)## [1] 3 2dim(O)## [1] 2 3dim(P)## [1] 3 3Npous calculons M+N,M-N, 3M,MO,OM,\(M^t\) et \(p^-{1}\)
M+N## [,1] [,2]
## [1,] 4 5
## [2,] 3 6
## [3,] 8 2M-N## [,1] [,2]
## [1,] -2 -3
## [2,] 1 0
## [3,] 0 03*M## [,1] [,2]
## [1,] 3 3
## [2,] 6 9
## [3,] 12 3M%*%O## [,1] [,2] [,3]
## [1,] 4 7 4
## [2,] 9 17 10
## [3,] 13 19 10O%*%M## [,1] [,2]
## [1,] 19 17
## [2,] 15 12t(M)## [,1] [,2] [,3]
## [1,] 1 2 4
## [2,] 1 3 1solve(P)## [,1] [,2] [,3]
## [1,] 1 1.480297e-16 -2
## [2,] -1 -1.000000e+00 4
## [3,] 1 2.000000e+00 -5Nous vérifions que \(PP^{-1}=I_{3}=P^{-1}P\)
P%*%solve(P)## [,1] [,2] [,3]
## [1,] 1.000000e+00 0 0
## [2,] 0.000000e+00 1 0
## [3,] 1.110223e-16 0 1solve(P)%*%P## [,1] [,2] [,3]
## [1,] 1.000000e+00 -8.881784e-16 -4.440892e-16
## [2,] 4.440892e-16 1.000000e+00 0.000000e+00
## [3,] 0.000000e+00 0.000000e+00 1.000000e+00x <- seq(1,4)
y <- seq(4,7)
outer(x,y,FUN="*")## [,1] [,2] [,3] [,4]
## [1,] 4 5 6 7
## [2,] 8 10 12 14
## [3,] 12 15 18 21
## [4,] 16 20 24 28kronecker(A,B)## [,1] [,2] [,3] [,4]
## [1,] 2 4 5 10
## [2,] 4 14 10 35
## [3,] 3 6 4 8
## [4,] 6 21 8 28Il est parfois utile de r ́ecup ́erer les sous-matrices triangulaires inf ́erieure et sup ́erieure d’une matrice. Cela est possible au moyen des fonctions lower.tri() et upper.tri().
M <- matrix(1:16,nrow=4)
lower.tri(M)## [,1] [,2] [,3] [,4]
## [1,] FALSE FALSE FALSE FALSE
## [2,] TRUE FALSE FALSE FALSE
## [3,] TRUE TRUE FALSE FALSE
## [4,] TRUE TRUE TRUE FALSEupper.tri(M,diag=TRUE)## [,1] [,2] [,3] [,4]
## [1,] TRUE TRUE TRUE TRUE
## [2,] FALSE TRUE TRUE TRUE
## [3,] FALSE FALSE TRUE TRUE
## [4,] FALSE FALSE FALSE TRUEM[lower.tri(M)] <- 0
M## [,1] [,2] [,3] [,4]
## [1,] 1 5 9 13
## [2,] 0 6 10 14
## [3,] 0 0 11 15
## [4,] 0 0 0 16L’op ́erateur matriciel vec appliaqu ́e a une matrice A consistea fabriquer lelong vecteur colonne vec(A) constitu ́e de l’empilement des colonnes de A. Il s’obtient en R au moyen de l’instruction suivante :
vec <- function(M) as.matrix(as.vector(M))
# ou de façon équivalente, mais à l’extérieur d’une fonction:
# dim(A) <- c(prod(dim(A)),1)
A## [,1] [,2]
## [1,] 2 5
## [2,] 3 4vec(A)## [,1]
## [1,] 2
## [2,] 3
## [3,] 5
## [4,] 4L’op ́erateur matriciel vech (pour vec half ou demi-vec) appliqu ́e a une ma-trice A consistea fabriquer le long vecteur colonne vech(A) constitu ́e de l’em-pilement des colonnes de A, mais en excluant les ́el ́ements au-dessus de la diagonale de A. Il s’obtient en R au moyen de l’instruction suivante :
vech <- function(M) as.matrix(M[lower.tri(M,diag=TRUE)])
vech(A)## [,1]
## [1,] 2
## [2,] 3
## [3,] 4La fonction det() permet de calculer le d ́eterminant d’une matrice.
det(A)## [1] -7Il n’existe pas de fonction R permettant de calculer directement la trace d’une matrice, mais il est tr`es facile de la calculer ainsi :
sum(diag(A))## [1] 6NB:Il ne faut pas utiliser la fonction trace() pour calculer la trace d’une matrice. En effet, cette fonction sert dans le d ́ebogage de code R.
Le nombre de conditionnement est le rapport de la plus grande sur la plus petite valeur singuli`ere non nulle. Une grande valeur du nombre de condition-nement indique de mauvaises propri ́et ́es num ́eriques de la matrice. Il s’obtientau moyen de la fonction kappa().
kappa(A,exact=TRUE)## [1] 7.582401La fonction scale() permet de centrer et/ou de r ́eduire une matrice. Le centrage consiste a retranchera chacune des colonnes de la matrice la moyenne de cette colonne. La r ́eduction consiste `a diviser chacune des colonnes de la matrice par son ́ecart type
NB:Notez que la fonction sd() calcule un ́ecart type avec un num ́erateur ́egal `a n − 1.
# Centre
scale(A,scale=FALSE)## [,1] [,2]
## [1,] -0.5 0.5
## [2,] 0.5 -0.5
## attr(,"scaled:center")
## [1] 2.5 4.5#Réductee
scale(A,center=FALSE,scale=apply(A,2,sd))## [,1] [,2]
## [1,] 2.828427 7.071068
## [2,] 4.242641 5.656854
## attr(,"scaled:scale")
## [1] 0.7071068 0.7071068NB: Si l’on veut obtenir une r ́eduction fond ́ee sur l’ ́ecart type de la popula- tion comme cela est pr ́econis ́e par exemple dans l’analyse de donn ́ees `a la fran ̧caise, on pourra utiliser l’instruction :
red <- sqrt((nrow(A)-1)/nrow(A))
scale(A,center=FALSE,scale=apply(A,2,sd)*red)## [,1] [,2]
## [1,] 4 10
## [2,] 6 8
## attr(,"scaled:scale")
## [1] 0.5 0.5# ou de façon équivalente: t(A/apply(A,2,sd))/redOn peut obtenir les valeurs propres et les vecteurs propres d’une matrice au moyen de la fonction eigen().
eigen(A)## eigen() decomposition
## $values
## [1] 7 -1
##
## $vectors
## [,1] [,2]
## [1,] -0.7071068 -0.8574929
## [2,] -0.7071068 0.5144958Notez que pour une matrice hermitienne C (c’est-a-dire une matrice complexe ́egalea sa transconjugu ́ee), la fonction eigen() permet d’obtenir la d ́ecomposition aux valeurs propres de cette matrice C, c’est-a-dire C =VDV∗(ou V∗ est la matrice adjointe de V) :
C <- matrix(c(1,1i,-1i,3),ncol=2)
e <- eigen(C,symmetric=TRUE)
V <- e$vectors
D <- diag(e$values)
all.equal(C,V%*%D%*%t(Conj(V)))## [1] TRUEUne racine carr ́ee d’une matrice d ́efinie positive C est toute matrice M quiv ́erifie M∗M = C, ou M∗d ́esigne la matrice adjointe (transpos ́ee conjugu ́ee)de M. On la note en g ́en ́eral C1/2 mˆeme si elle n’est pas unique. Lorsque C hermitienne (c’est-a-dire ou bien une matrice complexe ́egale `a sa transconju-gu ́ee, ou bien une matrice r ́eelle sym ́etrique), on peut calculer C1/2 de la fa ̧consuivante :est
e <- eigen(C,symmetric=TRUE)
V <- e$vectors
V %*% diag(sqrt(e$values)) %*% t(Conj(V)) # C1/2,## [,1] [,2]
## [1,] 0.9238795+0.0000000i 0.000000-0.3826834i
## [2,] 0.0000000+0.3826834i 1.689246+0.0000000i# qui est ici
# hermitienne.La matrice \(C^{−1/2}\) peut se calculer ainsi :
V %*% diag(1/sqrt(e$values)) %*% t(Conj(V)) # C−1/2,## [,1] [,2]
## [1,] 1.194478+0.0000000i 0.0000000+0.2705981i
## [2,] 0.000000-0.2705981i 0.6532815+0.0000000i# qui est ici
# hermitienne.On cherche a ́ecrire C = UDV∗ou D est la matrice diagonale des valeurs singulieres de C, U (respectivement V) est la matrice des vecteurs singuliers a gauche (respectivementa droite) de C. Pour cela, il faut utiliser la fonction svd().
res <- svd(C)
res## $d
## [1] 3.4142136 0.5857864
##
## $u
## [,1] [,2]
## [1,] -0.3826834+0.0000000i 0.9238795+0.0000000i
## [2,] 0.0000000-0.9238795i 0.0000000-0.3826834i
##
## $v
## [,1] [,2]
## [1,] -0.3826834+0.0000000i 0.9238795+0.0000000i
## [2,] 0.0000000-0.9238795i 0.0000000-0.3826834iD <- diag(res$d)
U <- res$u
V <- res$v
all.equal(C,U%*%D%*%t(Conj(V))) # A = UDV∗## [1] TRUEPour calculer l’inverse (g ́en ́eralis ́ee) de Moore-Penrose d’une matrice (non inversible), on peut utiliser la fonction suivante :
mpinv <- function(M,eps=1e-13) {
s <- svd(M)
e <- s$d
e[e>eps] <- 1/e[e>eps]
return(s$v%*%diag(e)%*%t(s$u))
}Pour une matrice r ́eelle sym ́etrique d ́efinie positive B, on cherche a ́ecrire B = UTU = LLT, ou U (respectivement L) est une matrice triangulaire sup ́erieure (respectivement inf ́erieure). Vous aurez not ́e que cela implique que U est une racine carr ́ee de B. Pour obtenir cette d ́ecomposition, il faut utiliser la fonction chol().
U <- chol(B) # C’est un autre moyen d’obtenir# B1/2.
L <- t(U)
U## [,1] [,2]
## [1,] 1 2.000000
## [2,] 0 1.732051 all.equal(B,t(U)%*%U) # B = UTU.## [1] TRUENotez ́egalement que vous pouvez utiliser la fonction chol2inv() pour cal-culer l’inverse \(B^{-1}\)d’une matrice carr ́ee sym ́etrique d ́efinie positive B, `a partirde sa d ́ecomposition de Cholesky.
B## [,1] [,2]
## [1,] 1 2
## [2,] 2 7chol2inv(U) # C’est $B^{-1}$.## [,1] [,2]
## [1,] 2.3333333 -0.6666667
## [2,] -0.6666667 0.3333333all.equal(chol2inv(U),solve(B))## [1] TRUENotez enfin que la fonction chol() peut ˆetre utilis ́ee pour calculer$ B^{−1/}2$dela fa ̧con suivante :
solve(chol(B))# C’est une version de B−1/2.## [,1] [,2]
## [1,] 1 -1.1547005
## [2,] 0 0.5773503On cherche a ́ecrire A = QR, ou Q est une matrice orthogonale (QQT =QTQ = I) et R est une matrice triangulaire sup ́erieure.
res <- qr(A)
Q <- qr.Q(res)
Q## [,1] [,2]
## [1,] -0.5547002 -0.8320503
## [2,] -0.8320503 0.5547002all.equal(I2,Q%*%t(Q)) # I2 = QQT## [1] TRUER <- qr.R(res)
R## [,1] [,2]
## [1,] -3.605551 -6.101702
## [2,] 0.000000 -1.941451all.equal(A,Q%*%R) # A = QR.## [1] TRUENotez que qr(A,tol=1e-07)$rank renvoie le rang de la matrice stock ́ee dans A, en utilisant l’argument de tol ́erance tol pour d ́etecter des d ́epen- dances lin ́eaires dans les colonnes de A.