Présentation sommaire des systèmes de recommandation. Exemple jouet de filtrage collaboratif traité avec R par décomposition en valeurs singulières, factorisation non négative de matrice ou NMF et complétion de matrice; exemple réaliste de recommandation de films (MovieLens). Les mêmes données sont traitées de façon plus performante dans un autre calepin utilisant Spark.
Tous les sites marchands mettent en place des systèmes de recommandation pour déterminer les produits les plus susceptibles d'intéresser les internautes / clients en visite. Lorsque ceux-ci sont basées sur les seules informations concernant les interactions clients x produits, ils sont nommés filtrage collaboratif. Parmi ces derniers, les systèmes les plus aboutis sont basés sur la recherche d'un modèle de quelques facteurs latents susceptibles d'expliquer en faible dimension les interactions entre clients et produits.
D'autres systèmes sont basés sur des connaissances clients (user-based) ou sur des connaissances produits (item-based)ou encore sur des approches mixtes. Ils ne sont pas abordés ici.
Les données se mettent sous la forme d'une matrice X, toujours très creuse, contenant pour chaque client i (ligne) le nombre d'achats du produit j (colonne) ou une note d'appréciation de 1 à 5 lorsqu'il s'agit de films (Netflix), musiques (itune), livres...
Attention, la valeur "0" a du sens lorsqu'il s'agit d'un nombre d'achats alors qu'elle doit signifier une donnée manquante dans le cas d'une notation. Dans le premier cas, (valeur 0) le filtrage collaboratif est obtenu par une factorisation de matrice. Dans le 2ème (0 données manquante), il s'agit d'une complétion.
La décomposition en valeurs singulières (SVD) d'une matrice ainsi que la Non Negativ Matrix Factorization (NMF) sont utilisées dans ce contexte pour rechercher les facteurs (matrices W et H) reconstruisant au mieux la matrice X # W.H avec une contrainte de parcimonie ou faible rang sur les matrices W et H. Contrairement à la SVD où les facteurs sont recherchés orthogonaux 2 à 2, la NMF impose la contrainte de non négativité des matrices pour construire les facteurs de la décomposition. Ces facteurs ne permettent plus de représentation comme en ACP ou en MDS mais au moins une classification non supervisée tant des objets lignes que des objets colonnes de la matrice initiale. Ces classifications sont respectivement basées sur les matrices W et H des facteurs dits latents.
Schématiquement, $w_{ij}$ dénote l'appétence du i-ème utilisateur pour le j-ème facteur latent, tandis que $h_{jk}$ décrit quelle part du k-ième item intervient dans le j-ème facteur latent; le modèle suppose que la note $x_{ik}$ est la somme, sur tous les facteurs latents j, des produits $w_{ij}\times h_{jk}$.
La SVD, est basée sur un critère de moindre carrés norme trace des matrices. La librairie NMF (Gaujoux et Seoighe, 2010)\ de R propose plusieurs algorithmes de factorisation non négative, principalement: Multiplicative update algorithms et Alternate least Square (ALS), adaptés à deux fonctions possibles de perte: divergence de Kullback-Leibler (KL) ou moindres carrés (norme trace).
Attention, les choix d'option: fonction objectif, algorithme, rang des matrices, influencent fortement les résultats obtenus et ce d'autant plus que les algorithmes (NMF) convergent (au mieux) vers des optimums locaux. La SVD bénéficie d'une convergence "globale" mais est moins adaptée au contexte car les solutions sont moins cohérentes avec l'objectif recherché: des notes ou comptages nécessairement positifs.
Lorsque, les données sont des notes d'appréciation, la valeur "0" signifie en principe une valeur manquante. L'usage de la NMF ou de la SVD est alors abusif. Cette situation a été largement popularisée avec le concours Netflix à 1M$. Il a été abordé de façon théorique par Candes et Tao (2010) comme un problème de complétion de matrice sous contrainte de parcimonie; problème difficile, dont de très nombreuses approximations et implémentations ont depuis été proposées. Une simple utilisation est proposée ici dont l'algorithme (Mazumder et al. 2010) implémenté dans R (softImpute) conduit également à une factorisation.
Des données fictives triviales sont testées afin d'illustrer la démarche. Elles contiennent des nombres d'achats de certains produits ou des notes d'appréciation et peuvent être complétées à loisir au gré de votre imagination.
Le fichier fictif recom-jouet.dat est dans le même dépôt que ce tutoriel. Il contient un nombre d'achats par client.
In [ ]:
jouet=read.table("recom-jouet.dat")
jouet
In [ ]:
options(repr.plot.width=6, repr.plot.height=3)
boxplot(jouet)
Les données sont bien creuses mais les variables s'expriment dans des unités et donc avec des variances très différentes. Une forme de normalisation peut s'avérer nécessaire. Elle concerne à la fois les produits (colonnes ou variables), car certains (chocolat) sont plébiscités plus que d'autres, ainsi que les clients qui peuvent avoir des échelles de notation ou des profiles de consommation très différents .
Chargement de la librairie et identification des algorithmes disponibles. Plusieurs initialisation sont possibles; seule celle aléatoire par défaut est utilisée.
In [ ]:
library(NMF)
nmfAlgorithm()
nmfAlgorithm("brunet")
nmfAlgorithm("lee")
nmfAlgorithm("snmf/l")
nmfAlgorithm("snmf/r")
Q Identifier la fonction perte; les deux derniers algorithmes sont issus de l'ALS.
Comparer les méthodes en exécutant pour chacune d'entre elles 10 factorisations de rang 5. Les exécutions sont répétées car la convergence locale dépend de l'initialisation.
In [ ]:
res.multi.method=nmf(jouet, 5,nrun=10, list("brunet","lee","snmf/l","snmf/r"),
seed = 111, .options ="t")
compare(res.multi.method)
In [ ]:
options(repr.plot.width=6, repr.plot.height=5)
consensusmap(res.multi.method,hclustfun="ward")
Q Choisir la méthode.
In [ ]:
estim.r=nmf(jouet,2:6,method="snmf/l", nrun=10,seed=111)
plot(estim.r)
In [ ]:
options(repr.plot.width=8, repr.plot.height=5)
consensusmap(estim.r)
Q Choix du rang des matrices de la décomposition?
Une fois méthode et rang déterminés, itérer plusieurs fois l'exécution pour retenir la "meilleure" puis extraction des "facteurs".
In [ ]:
nmf.jouet=nmf(jouet,4,method="snmf/l",nrun=30,seed=111)
w=basis(nmf.jouet)
h=coef(nmf.jouet)
Production de classifications non-supervisées et graphiques associés aux matrices w et h de la factorisation. Ceci permet d'identifier des groupes de clients au regard de leur consommation ou préférences comme de construire des classes de produits appréciés simultanément.
In [ ]:
options(repr.plot.width=4, repr.plot.height=3)
basismap(nmf.jouet,hclustfun="ward")
In [ ]:
options(repr.plot.width=5, repr.plot.height=3)
coefmap(nmf.jouet,hclustfun="ward")
Comme c'est logique, le dendrogramme produit dans les cartes précédentes est directement issu des classifications ascendantes hiérarchiques calculées à partir des distances euclidiennes entre les lignes de w d'une part et les colonnes de h d'autre part.
La classification des objets est représentable dans les coordonnées d'un MDS ou dans les composantes d'une ACP des "facteurs" de la NMF; c'est équivalent en considérant la distance euclidienne définie à partir de ces facteurs.
Les produits à plus forte occurrence ou note peuvent prendre trop d'importance, les facteurs sont réduits.
In [ ]:
distmod.h=dist(scale(t(h)), method="eucl")
mdjouet= cmdscale(distmod.h, k=2)
hclusmod.h=hclust(distmod.h,method="ward.D")
options(repr.plot.width=5, repr.plot.height=4)
plot(hclusmod.h)
In [ ]:
dN.h=dimnames(h)[[2]]
hclasmod.h = cutree(hclusmod.h,k=4)
plot(mdjouet, type="n", xlab="", ylab="",main="")
text(mdjouet,dN.h,col=hclasmod.h)
abline(v=0,h=0)
In [ ]:
distmod.v=dist(scale(w), method="eucl")
mdjouet= cmdscale(distmod.v, k=2)
hclusmod.v=hclust(distmod.v,method="ward.D")
plot(hclusmod.v)
In [ ]:
hclasmod.v = cutree(hclusmod.v,k=2)
dN.v=dimnames(w)[[1]]
plot(mdjouet, type="n", xlab="", ylab="",main="")
text(mdjouet,dN.v,col=hclasmod.v)
abline(v=0,h=0)
Il n'est pas possible comme en ACP ou AFCM de mettre en relation les deux représentations des lignes et colonnes, individus et variables de la matrice factorisée. Cela peut être fait de façon détournée à l'aide d'une heatmap qui intègre les classifications obtenues en réordonnant les lignes et colonnes de X.
In [ ]:
# intégration des deux classifications
aheatmap(jouet,Rowv=hclusmod.v, Colv=hclusmod.h,annRow=as.factor(hclasmod.v),
annCol=as.factor(hclasmod.h))
L'objectif de rechercher les produits les plus susceptibles d'intéresser les clients. Celui-ci est atteint en reconstruisant une approximation de la matrice x par produit des matrices des facteurs.
Les couples (client, produit) pour lesquels les valeurs reconstruites sont le plus élevées alors qu'il n'y a pas eu d'achat ou de notation, sont ceux qui sont ciblés afin de proposer le produit identifié au client de ce couple.
Attention, le choix du rang est déterminant. En utilisant le choix optimal précédent (r=4) la reconstruction est finalement "trop" bonne et aucune recommandation n'émerge de la reconstruction de x. Le choix r=2, sous-optimal, fait ressortir des couples candidats.
In [ ]:
# Exécution avec r=2
nmf.jouet=nmf(jouet,2,method="snmf/l", nrun=30,seed=111)
# Matrice reconstruite
xchap=w%*%h
Comparer avec les données initiales, identifier le plus fort score reconstruit par client et identifier le produit correspondant.
In [ ]:
prod=apply(xchap-10*jouet,1,function(x) which.max(x))
cbind(dN.v,dN.h[prod])
Remarques
La décomposition en valeurs singulières propose également une factorisation de la matrice X=UL V'. Celle-ci a de bien meilleure propriétés numériques dont l'unicité de la solution optimale qui est atteinte pour un rang fixé. Des contraintes de parcimonie ou de régularité peuvent également être associées à la fonction perte quadratique (sparse SVD).
La démarche est équivalente et découle directement de la SVD donc de l'analyse en composantes principales de X. Remarquer une forte similitude entre les représentions obtenus par MDS des facteurs de le NMF et celles de l'ACP:
In [ ]:
biplot(prcomp(jouet,scale=TRUE))
In [ ]:
res=svd(jouet)
# Matrice reconstruite
xchap=res$u[,1:2]%*%diag(res$d[1:2])%*%t(res$v[,1:2])
Comparer avec les données initiales, identifier le plus fort score reconstruit par client, identifier le produit correspondant.
In [ ]:
prod=apply(xchap-10*jouet,1,function(x) which.max(x))
cbind(dN.v,dN.h[prod])
Q Comparer aec les recommandations précédentes sur cet exemple trivial.
On considère que les données sont des notes d'appréciation: la valeur "0" signifie en principe une valeur manquante et le problème est formellement celui d'une complétion de matrice. La librairie softImpute de R en propose une solution par un algorithme de SVD seuillée (Hastie et al. 2010). Il s'agit donc d'approcher une matrice très creuse avec beaucoup de valeurs manquantes par une matrice de faible rang. L'algorithme est analogue à un algorithme EM pour imputation de données manquantes. La fonction R accepte la classe de représentation des grandes matrices creuses.
Les valeurs nulles sont remplacées par des valeurs manquantes.
In [ ]:
jouet.na=jouet
jouet.na[jouet==0]=NA
La recommandation se fait comme avec la SVD. Une étude plus approfondie de cet algorithme et de son usage s'avère nécessaire, notamment pour régler rang et paramètre de pénalisation.
In [ ]:
install.packages("softImpute")
library(softImpute)
res=softImpute(jouet.na,rank.max=2,type="svd",lambda=1)
# Matrice reconstruite
xchap=res$u%*%diag(res$d)%*%t(res$v)
Comparer avec les données initiales, identifier le plus fort score reconstruit par client, identifier le produit correspondant
In [ ]:
prod=apply(xchap-10*jouet,1,function(x) which.max(x))
cbind(dN.v,dN.h[prod])
Q Comparer les résultats.
Cette étude n'est qu'une brève introduction au problème du filtrage collaboratif. De nombreuses questions n'ont pas été abordées dont le cold start et surtout celle très importante de l'évaluation de tels systèmes. Un procédé simple consiste à volontairement supprimer des valeurs afin de voir si le système les retrouve et en quelle proportion; c'était le principe du concours Netflix retenue dans la suite.
Le filtrage collaboratif est appliqué aux données publiques du site GroupLens. L'objectif est de tester les méthodes et la procédure d'optimisation sur le plus petit jeu de données composé de 100k notes de 943 clients sur 1682 films où chaque client a au moins noté 20 films. Les jeux de données plus gros (1M, 10M, 20M notes) peuvent être utilisés pour "passer à l'échelle volume".
Les données initiales sont sous la forme d'une matrice très creuse (sparse) contenant des notes ou évaluations. Les0 de la matrice ne sont pas des notes mais des données manquantes, le film n'a pas encore été vu ou évalué.
La complétion de grande matrice creuse est obtnue par l'algorithme implémenté dans la librairie softImpute de R. Un autre calepin utilise la version de NMF de MLlib de Spark qui permet également la complétion.
En revanche,la version de NMF incluse dans la librairie Scikit-learn traite également des matrices creuses mais le critère (moindres carrés) optimisé considère les "0" comme des notes nulles, pas comme des données manquantes.
Dans la première partie, le plus petit fichier est partagé en trois échantillons: apprentissage, validation et test; l'optimisation du rang de la factorisation (nombre de facteurs latents) est réalisée par minimisation de l'erreur estimée sur l'échantillon de validation.
Ensuite le plus gros fichier est utilisé pour évaluer l'impact de la taille de la base d'apprentissage.
Néanmoins, comme les résultats obtenus ne concurrencent pas ceux de MLlib et sont même assez décevants (Besse et al. 2016), ce calepin ne fait que décrire a minima le code de mise en oeuvre de la librairie sur les données MovieLens.
Les fichiers de données sont dans le même dépôt que ce calepin.
In [ ]:
library(softImpute)
# lecture
dBrut=read.csv("ml-ratings100k.csv")
In [ ]:
# extraction d'un échantillon test et d'un échantillon d'apprentissage
dTestInd=sample(nrow(dBrut),nrow(dBrut)/10,replace=FALSE)
dTest=dBrut[dTestInd,1:3]
dTrain=dBrut[-dTestInd,1:3]
In [ ]:
# Mise au format d'une matrice creuse
dTrainSparse=Incomplete(dTrain$userId,dTrain$movieId,dTrain$rating)
# appel de la fonction
res=softImpute(dTrainSparse,rank.max=4,type="als",lambda=1,maxit=200)
# complétion
recom=impute(res,dTest[,1],dTest[,2])
In [ ]:
# Calcul de l'erreur (RMSE)
sqrt(mean((dTest[,3]-recom)**2))
In [ ]:
# Lecture de la matrice
dBrut=read.csv("ratings20M.csv")
In [ ]:
# Extraction de l'échantillon test
dTestInd=sample(nrow(dBrut),nrow(dBrut)/10,replace=FALSE)
dTest=dBrut[dTestInd,1:3]
nrow(dTest)
In [ ]:
# sous-échantillonage de l'échantillon d'apprentissage
dInter=dBrut[-dTestInd,1:3]
taux=0.1
dTrainInd=sample(nrow(dInter),nrow(dInter)*taux,replace=FALSE)
dTrain=dInter[dTrainInd,1:3]
# Matrice d'échantillonnage sparse
dTrainSparse=Incomplete(dTrain$userId,dTrain$movieId,dTrain$rating)
Attention. La démarche adoptée ici pour trouver les meilleurs valeurs des paramètres (rang, pénalisation) est abusive. Les valeurs ont été optimisées sur l'échantillon test sans pour autant atteindre des résultats concurrentiels avec la solution utilisant MLlib de Spark. Voir la comparaison détaillée par Besse et al. 2016.
In [ ]:
# Factorisation
t1=Sys.time()
res=softImpute(dTrainSparse,rank.max=10,type="als",lambda=20,maxit=200)
t2=Sys.time()
# Reconstruction
recom=impute(res,dTest[,1],dTest[,2])
#RMSE
sqrt(mean((dTest[,3]-recom)**2))
difftime(t2,t1)
Apprentissage avec fichier complet (taux=1): 20M de notes
| rang | lambda | temps (') | rmse |
|---|---|---|---|
| 4 | 1 | 5.6 | 1.07 |
| 10 | 10 | 12.6 | 1.02 |
| 10 | 20 | 12.2 | 1.033 |
| 15 | 10 | 19.4 | 1.016 |
| 20 | 1 | 26.9 | 1.02 |
| 20 | 10 | 26.1 | 1.016 |
| 20 | 15 | 24.4 | 1.018 |
| 20 | 20 | 27.0 | 1.016 |
| 30 | 20 | 40.1 | 1.02 |