Approximation de matrices

Dans le cours d’analyses des données, on va commencer l’analyse en composantes principales, en lien avec l’analyse spectrale de  ou de , ou directement l’analyse en valeurs singulières de .

Reprenons les données de l ‘enquête du quotidien L’Express,

> BD=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/securite.txt",
header=TRUE)
>  M=as.matrix(BD[,2:10])
>  colnames(M)=names(BD)[2:10]
>  rownames(M)=BD$dep

Rappelons que l’on a une matrice de taille

> dim(M)
[1] 94  9

La décomposition en valeurs singulières donne des matrices

> SVD=svd(M)
> dim(M)
[1] 94  9
> dim(SVD$u)
[1] 94  9
> dim(SVD$v)
[1] 9 9

telles que . La preuve, créons la matrice produit

> M1 = SVD$u %*% diag(SVD$d) %*% t(SVD$v)

Pour vérifier que les deux matrices sont (très très) proches, on a le choix de la distance. Par exemple

à partir de la norme de Frobenius, .

> sum( (M0-M)^2 )
[1] 2.160588e-25

On peut aussi utiliser

qui est la plus grande valeur singulière de .

> max(eigen(t(M0-M)%*%(M0-M))$value)
[1] 2.018236e-25

Pour (simplement) approcher la matrice de départ par une matrice de rang plus faible, on peut  utiliser une approximation de la forme

où seules les  premières colonnes sont gardées,

> k=4
> M4 = SVD$u[,1:k] %*% diag(SVD$d[1:k]) %*% t(SVD$v[,1:k])

Les deux matrices semblent proches. Par exemple, si on regarde la première ligne on a

> M[1,1:7]
infra  vols   eco  crim   vma   vvi  camb
44.11 29.51  4.14  2.79  0.08  0.27  5.91 
> M4[1,1:6]
[1] 44.07698034 29.57490613  4.18441060   
     2.86142461  0.06670419  0.50994705

Si on regarde la distance entre ces deux matrices, on est effectivement relativement faible

> sum( (M4-M)^2 )
[1] 240.5894
> max(eigen(t(M4-M)%*%(M4-M))$value)
[1] 130.3316

cette dernière étant à comparer avec

> max(eigen(t(M)%*%M)$value)
[1] 447176.6

Cette matrice est effectivement de rang plus faible (de range 4)

> svd(M4)$d[1:8]
[1] 6.687127e+02 3.520604e+01 1.860897e+01 1.328171e+01 3.017803e-14 9.930232e-15 8.511710e-15 5.410929e-15

On peut aussi regarder une approximation de rang 2,

> k=2
> M2 = SVD$u[,1:k] %*% diag(SVD$d[1:k]) %*% t(SVD$v[,1:k])
> svd(M2)$d
[1] 6.687127e+02 3.520604e+01 1.012972e-13 9.873132e-15 6.152361e-15 4.691758e-15 3.289106e-15 1.206215e-15 1.248952e-16
> sum( (M2-M)^2 )
[1] 763.2868
> max(eigen(t(M2-M)%*%(M2-M))$value)
[1] 346.2936

On continuera pendant le cours.

Ah oui, et pour avoir davantage d’illustration, on pourra utiliser une autre base, afin de travailler sur des rangs (et non pas des valeurs réelles)

> BD = read.table("http://freakonometrics.free.fr/acp-villes.csv",sep=";",header=TRUE)
>  M=as.matrix(BD[,2:12])
>  colnames(M)=names(BD)[2:12]
>  rownames(M)=BD$Agglo