Guillaume

Benoît Lou

Marie-Luce Asma

我还要特别感谢我的父母 。 他们不仅辛勤努力

Yohann De Castro

Yannig

Remerciement

Je remercie avant tout mes encadrants :

ii As in manufacture so in science-retooling is an extravagance to be reserved for the occasion that demands it. The significance of crises is the indication they provide that an occasion for retooling has arrived.

Thomas Kunh The Structure of Scientific Revolutions

2.1

The matrix representation of the variable of interest. Typically, an entry of the matrix is the electricity consumption of an individual (a column), during a period (a row). Measures taken by sensors on the network are often aggregations of individuals (sum of the entries in the green box). Measures read on user meters are often temporally aggregated (sum of the entries in the orange box). . . . . .

6.2

Mean RRMSE of the recovered matrices over three separate runs over the four datasets. On the samples with random observation periods, proposed methods (HALS and NeNMF, blue and purple lines, both penalized and unpenalized) outperforms the interpolation benchmark (solid red line). On the samples with regular observation periods, unpenalized HALS and NeNMF (solid blue and purple lines) are similar to the interpolation benchmark, whiled penalized HALS and NeNMF (dashed blue and purple lines) are an important improvement. The softImpute method (solid green line) only has comparable performance in two of the datasets, in the easiest task (50% sampling rate at random periods). In most cases, RRMSE of softImpute is larger than 100%. . . . . . . . . . . . . . .

6.3

The rank chosen by the cross validation procedure generally increases with the sampling rate, for the four datasets. This shows that the procedure is able to relax the rank constraint when more information is available in the data. . . . . . L'électricité est un aspect important des activités socio-économiques. Ainsi les données de consommation électrique sont une information intéressante en soi. Elles peuvent être mises en parallèle avec d'autres statistiques publiques sur la société. Si la consommation électrique est connue à la même échelle spatio-temporelle que d'autres informations socio-démographiques, croiser ces données peut permettre • d'une part, aux entreprises d'électricité de fournir un service de manière plus efficace et stable ;

• d'autre part, aux chercheurs en sciences sociales de révéler comment l'utilisation d'électricité influence les divers phénomènes sociaux, et vice versa.

Cependant, les données sont très souvent insuffisantes pour produire des estimations et prévisions à l'échelle spatio-temporelle exigée par les applications mentionnées. Deux sources de données sont souvent disponibles aux entreprises d'électricité :

• des mesures de capteurs installés sur le réseau ;

• des mesures venant des compteurs de clients.

Les deux sources de données ont des limitations importantes. Les données de capteurs sont souvent de fréquence temporelle assez haute (pas de temps de 10 ou 15 minutes), mais elles couvrent un nombre important de clients. En France, le dernier niveau du réseau de distribution équipé de capteurs de manière systématique est celui des départs Haute Tension A (HTA). Chaque départ HTA couvre en moyenne 1500 clients. D'autre part, les compteurs clients ne couvrent en général qu'un seul client. Par contre, la plupart des clients sont pour l'instant équipés de compteurs traditionnels qu'il faut relever de manière manuelle : un technicien doit se déplacer physiquement chez le client pour lire l'index sur le compteur. Ainsi, ces index sont souvent relevés une fois tous les quelques mois. Avec le déploiement de compteurs intelligents (Linky), les index peuvent être stockés localement toutes les 10 minutes. Néanmoins, à cause des coûts de transmission et de traitement, ainsi que des considérations de protection de la vie privée, ils ne seront accessibles aux entreprises d'électricité que sous forme de sommes journalières.

Dans cette thèse, on envisage des méthodes statistiques pour produire des estimations et prévisions sur une échelle spatio-temporelle assez fine. L'objectif est de développer des méthodes qui peuvent utiliser toute information disponible, pour avoir des estimations de qualité supérieure. À part les deux sources de données de consommation, les méthodes développées ici utilisent notamment des "informations auxiliaires", ou side information en anglais, pour enrichir les estimations. Ces informations auxiliaires peuvent être globales, comme l'autocorrélation temporelle des courbes de consommation. Elles peuvent aussi être locales : des variables exogènes correspondant à une période et un individu dans le jeu de données, comme la température d'une période ou l'information socio-démographique d'un client.

FORMALISATION DU PROBLÈME

Formalisation du problème

Dans cette thèse, on s'intéresse à la valeur d'une variable d'intérêt v ∈ R (typiquement la consommation électrique), pendant n 1 périodes consécutives, pour un groupe de n 2 individus1 . Pour un individu j, j ∈ {1, 2, ..., n 2 }, v i,j représente la consommation pendant la période i, pour i ∈ {1, 2, ..., n 1 }. On utilise des lettres majuscules en gras pour représenter des matrices et des lettres minuscules en gras pour représenter des vecteurs. Ainsi, V ∈ R n 1 ×n 2 est la matrice regroupant la variable d'intérêt, (v i ) T est le i-ième ligne de V (la consommation de tous les individus à la période i), et v j est la j-ième colonne (la consommation de l'individu j à travers les périodes).

On peut aisément formaliser le problème d'estimation mentionné ci-dessus avec cette représentation matricielle. La matrice représentant la variable d'intérêt est le tableau bleu dans figure 1

: Représentation matricielle de la variable d'intérêt. Dans cette thèse, une entrée de la matrice est souvent la consommation électrique d'un individu (une colonne), durant une période (une ligne). Les mesures de capteurs sur le réseau sont souvent des agrégations d'individus (somme d'entrées dans le rectangle vert). Les relevés de compteurs individuels sont souvent des agrégations temporelles (somme d'entrées dans le rectangle orange).

Les variables exogènes sont représentées par les deux tableaux verts. Leurs lignes et colonnes sont alignées avec celles du tableau bleu.

Les objectifs d'estimation ont une représentation naturelle dans cette formalisation matricielle. Dans la plupart du temps, on est intéressé par la valeur de la variable d'intérêt pour chaque période de chaque utilisateur. Ainsi, il faut estimer toutes les entrées de V. Selon des applications, on peut être intéressé par une agrégation des entrées de V. Par exemple, dans le problème de changement de support spatial, on est amené à estimer la variable d'intérêt sur une agrégation spatiale (entre individus) différente de celle des données. On peut aussi s'intéresser à la valeur maximale de la consommation électrique quotidienne, ce qui correspond à la valeur maximale de chaque colonne pendant un certain nombre de périodes.

La prévision correspond à estimer une ou plusieurs nouvelles lignes (nouvelles périodes) ou colonnes (nouvelles individus) de V dans laquelle aucune entrée n'a été couverte par les données. Quand aucune information à part les mesures de V n'est disponible, il est difficile de produire des prédictions. On considère deux cas de variables exogènes dans cette thèse :

• information implicite sur V, telle que les colonnes de V ont une structure auto-régressive, soit entre périodes temporelles, soit entre des individus qui ont une structure de corrélation spatiale;

• information explicite, comme par exemple la température de chaque période, ou des informations socio-démographiques sur les individus.

Cette thèse est organisée de la manière suivante : dans le reste de ce chapitre, les trois approches envisagées dans la thèse seront présentées brièvement. L'approche basée sur la factorisation de matrice, qui est la proposition principale de la thèse, sera introduite de manière un peu plus détaillée. À part ce chapitre, le reste de cette thèse sera en anglais. Le chapitre 2 reprendra le contenu de ce chapitre et fera une introduction en anglais. Les chapitres 3 et 4 forment la première partie, où l'on considère des données qui sont des agrégations spatiales ou des agrégations d'individus, et dont le but est l'estimation à des niveaux d'agrégation différentes de celles des données. Dans la deuxième partie, les chapitres 5-9, on présente l'approche basée sur la factorisation de matrice nonnégative, qui permet de résoudre le problème de manière plus systématique. Le lien entre les trois approches est résumé dans le chapitre 8 et les perspectives sont présentées dans le chapitre 10.

Solutions classiques et nouvelles 1.3.1 Statistique spatio-temporelle

Quand les positions géographiques des individus sont disponibles, un choix naturel est d'utiliser des outils de la statistique spatiale. En considérant des valeurs individuelles comme observations d'un processus aléatoire spatial ou spatio-temporel, cette classe de méthodes fournit des modèles qui sont adaptés pour des phénomènes assez divers. Dans le chapitre 3, on présentera la statistique spatiale. En particulier, on discutera du krigeage (kriging en anglais), une méthode qui permet d'interpoler et d'extrapoler la variable d'intérêt à des positions géographiques nouvelles. En modélisant la structure de corrélation du processus spatial, cette méthode peut fournir des estimations optimales. On appliquera cet outil à deux jeux de données, pour examiner sa capacité à résoudre des problèmes représentés sur la figure 1.1.

Prédiction de consommation électrique multiple

Dans une approche parallèle à celle de la statistique spatiale, on présentera une méthode qui modélise plutôt la structure de la moyenne de la variable d'intérêt dans le chapitre 4. Cette approche permet d'estimer la consommation d'électricité des zones cibles, où l'information sociodémographique est disponible, à partir des zones sources, où la consommation est disponible à une échelle temporelle petite. Cette méthode utilise le clustering pour relier ces deux types 1.4. CONTRIBUTIONS À LA FACTORISATION DE MATRICE NONNÉGATIVE de zones. Le transfert des modes de consommation électrique est réalisé par des modèles de régression utilisant des variables explicatives telles que la température ou le cycle annuel, hebdomadaire et journalier. Cette méthode est évaluée par sa capacité de prédiction sur deux jeux de données synthétiques et réelles.

Méthodes basées sur la factorisation de matrice

Comme suggéré par la figure 1.1, l'objet d'intérêt de cette thèse se prête bien à une représentation matricielle. Du chapitre 5 au chapitre 9, on propose une approche générique basée sur la factorisation de matrice nonnégative [4][5][6] qui résout le problème d'estimation. Nous allons présenter les intuitions derrière cette approche et la contribution de cette thèse sur ces méthodes dans la section suivante.

Contributions à la factorisation de matrice nonnégative 1.4.1 Introduction à la factorisation de matrice nonnégative

Récemment, la factorisation de matrice nonnégative (nonnegative matrix factorization en anglais, ou NMF) a connu un succès important dans des domaines très divers, de la segmentation des images [4] au système de recommandation [START_REF] Koren | Matrix Factorization Techniques for Recommender Systems[END_REF] en passant par le clustering des documents [START_REF] Xu | Document Clustering Based on Non-Negative Matrix Factorization[END_REF] et la décomposition du son [START_REF] Févotte | Itakura-Saito Nonnegative Factorizations of the Power Spectrogram for Music Signal Decomposition[END_REF]. L'idée principale derrière cette classe de méthodes est la suivante. Étant donné V ∈ R n 1 ×n 2 + , une matrice de grande taille, on peut souvent faire l'hypothèse que son rang k est faible devant sa dimension. Les colonnes de V sont ainsi à l'intérieur d'un sous-espace de R n 1 , de dimension k, bien inférieure à n 1 . Ceci équivaut à dire que V = F r F T c , avec

F r ∈ R n 1 ×k + , F c ∈ R n 2 ×k

+

. Les matrices F r et F c sont appelés les facteurs.

Le faible rang du problème a une interprétation très intuitive : dans le clustering de documents, le nombre de thèmes est souvent bien inférieur au nombre de documents et au nombre de mots dans le dictionnaire ; dans un système de recommandation, il y a peu de clients et marchandises "indépendants" parmi un grand nombre de clients et marchandises. Le choix d'imposer les contraintes de non-négativité est souvent justifié par les applications aussi : avec des facteurs, on peut aisément interpréter les documents comme un mélange de thèmes, et une marchandise comme un mélange de plusieurs types de marchandises.

Une interprétation du même type pour les applications dans le domaine d'électricité nous semble aussi naturelle, puisqu'on remarque qu'il y a un nombre réduit de manières d'utiliser l'électricité (appelé les profils), bien qu'il y ait un nombre important de clients. Tous les clients sont proches d'un mélange de ces profils. Cet aspect est d'ailleurs pris en compte dans des systèmes actuels de reconstruction de consommation électrique, qui utilisent des profils nationaux de consommation [10]. Ces profils sont actuellement calibrés par des méthodes statistiques qui nécessitent des traitements assez longs. Utiliser des méthodes basées sur la factorisation de matrice nonnégative permettrait de les estimer de manière beaucoup plus automatique, à partir des données partielles.

Pour obtenir une factorisation, on peut résoudre le problème d'optimisation suivant, min Fr∈R n 1 ×k ,Fc∈R n 2 ×k (V, F r F T c ) s.t. F r ≥ 0, F c ≥ 0, où est une fonction de perte qui mesure la différence entre F r F T c et V, et les inégalités sont CHAPTER 1. RÉSUMÉ EN FRANÇAIS élément par élément. La fonction peut être assez générale, mais on considère principalement la perte quadratique en norme Frobenius où (X, Y) = X -Y 2 F = i,j (x i,j -y i,j ) 2 . La factorisation de matrice n'est pas un nouveau problème. Dans le cas de l'Analyse en Composantes Principales (ACP), une factorisation de matrice est aussi calculée dans laquelle un des facteurs est contraint à être une matrice orthogonale (celui des composantes principales).

Malgré son interprétation très intuitive, la NMF n'est pas facile à résoudre. En fait, dans le cas général, la NMF est NP-difficile [11].

Une ligne de recherche théorique consiste à étudier si une matrice avec des entrées nonnégatives a une NMF "unique". Pour toute paire (F r , F c ), une NMF de V, et pour toute matrice inversible R ∈ R k×k , (F r R, F c (R -1 ) T ) est encore une factorisation de V (leur produit est égal à V). Mais cette nouvelle factorisation n'est pas nécessairement nonnégative. Par contre, pour certaines matrices, il est possible que les seules matrices inversibles R telles que (F r R, F c (R -1 ) T ) reste une factorisation nonnégative soient des compositions de dilatation et de permutation. Dans ce cas-là, on dit que la NMF est "unique". Décider si une matrice a une NMF unique est aussi un problème difficile. Deux conditions nécessaires et suffisantes de l'unicité ont été trouvées, mais étant donnée V, il est difficile de décider si elle vérifie ces conditions [5,12]. Des conditions suffisantes, plus faciles à vérifier, ont été proposées sur des classes particulières de NMF [12,[START_REF] Huang | Non-Negative Matrix Factorization Revisited: Uniqueness and Algorithm for Symmetric Decomposition[END_REF].

Des algorithmes avec garantie de convergence ont été proposés pour certaines classes de NMF, en particulier celles qui vérifient les conditions suffisantes d'unicité [START_REF] Arora | Computing a Nonnegative Matrix Factorization-provably[END_REF][START_REF] Gillis | Sparse and Unique Nonnegative Matrix Factorization through Data Preprocessing[END_REF][16][17]. Ces algorithmes sont souvent assez spécifiques, et la preuve de convergence utilise des outils complexes. Pourtant, en pratique, des algorithmes dits de "premier ordre" sont préférés. Ces algorithmes sont souvent basés sur la descente de gradient ou l'algorithme de descente par bloc. Même s'ils ont peu de garantie théorique, la simplicité d'implémentation et la performance empirique ont fait qu'ils sont devenus très populaires [18]. L'idée principale de ces algorithmes de descente est la méthode Gauss-Seidel : soit une fonction f (x 1 , ..., x d ) à minimiser sur un produit cartésien de convexe X 1 × ... × X d . L'algorithme minimise par rapport à chaque bloc de variables x 1 , ..., x d de manière itérative. Dans le cas de NMF, les algorithmes sont caractérisés par les blocs de variables : on peut utiliser des blocs scalaires et alterner entres les éléments de F r et F c , des blocs vectoriels (les colonnes et lignes de F r et F c ), ou des blocs matriciels (F r et F c ).

Grâce aux résultats d'optimisation, les algorithmes du type Gauss-Seidel convergent vers un point stationnaire sous des conditions faibles. Ainsi la plupart des algorithmes de NMF convergent vers un point stationnaire.

Appliquer la NMF pour estimer la matrice V (chapitres 5, 6)

La contribution principale de la thèse est l'application de la NMF au problème introduit dans la section 1.2. Sous l'hypothèse de faible rang, on peut formaliser le problème d'estimation de la figure 1.1 comme un problème similaire à la NMF. Dans ce problème, la matrice V est à estimer, à partir de mesures linéaires α : min (1.1) où A est un opérateur linéaire sur l'espace des matrices de taille n 1 × n 2 et α est un vecteur de 1.4. CONTRIBUTIONS À LA FACTORISATION DE MATRICE NONNÉGATIVE données de dimension N . Les opérateurs linéaires sur R n 1 ×n 2 sont caractérisé par N matrices de taille n 1 × n 2 , A 1 , ..., A N , et

Fr∈R n 1 ×k ,Fc∈R n 2 ×k ,V∈R n 1 ×n 2 V -F r F T c 2 F s.t. F r ≥ 0, F c ≥ 0, V ≥ 0, A(V) = α,
α i = Tr(VA T i ) = V, A i .
Ces matrices A 1 , ..., A N sont appelés des masques.

Dans le chapitre 6, nous proposons une extension de la NMF pour utiliser les agrégations temporelles comme données.

Ceci correspond à expliciter (1.1) dans le cadre où les mesures linéaires sont des agrégations temporelles. Concrètement, ceci veut dire que pour tout 1 ≤ i ≤ N ,

α i = V, A i = (t,n)∈I i v t,n ,
où I i est l'ensemble des indices des périodes couvertes par la i-ième mesure. Ces contraintes de somme, et la contrainte de nonnégativité, définissent la région admissible de V, qui est un simplex.

Dans le chapitre 6, nous proposons un algorithme qui résout (1.1). Pour ce faire, nous ajoutons, aux algorithmes classiques de NMF, une étape de projection dans ce simplex défini par des contraintes (Algorithm 1). Cet algorithme combine un algorithme efficace de projection dans un simplex, et deux algorithmes classiques de NMF [19][20][21]. C'est un algorithme convergent vers un point stationnaire comme la plupart d'algorithmes de NMF.

En plus de cet algorithme, nous proposons aussi de prendre en compte deux types d'informations supplémentaires.

Autocorrélation temporelle des individus (chapitre 6)

Le premier type d'information supplémentaire est l'autocorrélation temporelle des individus. En pratique, il est connu que les séries temporelles de consommation électrique sont souvent assez régulières au sens où il existe une certaine dépendance temporelle entre les observations. Nous prenons en compte cet aspect dans chapitre 6.

La régularité des facteurs a été envisagée dans de nombreuses études de factorisation de matrice ou de NMF. Souvent, une pénalisation est ajoutée à la fonction objectif pour obtenir des facteurs plus réguliers et/ou plus sparses [22][23][24][25][26]. Ici, nous avons un objectif différent. Nous supposons avoir de l'information a priori sur l'autocorrélation temporelle des individus, et essayons d'améliorer la reconstitution grâce à cette information supplémentaire. Cette information peut venir de sources différentes : elle peut être le résultat des observations historiques (sur des individus qui ont participé à des études expérimentales, par exemple), ou elle peut être issue de modélisation.

Nous modifions (1.1) en ajoutant un terme de pénalisation lié aux seuils d'autocorrélation. Nous montrons qu'ajouter la pénalisation est dans certains cas équivalent à contraindre que l'autocorrélation de chaque individu soit supérieure à ces seuils (Theorem 1). Nous proposons un algorithme pour ce problème modifié (Algorithm 2).

Prendre en compte des informations exogènes (chapitre 7)

Le deuxième type d'information supplémentaire que nous proposons de prendre en compte concerne les variables exogènes. Dans ce cadre, comme mentionné dans figure 1.1, des variables exogènes liées correspondant aux lignes ou colonnes de la matrice V sont disponibles. CHAPTER 1. RÉSUMÉ EN FRANÇAIS Dans la littérature de la NMF, les variables exogènes sont un sujet qui attire beaucoup d'attention. Ce principe est applicable aux diverses applications, notamment en système de recommandation, où la préférence d'un client pour une marchandise dépend naturellement des attributs du client (âge, sexe, profession, etc.) et de la marchandise (prix, utilité, gamme, etc.). En ce qui concerne l'application en électricité, ce principe est également pertinent, car la consommation électrique dépend d'attributs liés à la période (température, nébulosité, heure de la journée, jour de la semaine), et au client (résidentiel ou industriel, composition du ménage, type de logement).

Nous résumons l'état de l'art de la littérature sur ce sujet dans la table 7.1, sur la page 85. Remarquons que pour les applications qu'on envisage, il est intéressant de pouvoir traiter des mesures linéaires comme observations et des modèles de régression non-paramétriques pour gérer les variables exogènes. Comme les méthodes existantes ne sont pas assez flexibles sur ces deux points, nous proposons un cadre général de NMF avec variables exogènes (cf. section 7.1.1).

Connaissant des variables exogènes liées aux colonnes X c et aux lignes X r , en plus des mesures linéaires α, nous proposons de résoudre min

V,fr∈F k r ,fc∈F k c V -(f r (X r )) + (f c (X c )) T + 2 F s.t. A(V) = α, V ≥ 0, (1.2)
où (•) + est la fonction qui remplace toute valeur négative par 0, et F r et F c sont des espaces fonctionnels dans lesquels nous cherchons des fonctions de régression.

En variant l'opérateur de mesure A, l'espace fonctionnel de régression, et les variables exogènes, cette modélisation peut être adaptée à des problèmes assez divers. En particulier, pour modéliser les relations non-linéaires entre les variables exogènes et les facteurs, nous utiliserons des splines de régression et des méthodes à noyaux dans l'application de ce modèle. Sur l'aspect théorique de ce modèle, nous déduisons une condition suffisante d'unicité de NMF dans le cadre de NMF avec variables exogènes (Theorem 5). En pratique, garantir que la NMF est unique permet de mieux interpréter les facteurs obtenus par le modèle.

Algorithme convergent vers un point stationnaire

Nous proposons un algorithme pour résoudre (1.2) (Algorithm 4). Similaire aux algorithmes 1 et 2, dans cet algorithme appelé HALSX (Hierarchical Alternating Least Squares with eXogeneous variables), nous ajoutons une étape supplémentaire dans HALS, qui est un algorithme classique de NMF [27]. Cette étape consiste à estimer des modèles de régression avec comme variables explicatives les variables exogènes.

Nous démontrons que cet algorithme converge vers un point stationnaire pour une grande classe de modèles de régression (Theorem 7).

Applications aux jeux de données synthétiques et réelles

Les méthodes proposées sont évaluées sur plusieurs jeux de données synthétiques et réelles.

Les jeux de données synthétiques sont créés pour tester la pertinence des méthodes proposées dans un environnement contrôlé. Cela consiste à simuler des matrices de rang faible de manière exacte. Les détails de simulation sont présentés dans les sections 6.3 et 7.4. • Consommation électrique portugaise (Jeu de données publiques [28]) la consommation journalière de 370 clients individuels (anonymisés) portugais. Issu d'un jeu de données publique, la consommation est disponible au pas de 15 minutes pendant 4 ans. Pour les expériences aux chapitres 6 et 7, la consommation est agrégée au pas journalier.

Jeux de données de consommation électrique

• Consommation électrique irlandaise (Jeu de données publiques [29,[START_REF]Results of Electricity Coast-Benefit Analysis, Customer Behaviour Trials and Technology Trials Commission for Energy Regulation[END_REF]) la consommation journalière de 426 petites et moyennes entreprises (PME) irlandaises. Issu d'une expérimentation européenne, la consommation est disponible au pas horaire pendant l'an 2010. L'agrégation journalière de ce jeu est utilisée au chapitre 6.

Jeux de données de système de recommandation Pour évaluer la performance de la méthode sur un opérateur de mesure linéaire autre que l'agrégation temporelle, on utilise aussi un jeu de données standard en système de recommandation.

• MovieLens C'est un jeu de données publique de scores de films. En tout, 100,000 scores sont disponibles pour 943 utilisateurs et 1682 films.

Implémentation

Les algorithmes NMF proposés qui ont été développés, sont regroupés dans un package R appelé meterModels. Pour faciliter les tests, des fonctions d'expérimentation, d'évaluation des méthodes, ainsi que des méthodes de référence sont aussi incluses dans ce package.

La plupart du package est en R. Certaines fonctions intensives en calcul sont écrites en C++ avec une interface en R, pour accélérer l'exécution des programmes. Quand la bibliothèque est disponible, le package peut utiliser openMP pour distribuer certaines calculs. Quand le package est utilisé sur une machine multi-coeurs, l'estimation de plusieurs modèles peut aussi se faire en parallèle en utilisant le package parallel de R.

Représentation des observations

Les mesures d'agrégation temporelle et les mesures linéaires générales sont gérées par le package. Pour les agrégations temporelles, les observations sont représentées comme des dates d'observation et des index de compteurs. Pour les mesures linéaires générales, l'opérateur de mesure est représenté comme une grande matrice sparse, et les observations sont le produit de l'opérateur et le vecteur des entrées de la matrice.

Méthodes disponibles Les méthodes de NMF avec agrégations temporelles ou mesures linéaires sont disponibles, ainsi que les variantes proposées par la thèse : la prise-en-compte de l'autocorrélation et des variables exogènes. Les méthodes proposées ainsi que les méthodes de référence sont détaillées dans la table 9.1 sur la page 118. L'estimation de modèle se fait à travers une interface unifiée, où le choix des méthodes se fait par un paramètre. Il est aussi possible de comparer plusieurs méthodes sur un jeu de données avec des paramètres variables.

Conclusions et perspectives

Dans cette thèse, nous avons étudié un problème pratique en industrie : la modélisation de la consommation électrique à partir des données partielles issues de nombreux individus. Nous avons envisagé ce problème sous différents angles, avec des sources de données agrégées temporellement ou individuellement, avec comme objectif l'estimation de la consommation passée, et la prévision de la consommation future. Nous avons investigué des solutions dans plusieurs cadres différents : la statistique spatio-temporelle, le clustering des individus, ainsi que la factorisation de matrice nonnégative. Particulièrement pour la NMF, nous avons adapté la méthode pour la problématique, et proposé des améliorations méthodologiques par rapport à l'état de l'art. La méthodologie a été évaluée par des expériences nombreuses sur des jeux de données synthétiques et réelles, et des résultats prometteurs ont été obtenus.

Des perspectives de recherche sont envisagées. Outre les applications directes de la méthodologie étudiées dans cette thèse, il existe de nombreuses applications industrielles. Nous pouvons utiliser la NMF pour estimer des statistiques importantes du réseau électrique. Par exemple, par rapport à la consommation électrique de chaque moment, la consommation maximale instantanée est cruciale pour la maintenance du réseau. Il serait intéressant d'évaluer la méthode de NMF dans ce contexte, voire même changer la fonction objectif du problème d'optimisation pour s'adapter à la statistique qui nous intéresse. Dans la thèse, nous avons fait l'hypothèse que le dispositif d'acquisition des données est fixe. Or, il peut être envisageable de modifier la manière dont les données sont collectées, et donc de trouver une manière optimale de collecter les données, sachant que la méthodologie de matrice de rang faible est disponible pour le traitement de ces données. Nous avons étudié des méthodes adaptées séparément pour les agrégations temporelles et individuelles. Il serait intéressant de voir si combiner les deux types de données peut encore améliorer la performance de la méthode. Sur l'aspect théorique de la NMF, des questions importantes restent également ouvertes, comme l'analyse de la convergence globale des algorithmes de premier ordre pour certaines classes de matrices, et l'approfondissement des conditions d'unicité de la NMF.

Motivation for electricity consumption estimation and prediction

To ensure the safety and stability of electric power systems, the electricity injected into a network has to be met with the load, namely the consumption. This is called load balancing. Electricity consumption is determined by the client mix connected to the network and external conditions, on which utility companies have very little control. As a core characteristic of electricity is that it can not be stored efficiently, utilities mostly achieve load balancing by adjusting power generation and operating the electric network according to real-time estimation and predictions of the load1 .

Given these conditions, it is easy to see how electricity consumption data can be useful for a utility company, such as Electricité de France (EDF). In order to plan power generation and network development, it is important to have mid-term or long-term consumption prediction, starting from six-month ahead. To plan network operations, it is useful to have temporally finegrained short or mid-term prediction and estimation, for example hourly or daily prediction until several months ahead. To take fully advantage of the burgeoning decentralized renewable energy generation in load balancing, it can be interesting to have consumption data and prediction on a fine spatial grid.

Apart from its value for power planning and network operations, electricity consumption data can be valuable for other reasons. In the context of open electricity market, regulations often stipulate that each market participant, i.e. suppliers, utility traders, large consumers, be responsible for its own load balancing at all time. Any imbalance caused by a participant is typically billed by transmission system operators (TSO). To calculate the imbalance, TSOs have to produce an hourly or half-hourly estimation of the consumption of each participant.

Electricity is an important part of socio-economic activities, and thus can be put into connection with many other public datasets concerning the society. Should electricity consumption CHAPTER 2. INTRODUCTION data be available at the same spatio-temporal scale as other socio-demographic information, utility suppliers can provide energy service more efficiently and stably, and researchers in social science can shed new lights on various social phenomena by crossing electricity data with data from other domain.

However, there is rarely enough data to produce estimations and predictions that are spatiotemporally fine-grained enough for these applications. Two sources of consumption data are often available to utility companies:

• measures taken by sensors installed on the network;

• measures read on end users' electric meters.

Both data sources are somewhat limited. Sensor data on the electric network are of sufficiently high frequency (at an interval of ten or fifteen minutes), but covers a rather big area. In France, for example, the lowest level in the electric grid in which sensor data is available is at the mediumvoltage (MV) feeder level, where each MV feeder covers approximately 1500 clients. Individual electric meter, on the other hand, only covers one client. However, traditional meters typically can only be read manually by a worker who has to physically visit every client. Therefore, such meter data can only be at an interval of several months. With the development of smart meters, consumption can be recorded locally up to every minute. Nevertheless, due to transmission and processing costs and/or privacy issues, utility companies often have limited access to these data, for example daily aggregates of individual consumption.

In this thesis, we investigate statistical methods to produce estimations and predictions of spatio-temporally fine-grained electricity consumption. The objective is to use any information that is available, in order to produce high-quality estimation and prediction. Apart from the two sources of electricity consumption data considered above, the methods developed in this thesis use side information to enrich the estimations. Such side information could be global, such as the autocorrelation of the time series of electricity consumption. Or they could be local: exogenous variables corresponding to the periods and individuals in the dataset, for example, the temperature of a given period, or the socio-demographic information, or geographic location of the clients.

Given the importance of load information for the electric network, a plethora of authors have studied load modeling and forecasting. Traditionally, spatial load forecasting [31], focused on forecasting long-term trend in load over spatial regions, and mid-to-short term forecasting [32], focused on the national level, are two separate domains with differnt methodologies. This has changed in recent years, since the demands for load forecasting have changed due to the evolution of the electricity market as well as the availability of large-scale datasets.

On one hand, the spatial scale of short and mid-term load forecasting has drastly decreased [33][34][35][START_REF] Goude | Local Short and Middle Term Electricity Load Forecasting with Semi-Parametric Additive Models[END_REF]. A by-product of this change is that more attention is paid to seeing electric load as multivariate time series which needs to be processed as a whole. In this area, particular problems to be addressed are the similarity between load series [1,35], and the coherent structure between series from multiple components connected through the network [START_REF] Souhaib | Coherent Probabilistic Forecasts for Hierarchical Time Series[END_REF].

On the other hand, this change is also accompanied by the development of new forecasting methods, both from statistics and the machine learning community. As a testimony to this trend, winners for recent load forecasting competitions often used a variety of differnt methods. In Global Energy Forecasting Competition 2014 [START_REF] Gefcom | [END_REF], the winner of the load forecasting track used a quantile regression version of generalized additive model (GAM) to tackle probabilistic load forecasting [START_REF] Gaillard | Additive Models and Robust Aggregation for GEFCom2014 Probabilistic Electric Load and Electricity Price Forecasting[END_REF]. In a 2017 multivariate load forecasting contest organized by Réseau de Transport d'Electricité, the main TSO in France [START_REF] Datascience | Net -Prévision de La Consommation Électrique -RTE[END_REF], the winner used an ensemble of Xgboost 2.2. PROBLEM FORMALIZATION models to achieve an adaptive model with high performance. 2We believe the methods presented in this thesis are a useful complement to these recent developments, as they allow heterogenous data to be preprocessed, before state-of-art load modeling methods can be applied.

Problem formalization

In this thesis, we are interested in the value of a variable v ∈ R (typically electricity consumption), over n 1 consecutive periods, for a group of n 2 individuals. For an individual j, with j ∈ {1, 2, ..., n 2 }, v i,j represents its consumption at period i, with i ∈ {1, 2, ...n 1 }3 . We will use lower-case boldface letters to denote vectors and upper-case boldface letters to denote matrices, so that V ∈ R n 1 ×n 2 is the matrix of the variable of interest, and (v i ) T is the i-th row vector of this matrix, representing the i-th period, and v j is the j-th column vector, representing the j-th individual.

With this notation, we can easily represent both the data and the estimation target of this thesis. In Figure 2.1, a scheme of the variable of interest (blue table) is given. We can also represent the two types of data sources previously mentioned: temporally aggregated measure from a user meter (orange box), and measures taken on from a point on the electricity network, thus an aggregation of several individuals (green box).

The exogenous variables can be represented by matrices that can be put in parallel to the variables of interest (green tables in Figure 2.1).

The estimation targets can also be represented in this way. Most of the time, we aim at estimating the electricity consumption of all individuals for all periods. This is simply every entry of the matrix V. Depending on the applications, one can also be interested in some aggregations of the matrix V. These estimation tasks can be achieved by performing the aggregation on the estimated matrix V. For example, in a change-of-support problem, we would want to estimate the variable of interest on a spatial aggregation grid, that is different from that of the data. One can also be interested in the maximal value of electricity consumption throughout a day, which is the maximum of a number of periods for each user.

As for prediction, it corresponds to estimating the value of the variable for new rows (new periods) and new columns (new users) of V. When no information other than measures on the data matrix is available, it can be difficult to produce predictions. In this thesis, prediction is often achieved by considering • either implicit information about V, such as the temporal auto-regressive structure of the rows of V, or the spatial auto-regressive structure of columns of V;

• or exogenous variables per se, such as the temperature, or socio-demographic information of clients.

This thesis is organized as following: In Chapters 3-4, we present two methods that use individually aggregated data to estimate and predict electricity consumption for aggregations that is different from that of the data. In Chapters 5-9, we present nonnegative matrix factorization, a class of methods that can solve the task in a more general way. The three approaches are compared in Chapter 8, and perspectives are discussed in Chapter 10.

1: The matrix representation of the variable of interest. Typically, an entry of the matrix is the electricity consumption of an individual (a column), during a period (a row). Measures taken by sensors on the network are often aggregations of individuals (sum of the entries in the green box). Measures read on user meters are often temporally aggregated (sum of the entries in the orange box).

Classical and novel solutions 2.3.1 Spatio-temporal kriging

When geographic positions of the individuals are available, one natural way to model the data is to use spatial statistics. By considering the individual values of the dataset as observations of a spatial or spatio-temporal random process, this class of methods provide a toolbox for modeling a wide range of social and natural phenomena. In Chapter 3, we will briefly present the basics of spatial statistics. In particular, we will present kriging, a method to interpolate and extrapolate the value of variables of interest to new positions. Modeling the correlation structure of the spatial random process, this method can produce estimations that are optimal in a sense. By applying kriging to two datasets, we will examine the capacity of this model for solving the estimation problems of Figure 2.1.

Prediction of multiple electricity consumption series

In an approach parallel to spatial statistics, we will present a method to estimate the mean of electricity consumption of unobserved zones, in Chapter 4. This generic approach can estimate electric consumption from source zones, where fine-grained consumption data is available, to target zones where exogenous socio-demographic information is available. By using additional socio-demographic variables available both at source and target zones, this method creates a link between the two by clustering. The electricity consumption is then estimated by regression models. This method is validated on synthetic and real electricity consumption datasets.

CLASSICAL AND NOVEL SOLUTIONS

Matrix-factorization based methods

As suggested by Figure 2.1, the object of interest in this thesis can easily be represented in a matrix form. In Chapters 5 to 9, we propose a general framework based on matrix factorization, which solves the estimation problem of Figure 2.1.

To do this, we first extend nonnegative matrix factorization (NMF) to the matrix recovery context where the observations are general linear measurements of the matrix, instead of matrix entries. In Chapter 5, we describe this extension, and discuss several applications that can be fit into this framework. In addition, we also discuss the recent literature on the global convergence of matrix factorization in the nonnegative case.

In Chapters 6 and 7, we propose two ways to include additional information in NMF. In Chapter 6, we pay especially attention to temporal aggregate data. We propose an efficient algorithm for this data representation, which is guaranteed to converge to a stationary point. We also provide another version of this algorithm, which takes into account the autocorrelation of the consumption of each individual. The utility of these two algorithms are illustrated on synthetic and real electricity consumption datasets, in comparison to reference methods.

In Chapter 7, we propose an extension that takes into account side information: additional time and individual-dependent features. By adding a regression layer into the NMF framework, the algorithm proposed in this chapter can be applied to many application problems. In particular, it can be used to produce electricity consumption predictions for new individuals, from only partially observed consumption data. The experimental study on several datasets shows that this method has very good empirical performance, compared to reference methods.

Part I

State-of-art methods

Introduction to spatial statistics

In spatial statistics, we study the properties of a social or natural phenomenon indexed by its location. In formal terms, the subject is a possibly multivariate random process {Z(s)} with its index in a domain of interest D:

{Z(s) : s ∈ D},
where D is often a subset of R 2 or R 3 .

In general, both the position s and the variable of interest Z(s) can be random. Here is of primary interest the case where the locations are fixed. According to the categorization adopted CHAPTER 3. SPATIO-TEMPORAL CHANGE OF SUPPORT by [START_REF] Cressie | Statistics for Spatial Data: Wiley Series in Probability and Statistics[END_REF], with a fixed set of positions {s}, the data Z(s) emitted by a stochastic process is called geostatistical data.

When the value of Z on a partition of the domain D is of interest, the aggregation or average of Z on small areas is often studied. In this case, the stochastic process is studied on a partition of D, namely a collection of pairwise disjoint subsets {D i } i=1,2,...,n . This is called areal data analysis.

In this chapter, a particular class of questions addressed in spatial statistics, change-ofsupport problem (COSP), is considered. Areal datasets from different sources are often defined on overlapping spatial partitions. Therefore, in order to merge these datasets and study the mutual influence of different phenomena, it is necessary to transfer a dataset from one spatial support to another.

In this section, the basic assumptions of spatial statistics is presented. This section is brief, mostly following [START_REF] Cressie | Statistics for Spatial Data: Wiley Series in Probability and Statistics[END_REF], where a thorough presentation of spatial statistics can be found.

In Section 3.2, methods of univariate COSP are presented. The most important of these is a method called kriging, derived from the geostatistics literature. A few areal interpolation methods other than kriging are presented in Section 3.2.3. In Section 3.3 and Section 3.4, these methods are generalized in the multivariate and spatio-temporal context. Finally, Section 3.5 shows a small empirical study of the methods presented.

Basics in geostatistics

As mentioned above, geostatistics is focused on a random process {Z(s)} indexed by its fixed location s. The domain of interest D, of permissible values for s, is a subset of R d (in practice, d = 2 or 3). A realization of this random process is the observed value of Z, z 1 , z 2 , ..., z n , at a number of locations: s 1 , s 2 , ..., s n .

The random process model works well for applications in geology, mining, and environmental science, where a natural phenomenon, for example the ore reserve or air pollution, is measured by sensors at several locations, and something is to be said about the whole area.

The random distribution of a stochastic process is defined by a probability measure on the function space {f : D → R}, which can be very general. We will focus on Gaussian processes, namely, for any number of positions {s 1 , ..., s n }, the joint distribution of {Z(s 1 ), ..., Z(s n )} is a multivariate Gaussian. A Gaussian process is completely characterized by its first two moments:

• a mean function m : D → R, such that m(s) = E(Z(s)), ∀s ∈ D,
• and a covariance function

C : D×D → R, such that c(s 1 , s 2 ) = cov(Z(s 1 ), Z(s 2 )), ∀s 1 , s 2 ∈ D.
Similar to time series analysis, which is a statistical analysis based on a stochastic process index by an interval of R, reasonable assumptions of stationarity have to be made while analyzing geostatistical data. Several kinds of stationarity is defined in spatial process.

One useful assumption to consider is second-order stationarity. {Z(s)} is second-order stationary if

• E(Z(s)) = µ, which is a constant for all s ∈ D;

• there is a covariance function C(•) such that cov(Z(s 1 ), Z(s 2 )) = C(s 1 -s 2 ).

Additionally, if C(s 1 -s 2 ) is in fact only a function of s 1 -s 2 , then this process is isotropic.

INTRODUCTION TO SPATIAL STATISTICS

A less strict assumption is intrinsic stationarity. {Z(s)} is intrinsically stationary if • E(Z(s)) = µ, which is a constant for all s ∈ D;

• there is a variogram function γ(•) such that Var(Z(s 1 ) -Z(s 2 )) = 2γ(s 1 -s 2 ).

Additionally, if 2γ(s 1 -s 2 ) is in fact only a function of s 1 -s 2 , then this process is isotropic.

All second-order stationary process is intrinsically stationary because if the process has a covariance function C(•), then

Var(Z(s 1 ) -Z(s 2 )) = E[(Z(s 1 ) -Z(s 2 )) 2 ] = E(Z(s 1 ) 2 ) + E(Z(s 2 ) 2 ) -2E(Z(s 1 )Z(s 2 )) = C(0) + µ 2 + C(0) + µ 2 -2(C(s 1 -s 2 ) + µ 2 ) = 2C(0) -2C(s 1 -s 2 ) ≡ 2γ(s 1 -s 2 ).
Conversely, not all intrinsically stationary processes are second-order stationary. Consider the Lévy Brownian motion with d parameters {B(s)}. The variogram is well defined: Var(B(s 1 ) - 41, p. 68]). In univariate kriging, we mostly consider intrinsically stationary processes, and thus studying variograms and not covariances, while in multivariate and spatio-temporal context, covariance function is preferred.

B(s 2 )) = s 1 -s 2 . However, cov(B(s 1 ), B(s 2 )) = s 1 + s 2 -s 1 -s 2 , which is not a function of s 1 -s 2 ([

Variogram

For a function C or γ defined on D with values in R to be a covariance function or variogram, it must satisfy that for any n points in D: s 1 , ..., s n , the variance-covariance matrix (cov(s i , s j )) 1≤i,j≤n is positive semi-definite. For stationary processes, this is equivalent to conditional negative-definiteness on variograms:

∀a 1 , ..., a n ∈ R not identically zero, such that n i=1 a i = 0, n i=1 n j=1 a i a j γ(s i -s j ) ≤ 0.
We should notice that a function that is only semi-definite positive but not definite positive can still be a valid covariance function for a stationary Gaussian process. This means that, in general, the variance-covariance matrix (cov(s i , s j )) 1≤i,j≤n for n distinct points s 1 , ..., s n can be singular. However, it would be convenient (for kriging, for example) if we could restrict covariance matrix to be non-degenerate. Fortunately, this is the case for all covariance functions considered below ([42, p. 64] shows a sufficient condition for this to be satisfied).

Covariance functions or variograms are important in the specification of a spatial model, because their regularity implies the regularity of the random process (see [START_REF] Robert | The geometry of random fields[END_REF] for more thorough treatment). For instance, a stochastic process is continuous in mean square at a position s , if and only if c(s, s ) is continuous at s = s = s . As a reminder, mean-square continuity means that E(|Z(s n ) -Z(s )| 2 ) → 0 for all sequence (s n ) n convergent to s . For a stationary process, this is equivalent to the covariance function C being continuous at 0. In a similar fashion, the k-th order partial derivative ∂ k Z(s)/∂s i 1 ...∂s i k exists as a mean square limit, if and only if

∂ 2k C(s)/∂ 2 s i 1 ...∂ 2 s i k exists at s = 0.
Several classes of covariance functions are extensively studied and commonly used in the Gaussian process and spatial statistics literature. When studying an isotropic covariance function, two parameters are often considered • the sill, which is C(0), the value of the covariance function at 0;

• and a scale parameter, sometimes called the range, which is the minimal distance at which we can consider that the covariance is 0.

For a semivariogram, the sill is the maximal semivariance, and the range is the distance at which the semivariance stabilizes to the sill.

Below, a number of common isotropic stationary variograms are presented (see, for example, [44,Chp. 5] for other variograms). Note that all of these variograms are isotropic regardless of the dimension of the index space, which is not true for general covariance functions.

• Exponential variogram:

γ(h) = c(1 -exp(-h a )
), where c is the sill, a is a scale parameter, controlling how quickly γ(h) approaches the sill. The range of exponential variogram is +∞ in theory. Since its derivative at 0 is not zero, this variogram is appropriate for processes with high variation at a small scale.

• Gaussian variogram: γ(h) = c(1 -exp(-h 2 a 2 ))
, where c and a are the sill and the scale parameter. The major difference with exponential variogram is that, seen as a function defined on the real line, Gaussian variogram is infinitely differentiable at 0. Therefore, the associated process is also infinitely differentiable in mean square.

• Matérn variogram:

γ(h) = c(1 -2 1-2ν Γ(ν) ( √ 2νh
a ) ν K ν ( √ 2νh
a )) (see [START_REF] Michael L Stein | Interpolation of spatial data: some theory for kriging[END_REF] for more detail). Additionally to the sill and the scale parameter, Matérn class of variograms has a regularity parameter ν (Γ(•) and K ν (•) are the gamma function and the modified Bessel function of order ν). An interesting property of Matérn variogram is that the associated process is k-times mean square differentiable for all k < ν. Therefore, by setting the value of ν, the regularity of the stochastic process is controlled. Another feature is that when ν = p + 1 2 where p is an integer, the general formula becomes the product of an exponential function with a p-order polynomial.

• Nugget effect: γ(h) = b if h = 0 and γ(0) = 0. The nugget effect is often added to other variograms. If it is used alone, the random Gaussian process is a white noise, with no spatial correlation. The nugget effect is introduced into spatial modeling because empirically, an estimated variogram rarely has 0 value at h = 0. It is often interpreted as resulting from a measurement error on the data.

The discussion above established that variogram is a useful tool in geostatistics. However, in most empirical problems, the variogram is unknown. It has to be estimated from the data, before being used in an application.

Without imposing a functional form for the covariance, one can begin by estimating an empirical variogram. The semivariance at different distances is estimated on the data, while the analytical formula of the variogram is not explicit. In the pioneer work of Matheron [START_REF] Matheron | Traité de géostatistique appliquée[END_REF], the following estimator is proposed:

2γ(h) = 1 |N (h)| N (h) (Z(s i ) -Z(s j )) 2 , N (h) = {(i, j)|s i -s j = h}.
For irregularly spaced data, we can either be interested in the scatter-plot of semivariance with respect to distance, or group pairs of observations into bins of distance.

To conduct a spatial statistics study, the empirical variogram is estimated and visually examined first. This helps to decide on the family of covariance function, and whether or not to 3.2. KRIGING IN UNIVARIATE CHANGE OF SUPPORT PROBLEM adjust the anisotropy. In a second place, a chosen parametric covariance function model fitted on the data. This establishes an estimation of covariance between all potential points and can be used later in kriging.

Kriging in univariate change of support problem

Kriging

An important task in spatial statistics is the prediction of the variable of interest at locations not present in the data. Under different contexts, the interpolation of spatially-index data is called kriging, Gaussian process regression, or Wiener-Kolmogorov prediction. This and the next section present the basic techniques of kriging. More detailed description of these methods can be found in many reference books or lecture notes ( [START_REF] Cressie | Statistics for Spatial Data: Wiley Series in Probability and Statistics[END_REF], [START_REF] Bohling | Data Analysis in Engineering and Natural Science[END_REF]).

As in time series, the kriging predictor falls under the category of best linear unbiased predictor (BLUP). Namely, it is a weighted average of the data, with expected value equal to that of the process at this location, of minimal mean squared residuals.

Consider a Gaussian process {Z(s) : s ∈ D}, with mean function m and stationary covariance function C. This actually means that the residual component of the process, Z(s)-m(s) ≡ R(s) is a second-order stationary process of mean 0 and of covariance function C. Suppose that this process is observed at n different locations {s 1 , ..., s n }. An unbiased linear predictor Ẑ of Z at location s has the form:

Ẑ(s) -m(s) = n i=1 λ i (Z(s i ) -m(s i )),
where λ i is the weight of the i-th observation. Define the residual of the estimator by R(s) ≡ Ẑ(s) -m(s). The weights are then to be determined in order to minimize the quadratic risk,

E[( Ẑ(s) -Z(s)) 2 ] = Var( Ẑ(s) -Z(s)) = Var( R(s) -R(s)).
The three kriging methods presented in this section differ in their treatment of the mean function m(•). An empirical comparison of these methods on the log-density of population can be found in Section 3.5.

Simple kriging

The simplest case of kriging is called simple kriging. In this scenario, the intrinsically stationary process {Z(s) : s ∈ D} considered has a constant known expected value, m(s) = µ, for all s ∈ D. The quadratic risk can be written as:

Var( R(s) -R(s)) = Var( n i=1 λ i R(s i ) -R(s)) = Var( n i=1 λ i R(s i )) + Var(R(s)) -2cov( n i=1 λ i R(s i ), R(s)) = n i=1 n j=1 λ i λ j C(s i -s j ) + C(0) -2 n i=1 λ i C(s i -s).
The minimization of the quadratic risk is an unconstrained quadratic optimization that can be solved analytically. Note by λ = (λ 1 , ..., λ n ) T the column vector of weights, K = (C(s is j )) 1≤i,j≤n the matrix of mutual covariance between locations, and k = (C(s 1 -s), ..., C(s n -s)) T the column vector of covariance between s and the observations. Because of the positive-definite CHAPTER 3. SPATIO-TEMPORAL CHANGE OF SUPPORT property of C, the matrix K is also positive definite, and therefore invertible. Then it becomes clear that minimizing the quadratic risk is equivalent to minimizing

λ T Kλ -2λ T k = (λ -K -1 k) T K(λ -K -1 k) -k T K -1 k.
The optimal weights is therefore λ = K -1 k.

It can be noticed that the optimal weights do not depend on the value of observations Z(s i ), only the covariances.

In practice, if a large number of observations are available, the inversion of matrix K can be cumbersome. There might also be numerical problems with covariances calculated on an estimated variogram. In that case, a form of local kriging can be considered, where only the closest neighboring observations are used to provide prediction of a given location. The procedure remains essentially the same.

Ordinary kriging

In the ordinary kriging case, the mean function of the Gaussian process, m(s) is unknown. An additional assumption of constant local mean is added: for all neighboring locations involved in kriging, the expected value is unknown but constant. Namely, for all 1 ≤ i ≤ n, m(s i ) = µ = m(s), where µ is an unknown parameter. One way to specify a linear predictor

Ẑ(s) -µ = n i=1 λ i (Z(s i ) -µ)
without actually knowing µ is to impose that n i=1 λ i = 1, so that Ẑ(s) = n i=1 λ i Z(s i ). In this case, the quadratic risk is minimized under the linear equality constraint 1 T λ = 1. The minimizer of this problem is the solution of the linear system

K 1 1 T 0 λ α = k 1 ,
where α is a Lagrange multiplier to impose the constraint.

Universal kriging Now generalize the model to include a p-dimensional explanatory variable x(•) available everywhere in D, especially at observation locations and the kriging target. The mean function has the form

m(s) = x(s) T β, ∀s ∈ D,
where in x(•) there is potentially a constant component, and β is an unknown regression coefficient vector.

Similarly to ordinary kriging, we need to specify

Ẑ(s) -x(s) T β = n i=1 λ i (Z(s i ) -x(s i ) T β),
without knowing β. This is equivalent to specify Ẑ(s) = n i=1 λ i Z(s i ) under the constraint

x(s) = n i=1 λ i x(s i ) ≡ Xλ,
where X is a p-by-n matrix combining explanatory variable at all observations.

KRIGING IN UNIVARIATE CHANGE OF SUPPORT PROBLEM

This constraint leads to a p-dimensional Lagrange multiplier α. The minimizer of quadratic risk under the constraint is the solution of the linear system

K X T X 0 λ α = k x(s)
.

Variance of kriging predictors Because the processes are specified as Gaussian, the variance (which is equal to the quadratic risk in the case of unbiased predictor) of a linear interpolation method can be calculated with ease. In the case of kriging methods presented above, it only depends on the position of the observations and predictions, and not on the value. This provides a measure of uncertainty in the predictions obtained by kriging methods.

Block kriging

Kriging methods are the most elementary way to extrapolate observations to unobserved positions. However, to solve the change of support problem, one often needs to deal with average value of a process on areal support. To deal with this scenario, block kriging can be used.

Consider n disjoint areas in the domain, B 1 , ..., B 1 ⊂ D. The observations are the average of the process Z(•) on these are:

Z(B i ) = s∈B i Z(s)ds, ∀1 ≤ i ≤ n.
These integrals are defined with minimal regularity (for example Z(•) being continue on path and B i compact. This is the case for all covariance functions considered in this document). The target of kriging, Z(B), is the average on a new area B.

Because of the linearity of the averaging operation, the calculation in the previous section applies naturally on the case of block kriging, with minimal assumptions of regularity. Instead of covariance between points, integrals of the covariance function on areas are now used:

C(B i , B j ) = B i B j C(s i -s j )ds i ds j .
In universal block kriging (the most general case), the explanatory variables also need to be integrated. In that case, the optimal linear weights λ are the solution of the linear system

K X T X 0 λ α = k x(B)
,

where the elements of K, X, and k are all integrated.

In order to solve this kriging equation, it is necessary to insure that the block covariance matrix, K, is invertible. This can be difficult for general covariance models.

In practice, integrals of the covariance function can be approximated numerically by dividing the areas on a regular grid. To limit the time spent on these integrals, it is often interesting to use local kriging.

Bayesian estimation and prediction

Although the geostatistical methods are presented here in a frequentist point of view, Bayesian estimation is actually rather prevalent in the spatial statistics literature. Because of the hierarchical nature of spatial and spatio-temporal models, Bayesian estimation methods such Markov chain Monte Carlo (MCMC) algorithms are very popular in applied studies (see [START_REF] Cressie | Statistics for spatio-temporal data[END_REF], [START_REF] Banerjee | Hierarchical modeling and analysis for spatial data[END_REF] for a detailed description of these methods and numerous examples). However, such methods can be intricate to implement, and by definition MCMC is often expensive in computation time. This makes MCMC methods difficult to scale.

Areal interpolation beyond kriging

Areal interpolation Historically, the area-to-area change of support problem is solved by areal interpolation [START_REF] Michael F Goodchild | Dept. of Geography. Areal interpolation: a variant of the traditional spatial problem[END_REF], without explicit generative modeling assumptions. In this procedure, the interpolation weights are obtained proportionally to the area of the blocks considered.

This line of methods is further extended [START_REF] Flowerdew | Statistical methods for inference between incompatible zonal systems[END_REF] into areal regression models, where the relationship between a co-variable (known everywhere) and the variable of interest is estimated in a regression model. This relation is then used to produce estimation of the variable on an overlapped area. Areal regression models are generalized in a Bayesian context by [START_REF] Andrew S Mugglin | Fully model-based approaches for spatially misaligned data[END_REF]. For more historic context in general change of support problems, see [START_REF] Carol | Combining incompatible spatial data[END_REF]. Instead of deriving a block correlation matrix by integrating on an underlying Gaussian process, as in block kriging, in this approach, the relationship between area is derived directly from the spatial relationships of the blocks. The blocks are viewed as vertices of a graph, whose edges can be defined on either contiguity or by setting a distance threshold. The edges can also be weighted according to relevant criteria of the variable of interest [START_REF] Roger S Bivand | Applied spatial data analysis with R[END_REF].

Spatial statistical model on areal data

Once the covariance structure of blocks is specified, the joint distribution of Z(D 1 ), ..., Z(D n ) can be estimated from the data. Often, the Markov property is assumed for the process, so that the theory of Markov random field can be used in the model specification and estimation (see [START_REF] Cressie | Statistics for spatio-temporal data[END_REF] and references within for a more detailed exposition). Moreover, it is often possible to infer the distribution of the value on a target block, given the estimated model. This provides an interesting way to solve area-to-area change of support problem, especially when combined with Monte Carlo simulations [START_REF] Chiles | Geostatistics: modeling spatial uncertainty[END_REF]Chap. 7]).

Multi-scale models Additional structure can be added into a spatial model on areal data to fit multi-scale data. Tree models are particular adapted for integrating data from several resolutions ( [START_REF] Huang | Fast, resolution-consistent spatial prediction of global processes from satellite data[END_REF], [START_REF] Eric | Multiscale statistical models for hierarchical spatial aggregation[END_REF]). The connection between different grid levels is introduced to link multiple nested spatial areal models, under the same graphical model representation. Because of the complexity of the model, this line of research mostly uses Bayesian methods for estimation and prediction.

Multivariate kriging

In universal kriging, exogenous variables are used to explain the mean of the interpolated variable. One way to generalize this model is to include the spatial covariance between different variables. This technique, called cokriging, has a wider range of application as we will see in this section.

MULTIVARIATE KRIGING

Cross-covariance

Formally, cokriging studies a Gaussian process {Z(s) : s ∈ D} indexed s in the domain D, and with values in R d , where d is the number of variables in question. Similarly as in Section 3.1, we define a multivariate Gaussian process as second-order stationary if it has constant mean and a stationary cross-covariance function:

• E(Z(s)) = µ, for all s ∈ D; • there is a cross-covariance function C(•) such that cov(Z(s 1 ), Z(s 2 )) = C(s 1 -s 2 ), for all s 1 , s 2 ∈ D.
In this definition µ is a d-dimensional vector, and

C(s 1 -s 2 ) is a d-by-d matrix. Note that the covariance function C(•) is not restricted to be symmetric: in general, cov(Z i (s + h), Z j (s)) = cov(Z j (s + h), Z i (s)), for 1 ≤ i < j ≤ d.
Similar to the univariate case, not all matrix-valued function can be a valid cross-covariance. The following properties need to be verified:

• the limit of C(h) at h = 0 has to be positive definite;

• cov(Z i (s + h), Z j (s)) = cov(Z j (s), Z i (s + h)), in other words C(h) = C T (-h); • generally, for any n points in D, s 1 , ..., s n , ∀a 1 , ..., a n ∈ R d , n i=1 n j=1 a T i C(s i -s j )a j > 0.
Cross-covariance functions are a more complicated object than covariance. In particular, it is difficult to build matrix-valued functions satisfying the three properties mentioned above. Two methods are presented below. The rest of approaches are readily available in references (for example in [START_REF] Banerjee | Hierarchical modeling and analysis for spatial data[END_REF]Chap 9]).

Separable model C(h) = ρ(h)T

, where ρ(•) is a univariate correlation function (ρ(0) = 1), and T is a symmetric positive-definite matrix. This model is easy to understand and to estimate, but it can be too strong in practice: the cross-covariance is necessarily symmetric, and all variables have a single spatial-correlation range (that of ρ).

Linear model of coregionalization

We make the assumption that Z(s) = Aw(s), with A a d-by-d full rank matrix, and w(s) a d-dimensional Gaussian process with independent components and unit variance. The cross-covariance of Z(•) is then

C(h) = d i=1 ρ i (h)A i A T i ,
with ρ i (•) the correlation functions of w i . The cross-covariance matrix is still symmetric, but each component in Z(s) has a distinct covariance function.

The estimation of linear coregionalization model can be achieved either by an iterative least square method (alternatively vary the parameters in individual correlation functions, and in the matrix A, see [START_REF] Goulard | Linear coregionalization model: tools for estimation and choice of cross-variogram matrix[END_REF]), or integrated into a global Bayesian context [START_REF] Finley | spBayes: An R Package for Univariate and Multivariate Hierarchical Point-Referenced Spatial Models[END_REF].

In the multivariate case, intrinsically stationary processes are less considered. In fact, the natural generalization of variogram is the cross-variogram Var(Z(s 1 ) -Z(s 2 )). This matrix is CHAPTER 3. SPATIO-TEMPORAL CHANGE OF SUPPORT always symmetric whether the cross-covariance is symmetric. In order to have an object that contains as much information as the cross-covariance, a pseudo cross-variogram is defined as (Var(Z j (s 1 ) -Z i (s 2 ))) 1≤i,j≤d [START_REF] Myers | Pseudo-cross variograms, positive-definiteness, and cokriging[END_REF]. This quantity arises naturally in the variance of the cokriging predictor as the variogram in univariate kriging.

Cokriging

Once specification and estimation of cross-covariance is pinpointed, the cokriging equation itself can be derived in a similar fashion as in universal kriging.

Under constant mean, a linear predictor of one component Z 1 of Z at s, with data available at s 1 , ..., s n , is of the form

Ẑ1 (s) -µ 1 = n i=1 d j=1 λ ij (Z j (s i ) -µ j ).
In order to eliminate the unknown means and guarantee unbiasedness, the linear constraints,

n i=1 λ i1 = 1, n i=1 λ ij = 0, ∀2 ≤ j ≤ d, are imposed.
The quadratic risk minimization leads to a similar matrix equation as in ordinary and universal kriging with bigger matrices. Note by K the nd-by-nd cross-covariance matrix,

K =       K 11 K 21 • • • K d1 K 12 K 22 • • • K d2 . . . . . . . . . . . . K 1d K 2d • • • K dd      
, where K ij is a n-by-n cross-covariance matrix between the i-th and j-th component of the process, on the n positions of the data. The kriging coefficients and covariance between target position and available positions are stacked into nd-dimensional vectors in a corresponding fashion, namely,

λ =              λ 11 λ 21 . . . λ n1 λ 21 . . . λ nd              , k =              cov(Z 1 (s), Z 1 (s 1 )) cov(Z 1 (s), Z 1 (s 2 )) . . . cov(Z 1 (s), Z 1 (s n )) cov(Z 1 (s), Z 2 (s 1 )) . . . cov(Z 1 (s), Z d (s n ))             
.

By defining a d-dimensional Lagrange multiplier α, the kriging equation is

K I d ⊗ 1 I d ⊗ 1 T 0 λ α = k e 1 ,
where I d ⊗ 1 is the Kronecker product between the d-by-d identity matrix and a n-dimensional unit vector, and e 1 is the first vector of canonical basis in R d .

CHANGE OF SUPPORT IN SPATIO-TEMPORAL MODEL

Change of support in spatio-temporal model

Spatio-temporal covariance and kriging

As in the spatial case, the key to spatio-temporal kriging is the specification of covariance function model. [START_REF] Cressie | Statistics for spatio-temporal data[END_REF]Chap 6] is a recent review of methods in this vein. Here we will focus on the two principal ways to integrate the temporal component into the spatial Gaussian process.

Continuous-time spatio-temporal model

Continuously indexed time can be viewed as an additional dimension in the index of the Gaussian process. For an interval I of R, the spatiotemporal process is

{Z(s, t) : s ∈ D, t ∈ I}.
A second-order stationary spatio-temporal process has a constant mean, and a stationary covariance function

C(h, τ ) = cov(Z(s, t), Z(s + h, t + τ )), for all s ∈ D, t ∈ I.
As in the purely spatial case, the covariance function has to be positive-definite. The most direct way to achieve this is to treat time as an additional dimension of the index space, and use one of the spatial covariances mentioned in Section 3.1. The underlying assumption would be that time is a dilated dimension, which is often not the case.

A separable spatio-temporal covariance is one that verifies C(h, τ ) = C 1 (h)C 2 (τ ). This can be insufficient. Typically, if we want to obtain the value of the variable of interest at a moment t at a location s, and if this variable is observed at every other location except s, observations at other moments can simply be tossed away, since just doing spatial kriging already gives the optimal estimation. This condition is often too restrictive, for example, to model the propagation of a wave across the spatial domain.

To generalize a separable spatio-temporal covariance, we can notice that the sum or product of valid covariance functions are still valid covariance function ( [START_REF] Cressie | Statistics for spatio-temporal data[END_REF]Theorem 6.1]). This procedure allows us to produce non-separable but fully symmetric covariances easily. These covariance functions are called product-sum models. They are usually defined as

C(h, τ ) = C 1 (h) + C 2 (τ ) + kC 1 (h)C 2 (τ ), with a joint parameter k, given spatial covariance C 1 (h) and temporal covariance C 2 (τ ).
The most general way to specify spatio-temporal covariances is to start from a positive spectral measure. As seen in Section 3.1, there is a correspondence between positive-definite functions and positive measures in the frequency space. Therefore, covariance functions can be specified as inverse Fourier transform of positive measures. Moreover, [START_REF] Cressie | Classes of nonseparable, spatio-temporal stationary covariance functions[END_REF] showed that for a function α

: (ω, τ ) → α(ω, τ ) defined on R d × R, that is positive-definite in τ , the function C(h, τ ) = e i<h,ω> α(ω, τ )dω, h ∈ R d , τ ∈ R, is a positive-definite function in R d × R.
This is an easy way to construct very general covariance functions.

Fitting spatio-temporal functions to data can be done in a similar way as in cokriging, although given the complexity of the model, Bayesian methods are preferred in this context.

Discrete-time spatio-temporal model

Time can also be seen as discrete. For each position s i in the data, a discrete time series {Z t (s i ), t = 1, 2, ..., T } is observed. Often, such situation arises when the data is on a regular temporal rhythm.

Without assumption on temporal stationarity, the covariance structure can be treated in the same way as cross-covariance in cokriging. If temporal stationarity is assumed, then lagged auto-covariance can be estimated on averaging pairs of observations of the same lag.

More sophisticated dynamic models can be considered with an appropriate physical phenomenon, but that would be beyond the scope of the present document (see[ Chap 7] [48] and references there within for details).

Computational issues in spatio-temporal kriging

Once the covariance function is estimated and applied on the data, the kriging equations is exactly the same as in cokriging. The most important issue in this approach is scalability. Indeed, solving for the kriging weights needs the inversion of the matrix K. In a multivariate dataset, K is of dimension nd-by-nd, whereas in a spatio-temporal dataset, it is nT -by-nT .

Application

In order to test the geostatistical methods presented before, two illustrative examples (one purely spatial kriging, and one of spatio-temporal kriging) are presented next. The kriging methods used are implemented in the gstat package in R. The research article [START_REF] Edzer | Multivariable geostatistics in S: the gstat package[END_REF] is a presentation of this package. [START_REF] Bivand | CRAN Task View: Analysis of Spatial Data[END_REF] is a general presentation of spatial data analysis in R.

Variogram of the density of population in Lyon

In the following, the logarithm of population density (population divided by area) of IRIS in Lyon is used as an illustrating example (Figure 3.

1). An IRIS (Ilots Regroupes pour l'Information Statistique) is an administration zone division of approximately 2,000 habitants, used by the French official statistics department (INSEE). The logarithm of population density of an IRIS is calculated as log 10

Population in the IRIS Area in the IRIS (km 2 ) .

The domain of interest is included in an approximately 10km-by-10km square. Eache IRIS is colored with a color scale that correponds to its population density. As we can see, this data has some obvious discontinuity on the frontier of the IRISes. In the rest of this section, we suppose this density is attached to the geometric center of each IRIS, and try to find an estimation for each point in the area.

Suppose that the model is isotropic (the case of anisotropy is often resolved by applying a re-scale on some axes). Figures 3.2 and 3.3 shows the boxplots of semivariance of distance smaller than 5000 meters, grouped in 15 bins, and the empirical variogram, calculated for the 15 bins.

As can be seen in Figures 3.2 and 3.3, an isotropic variogram is often an increasing function on the positive real line. Actually, γ(•) is a symmetric function with respect to h = 0 defined on the real line. Therefore, when γ is differentiable at 0, its derivative at 0 has to be 0 too. It can also be noticed that the semivariance somewhat stabilizes after a period of rapid increase. This maximal semivariance is the sill in the geostatics literature. The distance h at which the semivariance stabilizes is the range. On the corresponding covariance function, this is the minimal distance at which the spatial correlation practically disappears. Once the geostatistician made a decision on which variogram to use, a fitting method is used to estimate the parameters in the variogram model. Several methods are studied and implemented in statistical softwares (for example gstats). These methods are mostly maximum likelihood methods, based on the Gaussian assumption of the data, or some weighted version of least square estimation, by minimizing a squared empirical risk. Figure 3.5 shows these variograms models on the log of population density, estimated by the default fitting method implemented in gstat (least square weighted by the number of observation points divided by squares of the distance). The estimated sill, nugget, and scale parameters are reported in Table 3.1.

Kriging of a univariate variable

A 100m-by-100m grid is created on the original spatial domain. Figure 3.6 shows the log-density of population assigned to the geometric center of each IRIS. The population density of other squares are estimated by universal kriging described in Section 3.2: the mean of the Gaussian process is assumed to be constant but unknown. The kriged value surface is smooth for all three variogram models, much smoother than linear interpolation. Although the exponential covariance is not differentiable at zero (therefore the kriged value surface is not differentiable either), it is only slightly less smooth than Gaussian and Matérn covariance models. At points where observations are available, there is sometimes a discrepancy even with its nearest neighbors, for example the geometric center of the most southern IRIS.

Figure 3.8 shows the kriging variance for the grid. The variance is evidently larger at the edge of the domain of interest: the distance with any observations being larger, little information is provided for these points.

Block kriging

The data used in this example is area-based. The reported log-density of population is the mean of the IRIS. Therefore, the point kriging method used in the previous section is a simplification of the actual data. In this section, block kriging and areal interpolation is compared to the last previous. The kriging target is the collection of blocks obtained by triangulation of the position of local substations.

As presented in Section 3.2, the major difference between block kriging and classical pointto-point kriging, is that the covariance function is evaluated and integrated on the source and/or target blocks. In gstat package, a numerical approximation is introduced to calculate point-toarea covariance, by dividing target blocks on a grid. For the moment, block-to-block covariance is not implemented in the kriging feature of the package.

To be the most accurate, with area-based data, the covariance model should be fitted considering the empirical covariance as block-based as well. Here the point-based fitted covariance models are used.

In Figure 3.9, the kriged mean value of target blocks are shown, so are the IRIS frontier (the source data are assumed to be point-based). Similar to the previous section, areal interpolation produces a value surface that is much less smooth than the kriged value surface. Exponential covariance is the least smooth of the three covariance functions, but the difference is rather small.

In Figure 3.10, the variance of kriged value is shown. Block kriging methods have smaller variance than point kriging, which is a natural, since it is the block averages that is estimated.

Kriging of spatio-temporal data

In this section, the simulated weekly electricity consumption of MV feeders with local substations overcrossing IRIS of Lyon for the 2010-2012 period is studied.

The simulations are produced in the following manner. We have the geographical positions of the low-voltage (LV) distribution transformers connected to the feeders. These transformers are therefore positioned inside the boundaries of each IRIS. For each transformer, we choose a client type according to socio-demographic information of the IRIS: the ratio between the number of apartments and contractual power (CP), and the ratio between the number of warehouses and CP. Using these two ratios, the transformer is assigned a client type whose hourly consumption time series is then simulated with regulated French profiles of typical consumers [64]. We plug real local temperature and calendar variables for the 2010-2012 period into the profiles to simulate hourly consumption shape of each LV transformer, and multiply it by the total contractual power of the transformer. The IRIS aggregation is calculated by summing over transformers insider its boundary, while the feeder aggregation is obtained by summing over transformers connected to it. The weekly consumption is obtained by taking the averages over hourly simulations. This dataset is also used in the next chapter. See Section 4.3.1 for more details on the simulation process.

The data spans 156 weeks (three years), and 114 MV feeders are included. Figure 3.16 shows the first four periods of data.

To understand the temporal profile of the data, the spatial average of the 122 time series is studied. This mean series has a strong periodic pattern (Figure 3.11a). Given the sample autocorrelation and partial autocorrelation (Figure 3.12), an AR(9) model with an annual cycle (Figure 3.11b) is chosen to estimate the spatial averages. A period of 9 weeks is used since it is large enough to cover seasonal dependency. We can see that the annual cycle explains quite well the time series of spatial averages.

We first apply purely spatial block kriging to each of the four periods (with the temporal trend removed), before comparing spatio-temporal kriging to this benchmark. Figure 3.13 shows the empirical variogram and a fitted order-1.5 Matérn variogram model for each of the four weeks. The parameters of the fitted variogram model are reported in Table 3.2. It can be observed that the variograms are different for the four weeks, which suggests a joint spatio-temporal modeling of such a covariance structure can be difficult.

An empirical spatio-temporal variogram is estimated on the data, minus the fitted time series model (the leftmost panel in Figures 3.14,3.15). The semi-variance increases quickly and plateaus out once the spatial distance exceeds a threshold of less than 1000 meters, suggesting a rather small spatial covariance. On the contrary, increase temporal lag does not increase much 3.3). Both model constructions are provided in gstat tool box. The separable model has a joint sill, while the spatial and temporal components have unit sills (as the sum of nugget effect and a Matérn covariance). The parameter k in sum-product model is the weight of the product term. In gstat, parameters are optimized using generic optimization function optim of R. Here the optimization routine is L-BFGS-B, a quasi-Newton optimization method for box domains.

The fitted variograms can be visually compared in Figure 3.14 and 3.15. It can be noticed that it is difficult to fit the decrease in semi-variance around 500m, either by separable or sumproduct model. The sum-product model seems more appropriate than the separable model, in that for large spatial distances, the increase of semi-variance when temporal lag increases is kept in fitted model. Spatio-temporal (point) kriging is applied using the fitted variograms. In Figure 3.17, the spatio-temporally krigged weekly IRIS electricity consumption is represented. The first and the third rows of Figure 3.17 are respectively the separable and product-sum model. Because the kriging computation becomes very expensive for non-separable variograms, only four weeks of data are used in the product-sum model. The second row shows the purely spatial kriged value of the same data. In the fourth row, the actual simulated consumption of the IRIS is plotted (the ground truth).

We observe that compared to the data (Figure 3.16 and the fourth row of Figure 3.17), the kriging values are much more smooth spatially. The temporal variation between the first week and next three weeks is also decreased in the kriging results. The separable covariance does the most spatial and temporal smoothing: the product-sum model seems more smooth spatially on the third date: there is the least contrast between different dates and areas. The product-sum model has the least smooth spatial contrast, while purely spatial kriging is somewhere in the middle.

Error rates of the kriging estimation are shown in Table 3.4. Two error metrics are used here

• the relative root-mean-squared error: RRMSE(x, x * ) = x-x * 2 2 * F ,
• and the mean absolute percentage error: MAPE(x,

x * ) = 1 n n t=1 |xt-x * t | x * t .
where x * is the true vector of dimension n to be estimated, and x is the estimated value. Indeed, the error rates of the kriging methods are all very high. By examining Figure 3.17, we realise that the kriging estimations are much too smooth compared to the true value, hence the high error rates.

In fact, by examining both the data (Figure 3.16) and the instant spatial variogram (Figure 3.13) that the stationary hypothesis, which we adopted from the beginning of this chapter, is not verified by this dataset. This indicates that although we can obtain an estimation by kriging, if we only model the variation by the covariance in a stationary process, the quality of the estimation can not be very convincing.

Conclusions

Open questions From the above examination of the spatio-temporal statistics methods, it is seen that:

• An extensive literature in geostatistics and related domain studies the problem of kriging, and the theoretical foundation of this technique has been greatly enriched in the past 20 years.

• Although there is a developed theory, not all of the methods have been implemented in a satisfactory way to solve every instance of kriging problems. In particular, block kriging in a spatio-temporal context, or block-to-block kriging has not really been implemented. There are two reasons for this lack: the first is the relatively high cost of doing block kriging, the other is the lack of application of such methods in geostatistics.

• The methods are general quite slow, especially when the models are complex and MCMC estimation has to be used for estimation.

CONCLUSIONS

Inspirations drawn from the spatial statistics literature for our problem In the rest of this thesis, the main objective to be achieved is a computationally feasible method for spatiotemporal interpolation estimating the consumption time series at the IRIS level, from data collected at another spatial division, probably of a similar level or below. From this preliminary examination of the literature, spatio-temporal kriging is not well suited for the problem introduced in the introduction, as tempting as this option sounds. Unlike in physical phenomena studied in geostatistics, electricity consumption is a phenomenon is mainly driven by human activities. Thus using appropriate exogenous variables to correctly model the mean function, probably has more impact to the performance of the methodology, than optimizing over the choice of spatial covariance structure. Moreover, the computational complexity of spatio-temporal kriging is generally too high for the applications considered in this thesis. In the previous section, we were only able to produce spatio-temporal kriging during four periods, whereas in the following chapters, experiments are generally conducted for data involving hundreds or thousands of periods.

However, a number of important lessons are to be kept in mind for the rest of the methods developed:

• Given an appropriate modeling of the trend, studying the spatio-temporal covariance structure of the residuals can probably enhance the model, to deliver better performance.

• The exogenous variables used in the trend modeling can be kriged to be estimated at places different than the data source. This is especially true for weather processes which, by their nature, has a spatio-temporal covariance structure.

• Even if not applied directly, understanding spatio-temporal kriging techniques is a plus for fully comprehending spatial data provided by external sources such as Météo France.

In this chapter, the electricity consumption data is available in the same way as the previous chapter: individual aggregations different from those that interest us. We'd like to obtain an estimation of electricity consumption in target zones, not with the auto-regressive structure of the consumption process, but with additional information of the source and target zones.

Introduction

Electric consumption at local level -a city, a village or a block -is an important concern for utilities and grid operators. A grid planner needs load data at a sufficiently small scale to perform spatial load forecasting to optimize maintenance and investments on the grid. With increasingly decentralized generation of renewable energy (for instance, wind or solar power), local consumption is becoming more important for managing the supply-demand balance at the distribution level. Crossed with telecommunication or other utility data at the same level, local electricity consumption can help authorities understand local human activities at a fine-grained temporal level.

In this chapter, we propose a method to estimate past and future electricity consumption at a small temporal scale for target zones with little or no historical consumption data. Two cases are of our interest here: 1. for new substations or feeders, historical data is insufficient for forecasting; 2. for a block or a village, measurements in substations or feeders are not CHAPTER 4. SPATIAL ESTIMATION OF ELECTRICITY CONSUMPTION USING SOCIO-DEMOGRAPHIC INFORMATION q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 45.71 spatially fine-grained enough to supply a direct estimate. To solve these two problems, we can use consumption data in source zones which are similar to the target zones in terms of socio-demographics. An illustration of this target/source zone distinction is given in Figure 4.1

Traditionally, spatial load forecasting is focused on the long term, with power system planning as primary objective [31]. In recent years, thanks to the increasing availability of data, various methods are proposed to model and to forecast electric load at various temporal and spatial scales [33]. For aggregated load forecasting from national down to substation level, stateof-art regression methods achieve satisfactory prediction error, with exogenous variables such as temperature and calendar variables [32,[START_REF] Fan | Short-term load forecasting based on a semi-parametric additive model[END_REF][START_REF] Pierrot | Short-term electricity load forecasting with generalized additive models[END_REF]. However, there is generally a decrease in accuracy of these methods when the aggregation level becomes smaller [START_REF] Mirowski | Demand Forecasting in Smart Grids[END_REF][START_REF] Sevlian | A model for the effect of aggregation on short term load forecasting[END_REF], even when historical data of the target zones are available. Consequently, more complex strategies are proposed to tackle small-area load forecasting. Reference [34] uses the hierarchy structure in electric grid to uncover patterns in the distribution network.

We propose to use socio-demographic information to transport consumption data from source zones to target zones with no or little historical consumption data. The idea is to form clusters of both source and target zones based on client information or socio-demographical information from publicly available datasets. We then estimate semi-parametric consumption models for each cluster, using consumption data only from source zones. These models are then transported to target zones to estimate consumption for past periods, or to forecast future periods. Clustering approaches [35,[START_REF] Alzate | Improved electricity load forecasting via kernel spectral clustering of smart meters[END_REF][START_REF] Rhodes | Clustering analysis of residential electricity demand profiles[END_REF][START_REF] Quilumba | Using Smart Meter Data to Improve the Accuracy of Intraday Load Forecasting Considering Customer Behavior Similarities[END_REF] are employed on smart meter time series to form consumer profiles to improve modeling accuracy. Contrary to these methods, we form clusters using sociodemographics or client information instead of consumption data. This allows us to produce estimates for target zones without historical data. The number of public datasets has been exploding for a few years. Many governments now provide detailed demographic, economic, and 4.2. METHODOLOGY sociological statistics collected at a local scale [72][START_REF]census for England and Wales -Office for National Statistics[END_REF][START_REF]INSEE -Database -Population census results[END_REF]. Our method takes advantage of these new data to produce new small-area load estimates.

The estimation procedure is explained in Section 4.2. The dataset, the validation procedure, and an application using real and simulated consumption is presented in Section 4.3.

Methodology

Our method has two main components: clustering of the geographical zones and consumption models. In the clustering part, we form clusters of both source and target zones using sociodemographic and client information. In the consumption model part, for each of the estimated clusters, we estimate regression models to explain the consumption behavior of the source zones included in this cluster.

Consumption models

The basic building block of our method is a regression model which explains electricity consumption using temperature and calendar variables. The general model is the following:

E(Y t ) = f (Temperature t ) + g(TOY t ) + α SR I SRt + α Xmas I Xmast + 9 h=1 α h I DayType t =h . (4.1)
In this regression model, the expected value of electric consumption at time t, E(Y t ), is the sum of several terms:

• f (Temperature t ): a non-linear function of temperature at time t;

• g(TOY t ): a non-linear function of time of the year (TOY), which is the proportion of the year between January 1 and the day in question (0 on January 1 st and 1 on December 31 st );

• α SR I SRt , α Xmas I Xmast , α h I DayType t =h : three additional dummy variables for periods with special electricity rates (SR), periods included in Christmas holidays (Xmas), and a 9category day type (seven days of the week, plus two types of public holiday). The boldface I is a variable that is equal to 1 when the index is true, and 0 otherwise.

To fit the consumption model on data, we need to specify how the non-linear functions f (Temperature t ) and g(TOY t ) are estimated. We consider two specifications:

• Parametric specification:

-f (Temperature t ) = a 1 h 1 (Temperature t ) + a 2 h 2 (Temperature t )
where h 1 and h 2 are two functions of temperature specific to heating and air conditioning, previously estimated on the French national data [START_REF] Bruhns | A non linear regression model for mid-term load forecasting and improvements in seasonality[END_REF]. Parameters a 1 and a 2 are to be estimated.

g(TOY t ) = m=4 m=1 (b m cos(2πmTOY t ) + c m sin(2πmTOY t )) is the first four terms of a Fourier series to model the yearly pattern of the consumption. Parameters b m and c m , for 1 ≤ m ≤ 4, are to be estimated.

• Generalized additive model (GAM, [START_REF] Wood | Generalized additive models: an introduction with R[END_REF]) specification: we estimate f (Temperature t ) and g(TOY t ) as linear combinations of spline functions. We use thin plate regression splines to estimate f , the temperature function. For g, the time-of-year function, we use cyclic cubic regression splines, which have identical values for up to the second derivative at both ends of a year. We estimate the time-of-year functions for weekdays and weekends separately. In this section, we present the datasets, the validation procedure, and finally the test results of our estimation procedure.

Datasets used

Real electricity consumption We use the hourly consumption during three years (2010-2012) of 473 MV feeders located in the region of interest. Fig. 4.2 shows an example of the time series. For confidentiality reasons, the y-axis is hidden in these figures. The hourly average load of all feeders included is around 1,500 kW over the three years. There is a strong presence of annual, weekly, and daily cycles on these time series, as we can see in the example.

Simulated electricity consumption

To apply our method to IRIS consumption estimation, we generate consumption simulations for the 1,094 IRIS, and the feeders in the region of interest.

The simulations are produced in the following manner. We have the geographical positions of the low-voltage (LV) distribution transformers connected to the feeders. These transformers are therefore positioned inside the boundaries of each IRIS. For each transformer, we choose a client type according to two ratios: the ratio between the number of apartments and contractual power (CP), and the ratio between the number of warehouses and CP. Crossing these two ratios, the transformer is assigned one of the 4 client types (Table 4.1). In addition, we assign a daily peak/off-peak period pattern randomly to each IRIS, from 4 realistic patterns. For each transformer, hourly consumption time series is simulated with regulated French profiles of typical consumers [64]. These profiles are a set of coefficients specifying, for each typical consumer group, the relative hourly consumption throughout a year, and adjustments linked to the temperature. We plug real local temperature and calendar variables for the 2010-2012 period into the profiles to simulate hourly consumption shape of each LV transformer, and multiply it by the total contractual power of the transformer. The IRIS aggregation is calculated by summing over transformers insider its boundary, while the feeder aggregation is obtained by summing over transformers connected to it. SOCIO-DEMOGRAPHIC INFORMATION Socio-demographic characteristics A number of variables are available at IRIS level from the French census [START_REF]INSEE -Database -Population census results[END_REF] and from the SIRENE registry1 . We use the age structure of the population and the number of shops, offices, factories, apartments and warehouses in each IRIS. We distribute the value of each variable for an IRIS in equal parts to the transformers located in it, and then re-aggregate to MV feeders.

Contractual power

We use the contractual power (CP) at each transformer by client class as a proxy to client information. Four types of clients are considered: residential, small business, corporate clients supplied on low voltage, and industrial clients on medium voltage. For each transformer, the sum of CP for each of these four categories is known. The CP at an MV feeder (or an IRIS) is obtained by aggregating that of connected transformers (or within its boundaries).

Temperature Hourly historical temperature is obtained from www.climate.gov for weather stations in the region of interest. Each MV feeder and IRIS is associated to the weather station closest to its served zone.

Validation procedure

We test the estimation procedure by its prediction performance on target zones. The prediction is done in a "batch learning" fashion: the models are estimated using the first 30 months of data at once (training periods) on source zones for each of the clusters. Once the models are estimated, we produce forecasts for the last 3 months (test periods) on target zones, by using the temperature and calendar variables corresponding to test periods. For each target zone the mean absolute percentage error is calculated, according to the definition

M AP E j = 1 n n t=1 |Y j,t -Ŷj,t | Y j,t
, where Y j,t is the consumption of Zone j at Period t, and Ŷj,t the prediction of the procedure. Since the forecasting horizon varies from one hour to three months ahead, the M AP E measures the average error of prediction of the procedure across forecasting horizons. As a comparison, we also predict consumption on test periods for the source zones using the same estimated models. The difference of average M AP E between target and source zones is a measure of the transportability of the models, hence the performance of the procedure.

Application on real MV feeder consumption

On the real consumption dataset, which consists of 473 feeders, we randomly select 100 feeders as source zones, and use the rest as target zones. We conduct 5 independent random selections, and the errors reported are averaged over these 5 runs. On source zones, we run the k-means algorithm on two sets of variables:

• the age structure and number of offices, shops and factors;

• percentage of contractual power in client types.

This results in two sets of clusters (IRIS clusters and CP clusters). For each cluster, we use the 30 months of source zone consumption data to estimate one consumption model. Predictions are made for the last three months on source zones and target zones alike. Jan Mar May Jul Sep Nov Jan -1.0 -0.5 0.0 0.5

IRIS clusters

Effect q q Jan Mar May Jul Sep Nov Jan -1.0 -0.5 0.0 0.5

CP clusters

Effect

IRIS clusters

Effect q q Jan Mar May Jul Sep Nov Jan -0.6 -0.2 0.2 0.6 Figs. 4.3 and 4.4 show the GAM estimation of the annual cycle for consumption at 10 a.m. for each cluster of both sets of clusters, respectively for weekdays and weekends 2,3 . Whether on weekdays or weekends, and across the clusters, the annual cycle is that consumption is higher in winter than in summer, with a large decrease in August due to vacations, and a smaller decrease in May due to the bank holidays. In both sets of clusters, several clusters have a noticeably smaller decrease in August (Clusters 3 and 5 in IRIS data clusters, Clusters 1 and 7 in CP clusters). This is linked to the fact that these clusters are generally located in rural area (i.e. they have a lower than average percentage of apartments), where less people go on vacations in summer. These differences in clusters are also visible in functions estimated in the parametric specification (not shown here). Fig. 4.5 shows the effect of special rate (SR) days estimated by the GAM specification. The higher price starts at 6 a.m. on the SR day, and ends at midnight. This is why these effects start at 6 a.m.. Often clients reduce their consumption during high rate hours and increase consumption just after midnight, when the rates come back to normal, as can be seen in the figures. Clusters 1 and 5 in the IRIS clusters and Cluster 1 in CP clusters have much less drastic SR effect. This shows that the clusters we obtained are able to capture some differences in electricity consumption of the clients, although no consumption data is used in the clustering step. The parametric specification has qualitatively similar results.

CP clusters

Effect

Prediction error rates for the different model specifications are shown in Table 4.2. In forecasts with observed levels (first two rows in each half of the table), the difference of error rate between source and target zones is small. This suggests that the consumption models SOCIO-DEMOGRAPHIC INFORMATION q q 06:00 16:00 02:00 -0.2 0.0 0.2

IRIS clusters

Hour Effect q q 06:00 16:00 02:00 -0.15 0.00 0.15 obtained by the estimation procedure are transportable to target zones. In both cases, the error rate is less than 20%, comparable to the state of art at similar aggregation level with full observations [START_REF] Mirowski | Demand Forecasting in Smart Grids[END_REF].

CP clusters

Hour

As a reference, a preliminary study is done to evaluate of predictions of individual feeder GAM models with the similar specification. The error rates Table 4.2 are higher than these individual models (11.7% for the 473 feeders on average, compared to the around 18% for source zones). This is not surprising, since the cluster models are much less specialized than individual feeder models. This decrease in model accuracy is a small price to pay for the increased model transportability, which enables us to estimate hourly consumption for target zones without using its historical data.

The CP clustering model has a lower error rate, indicating that direct client information has greater importance than indirect socio-demographic variables. Parametric models have similar performances in prediction as GAMs.

We also tried to estimate consumption levels by linear regression using CP variables. The error rates are in the third row of each half of Table 4.2. These predictions have an average error rate between 35% and 40%. This is because the estimation of consumption level by CP is not accurate. If yearly aggregated consumption of target zones is available, this error rate can be reduced.

Application on synthetic consumption of feeders and IRIS

On the synthetic consumption dataset, we keep the same five sets of randomly selected 100 feeders as source zones. As target zones, we use the rest of the feeders, and the 1,094 IRIS.

As in the real consumption application, models are estimated on source zone load of the first 30 months. Predictions are made for the last 3 months on both source and target zones. We only report the results on GAM specification with observed levels, since the relative difference is similar to the real consumption case. To produce clusters, we use either the apartment, warehouse variables used in simulation, or the percentage of CP in client type. The two sets of clusters are similarly called IRIS clusters and CP clusters.

Prediction error rates are shown in Table 4.3. The difference between source and target zones of feeder service area is small, as in the real consumption case. This confirms that consumptions models obtained by the process are indeed transportable to the target zones. The prediction error is generally lower in this case, as consumption is smoother in simulations than in real data, and the clusters are more homogeneous. We observe a similar slight decrease in prediction error when the percentage of CP in client type is used. This is also reassuring, for it suggests a similarity between real data and simulations.

The third row in Table 4.3 shows the prediction error on IRIS consumption. We observe an average error rate of 16.93% and 14.09% for IRIS clusters and CP clusters. This is higher, but of the same order as short-term forecasting at this aggregation level with full data [START_REF] Mirowski | Demand Forecasting in Smart Grids[END_REF]. Such precision is reasonably good, given that the estimation requires only the average consumption level for each IRIS. As a benchmark, we also used an interpolation proportional to the population from the 100 source zones to produce an estimation for aggregated consumption of the IRIS. This method does very poorly on aggregated IRIS consumption, with 49.57% of mean error. Our method is thus significantly better than the interpolation approach.

Conclusion

In this chapter, we propose a method to estimate hourly electricity consumption for target zones with insufficient or no historical consumption data using measurements from source zones and public socio-demographic information. The procedure provides state-of-art prediction performance on real and simulated datasets for target zones with very little data, and significantly outperforms the benchmark interpolation method in simulations. Introduction

Introduction

As is suggested by Figure 2.1 on page 14, the object of interest in this thesis can easily be represented in a matrix form. Based on this data representation, we will use nonnegative matrix factorization (NMF, [5]) methods to tackle the estimation problems introduced in Chapter 2. Matrix factorization methods leverage the fact that big matrices observed in real datasets are often of rank much smaller than their dimension, to facilitate the estimation, compression, and processing of such datasets.

In this chapter, after a short introduction to NMF (Section 5.2.1), we will reformulate the estimation problems discussed in Chapter 2 in terms of matrix factorization (Section 5.2.2). In Section 5.3, we will examine some recent progress in the non-convex optimization literature, particularly, how non-convex problems such as matrix factorization can have global convergence with an iterative algorithm, and why it is difficult to achieve the same in NMF.

In Chapters 6 and 7, we will study two specific cases of the general model to better adapt it to our estimation problems. In Chapter 8, we will compare the NMF framework with kriging and socio-demographic clustering introduced previously in Part I (Chapters 3 and 4). We will see that these two methods can be reformulated as special cases of the matrix factorization framework. In Chapter 9, we will review some particularities in the implementation of the CHAPTER 5. INTRODUCTION TO NONNEGATIVE MATRIX FACTORIZATION algorithms developed in this part, and present meterModels, the R package developed in this thesis.

Multivariate time series estimation and prediction as non-

negative matrix recovery

Nonnegative matrix factorization

Low-rank approximation and factorization methods have been extensively studied in recent years. This class of methods uses the low-rank hypothesis: that a matrix V ∈ R n 1 ×n 2 is of rank k much smaller than its dimensions, with k n 1 , n 2 . This means that V can be rewritten in a factorized form:

V = F r F T c , ( 5.1)
where F r ∈ R n 1 ×k and F c ∈ R n 2 ×k are called factors. In particular, by doing a singular value decomposition (SVD) of V, one can easily find such a low-rank factorization.

It is in general interesting to use F r and F c instead of V. Since they have lower dimensions, they can be useful for a wide range of applications such as estimation or compression. Even if V is not exactly of a small rank k, it can still be interesting to find a rank-k approximation for it, by minimizing min

Fr∈R n 1 ×k ,Fc∈R n 2 ×k V -F r F T c 2 F , ( 5.2)
for example, where for any matrix X, X F is its Frobenius norm (Euclidean norm of its entries).

In this case, the minimization problem (5.2) is called low-rank approximation.

In this thesis, we are especially interested in nonnegative matrix factorization, where the matrix to be factorized and the factors only have nonnegative entries. Note that a matrix of rank k does not necessarily have a nonnegative matrix factorization with k columns in the factor matrices. The smallest integer k for which there exist

F r ∈ R n 1 ×k + and F c ∈ R n 2 ×k + so that V = F r F T c is called the nonnegative rank of V.
The nonnegative rank is always equal or larger than the rank.

For a given integer k, the nonnegative factorization/approximation problem then becomes min

Fr∈R n 1 ×k ,Fc∈R n 2 ×k V -F r F T c 2 F s.t. F r ≥ 0, F c ≥ 0, (5.3)
where X ≥ 0 for a matrix X is an entry-wise inequality.

Nonnegative matrix factorization is generally more interpretable, when the underlying quantites are indeed nonnegative. In our application case, the entries of the matrix to be recovered are the quantity of electricity consumption, which is positive. In the convention adopted in this thesis, the rows are periods, and columns are individuals (Figure 2.1). The factor matrices have a naturel interpretation in this convention: the left factor F r has columns which are interpreted as typical consumption profiles, and the entries of the right factor matrix F c are weights of the individuals in each of the profiles (Figure 5.1).

General model definition

To solve the estimation problems considered in Figure 2.1, we are interested in reconstructing a nonnegative matrix

V * ∈ R n 1 ×n 2 + , from N linear measurements, α = A(V * ) ∈ R N , ( 5.4)
where A : R n 1 ×n 2 → R N is a linear operator. Formally, A can be represented by A 1 , ..., A N , N design matrices of dimension n 1 × n 2 , and each linear measurement can be represented by We suppose that the matrix of interest, V * , is of nonnegative rank k: we can find two nonnegative matrices

α i = Tr(V * A T i ) = V * , A i . ( 5
F r ∈ R n 1 ×k + and F c ∈ R n 2 ×k + so that V * = F r F T c , with k n 1 , n 2 .
Note that this implies that V * is of rank at most k, and therefore is of low rank.

In the nonnegative matrix recovery problem, we try to recover the true matrix V * as well as the factor matrices F r and F c , given A, the measurement operator, and α, the measurements.

To obtain an estimation, we minimize the quadratic error of the matrix factorization. In the two following chapters, we will propose algorithms for optimization problems of the following form: min

Fr∈R n 1 ×k ,Fc∈R n 2 ×k ,V∈R n 1 ×n 2 V -F r F T c 2 F s.t. F r ≥ 0, F c ≥ 0, V ≥ 0, A(V) = α.
(5.6)

Note that in this thesis, α, the data used in matrix recovery are generated by a measurement operator, hence are "coherent". Since there is actually a ground truth matrix V * that verifies A(V * ) = α, the feasible set is obviously non-empty. In order to recover a low-rank matrix, the following alternative optimization problem can also be considered: min

Fr∈R n 1 ×k ,Fc∈R n 2 ×k α -A(F r F T c ) 2 2 s.t. F r ≥ 0, F c ≥ 0.
(5.7)

Instead of minimizing the low-rank approximation error for a matrix that verifies the data constraint in (5.6), (5.7) minimizes the sampling error of an exactly low-rank matrix. Both have been studied in the literature. For example, objective functions similar to (5.6) have been considered in [78], and ones similar to (5.7) in [79]. We will compare the two in Chapter 7 to argue that (5.6) is a more efficient strategy than (5.7).

Applications of matrix recovery

By specifying different masks, matrix recovery has a number of interesting applications, old and new.

• Complete observation:

N = n 1 n 2 , A i 1 ,i 2 = e i 1 e T i 2
, where e i is the i-th canonical vector. This means every entry of V * is observed. Obviously, in this case, there is no need to estimate V * . The main reason to do this is for dimension reduction or denoising. With nonnegative factors, this is used in image segmentation [5], topic extraction [17], and clustering [START_REF] Ding | Orthogonal Nonnegative Matrix T-Factorizations for Clustering[END_REF], for example.

• Matrix completion -recommender systems: N < n 1 n 2 , the set of design matrix is a subset of complete observation masks. The objective is then to infer the unobserved entries. This is the standard task for recommender systems [START_REF] Yu | Scalable Coordinate Descent Approaches to Parallel Matrix Factorization for Recommender Systems[END_REF], whether with nonnegative factors or real factors.

• Matrix sensing: the design matrices A i are i.i.d. (independent and identically distributed) random variables from a certain random matrix probability distribution. Typi-cally, the probability distribution needs to verify certain conditions, so that with a large probability, A verifies the Restricted Isometry Property [78]. This class of masks are very well studied theoretically. For a more recent review, see [START_REF] Davenport | An Overview of Low-Rank Matrix Recovery from Incomplete Observations[END_REF]. An application of this class of masks with nonnegative factors in imaging can be found in [83].

• Rank-one projections: the design matrices are random rank-one matrices, that is

A i = α i β T i
, where α i and β i are respectively random vectors of dimension n 1 and n 2 . The main advantage to this setting is that much less memory is needed to store the masks, since we can store α i and

β i (dimension-(n 1 + n 2 )) instead of A i (dimension-(n 1 × n 2 )).
In [84,85], theoretical properties are proved for the case where α i and β i are vectors with independent Gaussian entries and/or drawn uniformly from the vectors of the canonical basis.

• Temporal aggregate measurements: in this case, the matrix is composed of n 1 time series concerning n 2 individuals, and each measure is a disjoint temporal aggregate of the time series of an individual. The design matrices are defined as

A i = t 0 (i)+h(i)
t=t 0 (i)+1 e t e T s i , where s i is the individual concerned by the i-th measure, t 0 (i)+1 the first period covered by the measure, and h(i) the number of periods covered by the measure. This class of masks is the main topic of Chapter 6, since it models the meter readings of resident clients with traditional electricity meters. A similar type of masks has also been studied in network traffic matrix estimation [79].

Convergence of nonnegative matrix factorization algorithms

Global convergence of alternating minimization algorithms in non-convex problems

The general matrix factorization problem (5.2) has a simple form, but finding its solutions is not simple. Methods based on nuclear norm minimization (by minimizing the trace of V) for matrix recovery and matrix completion problems often have theoretical guarantees [78,86]. However, first-order algorithms based on the factorized form (alternating minimization, gradient descent, etc.) are much more efficient, and more preferred by practitioners, because they work on a smaller parameter space. This class of methods try to solve a non-convex problem, and therefore rarely had global convergence guarantees.

It is only recently that theoretical guarantees for general matrix factorization algorithms began appearing [START_REF] Jain | Low-Rank Matrix Completion Using Alternating Minimization[END_REF][START_REF] Sun | Guaranteed Matrix Completion via Nonconvex Factorization[END_REF][START_REF] Ge | Matrix Completion Has No Spurious Local Minimum[END_REF][90][91][START_REF] Park | Non-Square Matrix Sensing without Spurious Local Minima via the Burer-Monteiro Approach[END_REF][START_REF] Li | Symmetry, Saddle Points, and Global Geometry of Nonconvex Matrix Factorization[END_REF]. This progress is connected to the literature on global convergence results for other non-convex problems [16,[START_REF] Choromanska | The Loss Surfaces of Multilayer Networks[END_REF][START_REF] Dauphin | Identifying and Attacking the Saddle Point Problem in High-Dimensional Non-Convex Optimization[END_REF][START_REF] Afonso | On the Low-Rank Approach for Semidefinite Programs Arising in Synchronization and Community Detection[END_REF][START_REF] Mei | The Landscape of Empirical Risk for Non-Convex Losses[END_REF]. Motivated by the superior empirical performances, this trend of analysis shows that although a specific class of problems is not convex, nice qualities about the landscape can still guarantee the convergence of descent algorithms.

Of particular interest is the reference [90], in which the authors showed that the general matrix factorization/recovery problem has no spurious local minima. Combined with results on escaping saddle points, stochastic gradient descent algorithms therefore can converge to the global optimum. Contrary to a number of previous works, this result does not rely on a particular initialization of the algorithm.

In general, NMF is an NP-hard problem [11]. Therefore, there cannot be similar results to the general NMF problem. However, in [91] and [98], guarantees of the NMF have been given for alternating direction algorithms, for cases where one of the factors is initialized at a point close to the ground truth.

In the rest of this section, we will go through the steps by [90] to show why this analysis cannot be directly applied to NMF.

What would global convergence mean?

In order to simplify the theoretical analysis, we focus temporarily on the factorization case (5.3), where the data matrix V * is completely observed. We note by (F r ,

F c ) = V * -F r (F c ) T 2
F , the loss function of (5.3).

A pair of matrices (F * r , F * c ) is called a local minimizer of (5.3), if there exists > 0, so that

∀(Z 1 , Z 2 ) ∈ R n 1 ×k + × R n 2 ×k + with Z 1 -F * r 2 F + Z 2 -F * c 2 F ≤ , we have V * -F * r (F * c ) T 2 F ≤ V * -Z 1 (Z 2 ) T 2 F .
Any local minimizer satisfies the Karush-Kuhn-Tucker (KKT) necessary conditions [START_REF] Bertsekas | Nonlinear Programming[END_REF]. Noting by ∂ Fr (F r , F c ) and ∂ Fc (F r , F c ) the partial derivatives of the loss function, and by ∇ 2 (F r , F c ) the second order derivative, the KKT conditions are,

∂ Fr (F * r , F * c ) ≥ 0, ∂ Fr (F * r , F * c ) • F * r = 0, ∂ Fc (F * r , F * c ) ≥ 0, ∂ Fc (F * r , F * c ) • F * c = 0, (5.8)
and

Z 1 Z 2 , ∇ 2 (F * r , F * c ) Z 1 Z 2 ≥ 0, ∀ Z 1 Z 2 ∈ R (M +N )×K , that satisfies Z 1 • 1 Fr=0 = 0, Z 2 • 1 Fc=0 = 0, (5.9)
where • denotes the element-wise matrix product, and 1 X=0 is a matrix of the same dimension of X for any matrix X, with an entry 1 where X has an entry 0, and 0 otherwise. This means we only need to check for second-order increments in directions in which F r and F c are not on the boundary of the first orthant.

If (F * r , F * c ) = 0 is true for all matrix pairs (F * r , F * c ) that verify the KKT conditions, then all local minima of the optimization problem above is a global optimum. If, additionally, the matrix has a unique NMF up to scaling and permutation, then all local minima of the optimization program is the NMF of the matrix V * .

Therefore, what we would like to know is if there are sufficient conditions on V * that an algorithm could check with relative ease, for the KKT conditions to imply zero residual ( (F * r , F * c ) = 0).

Link with general matrix factorization

The question of global optimality of local NMF algorithms arises naturally from the work of [90], who proved that the general matrix factorization problem (namely minimization of the same loss function without non-negativity constraint on the factors) has no spurious local minima. Their work is in a more general context where the data are linear measures of V under a measurement operator which satisfies a restricted isometry property (RIP).

There are two-ways to view the link between general matrix factorization and NMF:

• NMF is the optimization of the same objective function as general matrix factorization, in a constrained domain (local minima could therefore exist on the boundary of the feasible domain of NMF). ) is also a global optimum. NMF is about finding out a specific invertible matrix R so that both F r R and F c (R -1 ) T are non-negative. Intuitively, this could be quite difficult, since there is no information on R in the data.

CONVERGENCE OF NONNEGATIVE MATRIX FACTORIZATION ALGORITHMS

Empirical convergence

Although there is currently no prove of global convergence for local NMF algorithms, empirically there seems to be global convergence on certain class of NMF problems, especially if the true factors have many zero entries.

In Tables 5.1 and 5.2, we report the NMF objective function value and the relative error on the obtained factors F c and F r , on matrices with unique NMFs with rank ranging from 4 to 14. The algorithms used (NeNMF and HALS) are classical NMF algorithms [20,27].

The original and estimated factors, and the error on each element, are plotted in Figures 5.3, 5.4, 5.5, 5.6. The relative error on factors is of several percentage points, and we can see that estimated and original matrices are quite close.

Local optimality conditions

By differentiating the loss function, we obtain that

∂ Fr (F r , F c ) = (F r F T c -V * )F c , ∂ Fc (F r , F c ) = (F c F T r -(V * ) T )F r .
By differentiating a second time, we obtain that

∀(Z 1 , Z 2 ) ∈ R n 1 ×k × R n 2 ×k , ∇ 2 (F r , F c ) Z 1 Z 2 = lim →0 1 ∂ Fr (F r + Z 1 , F c + Z 2 ) -∂ Fr (F r , F c ) ∂ Fr (F r + Z 1 , F c + Z 2 ) -∂ Fc (F r , F c ) = (F r F T c -V * )Z 2 + (F r Z T 2 + Z 1 F T c )F c (F c F T r -(V * ) T )Z 1 + (F c Z T 1 + Z 2 F T r )F r .
This leads to

Z 1 Z 2 , ∇ 2 (F r , F c ) Z 1 Z 2 = Z 1 , (F r F T c -V * )Z 2 + (F r Z T 2 + Z 1 F T c )F c + Z 2 , (F c F T r -(V * ) T )Z 1 + (F c Z T 1 + Z 2 F T r )F r =2 Z 1 Z T 2 , (F r F T c -V * ) + Z 1 F T c + F r Z T 2 , F r Z T 2 + Z 1 F T c =2 Z 1 Z T 2 , (F r F T c -V * ) + Z 1 F T c + F r Z T 2 2 F .
Therefore, the KKT conditions for the NMF problem are

(F r F T c -V * )F c ≥ 0, F r ≥ 0, (F r F T c -V * )F c • F r = 0, (F c F T r -(V * ) T )F r ≥ 0, F c ≥ 0, (F c F T r -(V * ) T )F r • F c = 0, 2 Z 1 Z T 2 , (F r F T c -V * ) + Z 1 F T c + F r Z T 2 2 F ≥ 0, ∀ Z 1 Z 2 ∈ R (n 1 +n 2 )×k , that satisfies Z 1 • 1 Fr=0 = 0, Z 2 • 1 Fc=0 = 0.
Define the linear function f Fr,Fc : R

n 1 ×k × R n 2 ×k → R n 1 ×n 2 , where f Fr,Fc (Z 1 , Z 2 ) = Z 1 F T c + F r Z T 2 . Consider the following linear subspace of R n 1 ×k × R n 2 ×k : e(F r , F c ) ≡ {(Z 1 , Z 2 )|Z 1 • 1 Fr=0 = 0, Z 2 • 1 Fc=0 = 0}.
and the following linear subspace of R n 1 ×n 2 :

E(F r , F c ) ≡ {Z 1 F T c + F r Z T 2 |(Z 1 , Z 2 ) ∈ e(F r , F c )},
Obviously, f Fr,Fc (e(F r , F c )) = E(F r , F c ), and e(F r , F c ) ⊂ f -1 Fr,Fc (E(F r , F c )). Essentially, if (Z 1 , Z 2 ) ∈ e(F r , F c ) that Z 1 and Z 2 have zero entries, wherever F r and F c have zero entries. In particular, this means

(F r F T c -V * )F c • Z 1 = 0, (F c F T r -(V * ) T )F r • Z 2 = 0.
Note that f Fr,Fc is not injective on e(F r , F c ). Especially, for all 1 ≤ i ≤ k, note by x i the matrix with the i-th column equal to that of X and 0 everywhere else, then (f c,i , -f r,i ) ∈ Ker(f Fr,Fc ).

Projections

If (F r , F c ) satisfies the first-order optimality condition, then F r F T c is the projection of

V * in E(F r , F c ). This is true, because F r F T c ∈ E(F r , F c ), and V * -F r F T c ∈ E(F r , F c ) ⊥ . Note this fact by P E(Fr,Fc) (V * ) = F r F T c .
In fact, if (F r , F c ) satisfies the first-order optimality condition, then F r F T c is the projection of V * in a even larger set. The following closed convex cone contains E(F r , F c ):

F (F r , F c ) ≡ {Z 1 F T c + F r Z T 2 |(Z 1 , Z 2 ) ∈ R n 1 ×k × R n 2 ×k , Z 1 , (F r F T c -V * )F c ≥ 0, Z 2 , (F c F T r -(V * ) T )F r ≥ 0}. By definition of F (F r , F c ), if Z 1 F T c + F r Z T 2 ∈ F (F r , F c ), then Z 1 F T c + F r Z T 2 , V * -F r F T c ≤ 0. This means that F r F T c is in the polar cone of F (F r , F c ) (noted as F r F T c ∈ F (F r , F c ) , [100, Chap.6]). Together with F r F T c ∈ F (F r , F c ), and V * -F r F T c ⊥F r F T c , this constitutes a charac- terization of projection of V * in F (F r , F c ). That is P F (Fr,Fc) (V * ) = F r F T c . Define Π(F * r , F * c ) = P e(Fr,Fc) (F * r Π, F * c Π)
, for any K-by-K permutation matrix Π ( Π is the composition of the permutation with projection onto the support of (F r , F c ). I will also use Π(F * r ) (and Π(F * c )) to note the projection of F * r (and F * c ) onto the appropriate support, by an abuse of notation.

Consider the following subset of R n 1 ×n 2 :

{Z 1 Z T 2 |(Z 1 , Z 2 ) ∈ e(F r , F c )}.
This set contains Π(F * r , F * c ), but is in general not convex. This makes it difficult to consider the projection of V * in it. On the other hand,

G(F r , F c ) ≡ conv{Z 1 Z T 2 |(Z 1 , Z 2 ) ∈ e(F r , F c )}, 64 5.3

. CONVERGENCE OF NONNEGATIVE MATRIX FACTORIZATION ALGORITHMS

the convex hull of the previous subset, is a convex cone. In fact

G(F r , F c ) ⊂ {X ∈ R n 1 ×n 2 |X • (F r F T c -V * )F c F T r (F r F T c -V * ) = 0}. Since both (F r F T c -V * )F c and F T r (F r F T c -V * ) are nonnegative, the support of the product is {(i, j)|supp((F r F T c -V * ) i F c ) supp((F c F T r -(V * ) T ) j F r ) = ∅}. Therefore, any element of G(F r , F c ) has a support of empty intersection with (F r F T c -V * )F c F T r (F r F T c -V * ). We have P G(Fr,Fc) (V * ) -Π(F * r ) Π(F * c ) T ≥ 0.

Applying second order condition

Now suppose that (F r , F c ) is a local minimum, i.e., it verifies both (5.8) and (5.9).

Note that if P

G(Fr,Fc) (V * ) = Π(F * r ) Π(F * c ) T , then Π(F * r ) Π(F * c ) T -V * = P G(Fr,Fc) ⊥ (V *
). If this is the case, we can prove the following inequality, which is similar to the second inequality in [90, Corollary4.1]:

k j=1 (F r -Π(F * r ))e j e T j F T c + F r e j e T j (F c -Π(F * c )) T 2 F ≥ 2 F r F T c -Π(F * r ) Π(F * c ) T 2 F . (5.10)
To see why (5.10) is true, consider the second-order condition with

Z 1 = (F r -Π(F * r ))e j and Z 2 = (F c -Π(F * c ))e j , for 1 ≤ j ≤ k. Therefore, 2 (F r -Π(F * r ))e j e T j (F c -Π(F * c )) T , (F r F T c -V * ) + (F r -Π(F * r ))e j e T j F T c + F r e j e T j (F c -Π(F * c )) T 2 F ≥ 0.
We know that

(F r -Π(F * r ))e j e T j (F c -Π(F * c )) T , F r F T c -V * = F r e j e T j F c + Π(F * r )e j e T j Π(F * c ) T , F r F T c -V * = Π(F * r )e j e T j Π(F * c ) T -F r e j e T j F c , F r F T c -V * = Π(F * r )e j e T j Π(F * c ) T -F r e j e T j F c , F r F T c -Π(F * r ) Π(F * c ) T + Π(F * r )e j e T j Π(F * c ) T -F r e j e T j F c , Π(F * r ) Π(F * c ) T -V * .
Since Π(F * r )e j e T j Π(F * c ) T -F r e j e T j F T c ∈ G(F r , F c ), the second term is zero. In this case, we have

K j=1 (F r -Π(F * r ))e j e T j F T c + F r e j e T j (F c -Π(F * c )) T 2 F -2 K j=1 F r e j e T j F c -Π(F * r )e j e T j Π(F * c ) T , F r F T c -Π(F * r ) Π(F * c ) T = K j=1 (F r -Π(F * r ))e j e T j F T c + F r e j e T j (F c -Π(F * c )) T 2 F -2 F r F T c -Π(F * r ) Π(F * c ) T 2 F ≥ 0,
which then leads to (5.10).

For any local minimum, (5.10) bounds the approximation error (right-hand side) by a quantity connected to the error in the factor space (left-hand side), regardless of the permutation of columns in the factor space.

In [90], once the similar inequality [90, Corollary4.1] is established, an upper bound is then proposed for the left-hand side, so as to prove that the global error is actually zero. Therefore the local minimum is actually a global minimum as well. This is achieved using the properties of the orthogonal matrix which minimizes the left-hand side. Here instead of considering an orthogonal matrix, we use a permutation matrix Π, since we are only allowed to consider points in the nonnegative orthant. It is not clear how we can translate this into the nonnegative case, since we would need Π(F * r )R to be nonnegative as well, therefore the properties of R would change.

Introduction

In this chapter, we propose a new matrix recovery method using nonnegative matrix factorization (NMF, [4]) where matrix columns represent time series at a fine temporal scale. Moreover, only temporal aggregates of these time series are observed.

As explained in previous chapters, the motivation comes from the context of electricity load balancing, where time series represent electric power consumption. To avoid failure in the electricity network, suppliers are typically required by transmission system operators (TSO) to supply as much electricity as their consumers consume at every moment. This mechanism is called balancing. In the context of an open electricity market, all market participants, such as suppliers, utility traders, and large consumers, have a balance responsibility: any imbalance caused within the perimeter of a participant is billed by the TSO. To calculate the imbalance caused by a market participant, one needs an estimation of the consumption and production within its perimeter at a small temporal scale, for example, half-hourly ( [10], [START_REF] Svk | Balance Responsibility[END_REF], [START_REF] Ree | Balance Responsible Party[END_REF]).

CHAPTER 6. NONNEGATIVE MATRIX FACTORIZATION FOR TIME SERIES RECOVERY FROM

A FEW TEMPORAL AGGREGATES However, for many customers (for instance residential) within a perimeter, electricity consumption is not recorded at that scale. Although smart meter readings may be recorded locally up to every minute, utility companies often have very limited access to such data, due to data transmission and processing costs and/or privacy issues. Following a fixed schedule, cumulative consumption of each meter is recorded by the utility company, for instance every day or every month. By differentiating consecutive readings, the utility obtains the consumption of a customer between two reading dates. Currently, TSOs use proportional rules to reconstitute consumption from these measurements, based on national consumption profiles adjusted by temperature. In this chapter, we develop an NMF-based matrix recovery method providing a solution to consumption reconstitution from such temporal aggregates.

Recent advances in matrix completion have made it clear that when a large number of individuals and features are involved, even partial data could be enough to recover much of lost information, thanks to the low-rank property [86]: although the whole data matrix

V * ∈ R n 1 ×n 2 is only partially known, if V * = F r F T c , where F r ∈ R n 1 ×k , F c ∈ R n 2 ×k
, with k much smaller than both n 1 and n 2 , one could recover V * entirely under some conditions over the sampling process.

Temporal aggregated electricity consumption

In this chapter, we address electricity consumption reconstitution as a matrix recovery problem. Consider the electricity consumption of n 2 consumers during n 1 periods. Since consumption is always positive, the n 2 time series are organized into a nonnegative matrix V * ∈ R n 1 ×n 2

+

. An entry of this matrix, v * t,n represents, for example, the electricity consumption of Consumer n for Period t.

Information about consumption is revealed as meter readings which do not correspond to matrix entries but to cumulative sums of each column of V * : at a meter-reading date t, we observe that Consumer n has consumed t i=1 v * i,n since the first period. Several readings are available for each consumer.

An alternative matrix representation could be to define entries directly as the cumulative consumption since the first period. Again, this matrix has missing values and a matrix completion algorithm can be applied. However, this cumulative matrix has increasing columns, which is quite different from matrices considered in the standard matrix completion literature, where matrix completion error is typically bounded by the largest entry of the matrix.

We represent meter readings as linear measures on the consumption matrix V * . Temporal aggregates are derived from meter readings by differentiating consecutive meter readings: we will consider these temporal aggregates as "observations" in the rest of the chapter. We consider D scalar observations, represented by a data vector α ≡ A(V * ) ∈ R D + , where A is a D-dimensional linear operator. To recover V * from α, we look for a low-rank NMF of

V * : F r F T c V * , where F r ∈ R n 1 ×k + , F c ∈ R n 2 ×k

+

. The columns of F r are k nonnegative factors, which can be interpreted as typical profiles of the time series, and the columns of F c as the weights of each individual on each profile. The problem is formalized as the minimization of a quadratic loss function under nonnegativity and data constraints, as in (5.6): min

V, Fr, Fc (V, F r , F c ) = V -F r F T c 2 F s.t. V ≥ 0, F r ≥ 0, F c ≥ 0, A(V) = α, (6.1)
where X ≥ 0 (or x ≥ 0) means that the matrix X (or the vector x) is element-wise nonnegative.

INTRODUCTION

Note the difference between (6.1) and another potential estimator, as in (5.7), min

Fr, Fc ˆ (F r , F c ) = A(F r F T c ) -α 2 2 s.t. F r ≥ 0, F c ≥ 0, (6.2)
studied in [79]. If V is a solution to (6.1), it satisfies exactly the measurement constraint, but is approximately low-rank, while F r F T c , a solution to (6.2) is exactly of low rank, but only matches the measurements approximately. Since in our application, the estimated time series matrix is to be used for billing, the match to metering data is essential. Therefore, we use (6.1) in this chapter.

Autocorrelation

Electricity consumption time series has the nice property of having high temporal autocorrelation. In order to better recovery consumption from temporal aggregates, we propose to exploit this information.

We suppose that the temporal autocorrelation of each individual considered in the dataset is available to us. In practice, this information could be obtained by exploiting historical data of each client. It is also possible to use the results of some prior modeling, which estimates the autocorrelation of each class of clients.

We add a penalty term to the base algorithm to take advantage of this information. Efficient computation is then deduce by using recent convex relaxation of quadratically constrained quadratic programs [START_REF] Ben | Hidden Conic Quadratic Representation of Some Nonconvex Quadratic Optimization Problems[END_REF].

Prior works

The measurement operator A considered here is a special instance of the trace regression model [104] which generalizes the matrix completion setting. In matrix completion, each measurement is exactly one entry. Various forms of linear measurements other than matrix completion have been considered for matrix recovery without nonnegativity [78,85,105].

The NMF literature is generally focused on full observation [6,[START_REF] Alquier | An Oracle Inequality for Quasi-Bayesian Non-Negative Matrix Factorization[END_REF], or on matrix completion [START_REF] Gillis | Low-Rank Matrix Approximation with Weights or Missing Data Is NP-Hard[END_REF][START_REF] Xu | A Block Coordinate Descent Method for Regularized Multiconvex Optimization with Applications to Nonnegative Tensor Factorization and Completion[END_REF] Random projection measurements are used in an NMF context in [START_REF] Pnevmatikakis | Sparse Nonnegative Deconvolution for Compressive Calcium Imaging: Algorithms and Phase Transitions[END_REF], where a maximum likelihood estimator is developed based on a specific generative model in neural imaging. The particular form of measurement operator considered here arises from meter reading, and can be used in other fields, such as Internet traffic matrix estimation [79]. Because of our choice of estimator (6.1) over estimator (6.2), we derive a novel algorithm for this measurement operator, which has a smaller time complexity than previously studied ones (more details in Section 6.2.1).

The idea of imposing time series structure on matrix factorization is not new. Previous approaches combining matrix factorization and autoregressive structure are often focused on obtaining factors that are more smooth and/or sparse, both in NMF [22][23][24] and without nonnegativity [25,26]. Our objective is different from these studies: we try to further improve the matrix recovery by constraining temporal correlation on individual time series (not factors). We use a recent convex relaxation of quadratically constrained quadratic programs [START_REF] Ben | Hidden Conic Quadratic Representation of Some Nonconvex Quadratic Optimization Problems[END_REF] to deduce a closed-form projection step in this case.

We propose an algorithm to solve (6.1) in Section 6.2.1. To take into account individual autocorrelation, a second algorithm is proposed in Section 6.2.2. In Section 6.3, both algorithms are validated on synthetic and real electricity consumption datasets, compared to a linear benchmark and a state-of-art matrix completion method.

Reconstitution of time series with NMF

Iterative algorithm with simplex projection

We represent temporal aggregation by a linear operator A. For each 1 ≤ d ≤ D, the d-th measurement on X, A(X) d , is the sum of several consecutive rows on one column of X, that is,

A(X) d = (t,n)∈I d x t,n ,
where Note that although we are focused on temporal aggregates in this chapter, when applicable, the results announced can also be generalized to general linear measurements (see Chapter 7 for more details).

I d = {(t, n)|t 0 (d) + 1 ≤ t ≤ t 0 (d) + h(d), n = n d },
A block Gauss-Seidel algorithm (Algorithm 1) is used to solve (6.1). We alternate by minimizing (V, F r , F c ) over F r , F c or V, keeping the other two matrices fixed. Methods from classical NMF problems are used to update F r and F c [18]. In the implementation, we use two variants that seem similarly efficient (more details in Section 6.3): Hierarchical Alternating Least Squares (HALS, [19]), and a matrix-base NMF solver with Nesterov-type acceleration (NeNMF, [20]).

When F r and F c are fixed, the optimization problem on V is equivalent to D simplex projection problems, one for each scalar measurement. For 1 ≤ d ≤ D, we have to solve

min v I d v I d -(F r ) t 0 (d)+1≤t≤t 0 (d)+h(d) (F c ) n d 2 s.t. v I d ≥ 0, v T I d 1 = α d , (6.3)
where (F r ) t 0 (d)+1≤t≤t 0 (d)+h(d) is the h(d) rows of F r which correspond to the periods of the measurement, and (F c ) n d is the n d -th row of F c stacked as a column vector. The simplex projection algorithm introduced by [21] solves this subproblem efficiently. Define the operator, P A , as the orthogonal projection into the simplex A ≡ {X ∈ R n1 ×n 2 + |A(X) = α}. The simplex A is the intersection of the affine subspace {X ∈ R n 1 ×n 2 |A(X) = α} and the first orthant. Projector P A encodes the measurement data α = A(V * ). In Algorithm 1, we apply P A to a working copy of V in order to obtain its projection in A.

Contrary to previously studied algorithms [79], by choosing estimator (6.1) over (6.2), the simplex projection step is separated from the classical NMF update steps in our algorithm. Instead of multiplying the rank and the complexity introduced by the number of measurements, we have an algorithm whose complexity is the sum of the two. In cases where the number of measurements is large, this difference can be crucial 1 .

A classical stopping criterion in the NMF literature is based on Karush-Kuhn-Tucker (KKT) conditions on (6.1) [6,Section 3.1.7]. We calculate

R(F r ) i,j = |(F r F T c -V)F T c ) i,j |1 (Fr) i,j =0 , and R(F c ) i,j = |((F c F T r -V T )F T r ) i,j |1 (Fc) i,j =0 . The algorithm is stopped if R(F r ) 2 F + R(F c ) 2
F ≤ , for a small threshold > 0.

RECONSTITUTION OF TIME SERIES WITH NMF

Algorithm 1 Block coordinate descent for NMF from temporal aggregates

input P A , 1 ≤ K ≤ min{T, N } Initialize F 0 r , F 0 c ≥ 0, V 0 = P A (F 0 r F 0 c ), i = 0 while Stopping criterion is not satisfied do F i+1 r = Update(F i r , F i c , V i ) F i+1 c = Update(F i+1 r , F i c , V i ) V i+1 = P A (F i+1 r F i+1 c ) i = i + 1 end while output V i ∈ A, F i r ∈ R n 1 ×k + , F i c ∈ R n 2 ×k

+

Convergence to a stationary point has been proved for past NMF solvers with the full observation or the matrix completion setting [18,20]. Although the subproblems on F r and F c do not necessarily have unique optimum, the projection of V attains a unique minimizer. By [110,Proposition 5], all limiting points of Algorithm 1 is a stationary point.

Algorithm 1 can be generalized to other types of measurement operators A, as long as a projection into the simplex defined by the data constraint A(X) = α and the positivity constraint can be efficiently computed.

From autocorrelation constraint to penalization

In addition to the measurements in α, we have some prior knowledge on the temporal autocorrelation of the individuals. To take into account information about autocorrelation, we add a penalization term to the original matrix recovery problem, replacing (6.1) by: min (6.4) where λ ≥ 0 is a single fixed penalization parameter, and ∆ ρn is a symmetric matrix described shortly after. In the rest of this section, we first show by Theorem 1, that with an appropriately chosen value of λ, adding the penalization term v T n ∆ ρn v n is equivalent to impose that the temporal autocorrelation of v n to be at least equal to ρ n , a prior threshold. Then we modify Algorithm 1 to solve this penalized problem.

V,Fr,Fc V -F r F T c 2 F -λ N n=1 v T n ∆ ρn v n s.t. V ≥ 0, F r ≥ 0, F c ≥ 0, A(V) = α,
For 1 ≤ n ≤ n 2 , suppose that the lag-1 autocorrelation of Individual n's time series is at least equal to a threshold ρ n (e.g. from historical data, excluded from observed temporal aggregates), that is,

n 1 -1 t=1 v t+1,n v t,n ≥ ρ n n 1 t=1 v 2 t,n . (6.5)
Notice that with the lag matrix,

∆ =          0
         , we have n 1 -1 t=1 v t+1,n v t,n = v T n ∆v n . Define ∆ ρ ≡ ∆ + ∆ T -2ρI, for a threshold -1 ≤ ρ ≤ 1. Inequality (6.5) is then equivalent to v T n ∆ ρn v n ≥ 0. (6.6)
Imposing (6.6) would require one to solve, at each iteration, n 2 quadratically constrained quadratic programs (QCQP) of the form:

min x ||x -x 0 || 2 s.t. x T Sx ≥ 0, (6.7)
where S is a general symmetric matrix, not necessarily semi-definite positive. This means that the QCQP is in general a non-convex problem. Let δ be the vector of eigenvalues of ∆. By eigendecomposition, S = U T DU. where the matrix U is orthogonal. The entries of δ are the diagonal entries of D. The following theorem justifies the choice of penalization term in (6.4), by showing with an appropriate λ, adding this penalization term is equivalent to imposing the autocorrelation constraint (6.5).

Theorem 1. Suppose that δ 1 , the largest eigenvalue of S, is strictly positive. Suppose that z 0 ≡ Ux 0 has no zero component. Then there exists 0

≤ λ < 1 δ 1 , that verifies n 1 t=1 δ t z 2 0,t
2(1-λδt) 2 = 0, so that x * ≡ (I -λS) -1 x 0 is an optimal solution of (6.7).

Proof. We follow [START_REF] Ben | Hidden Conic Quadratic Representation of Some Nonconvex Quadratic Optimization Problems[END_REF] to obtain a convex relaxation of (6.7).

Define z ≡ Ux, z 0 ≡ Ux 0 , y t ≡ 1 2 z 2 t , ∀1 ≤ t ≤ n 1 . Recall that δ 1 > 0, and that ∀t, 1 ≤ t ≤ n 1 , z 0,t = 0. Problem (6.7) is equivalent to the non-convex problem min y,z

1 T y -z T 0 z s.t. -δ T y ≤ 0, 1 2 z 2 t = y t , ∀1 ≤ t ≤ n 1 . (6.8)
Now consider its convex relaxation min y,z

1 T y -z T 0 z s.t. -δ T y ≤ 0, 1 2 z 2 t -y t ≤ 0, ∀1 ≤ t ≤ n 1 .
(6.9)

By [START_REF] Ben | Hidden Conic Quadratic Representation of Some Nonconvex Quadratic Optimization Problems[END_REF]Theorem 3], if (z * , y * ) is an optimal solution of (6.9), and if 1 2 (z * t ) 2 = y * t , ∀1 ≤ t ≤ n 1 , then (z * , y * ) is also an optimal solution of (6.8), which makes x * = U T z * an optimal solution of (6.7).

We will look for such a solution to (6.9) by examining its first-order conditions of optimality. Problem (6.9) is convex, and it verifies the Slater condition: ∃(ŷ, ẑ), -δ T ŷ < 0, 1 2 ẑ2

t < ŷt , ∀1 ≤ t ≤ n 1 . This is true, because δ 1 > 0. We could choose an arbitrary value of ŷ1 > 0 and strictly positive but small values for other components of ŷ so as to have -δ T ŷ < 0, and ẑ = 0. Thus, Problem (6.9) always has an optimal solution, because the objective function is coercive over the constraint. This shows the existence of (z * , y * ).

Now we show that

1 2 (z * t ) 2 = y * t , ∀1 ≤ t ≤ n 1 .
The KKT conditions of (6.9) are verified by 6.2. RECONSTITUTION OF TIME SERIES WITH NMF (z * , y * ). In particular, there is some dual variable λ ≥ 0, µ ∈ R T + that verifies, 1 -λδ -µ = 0, (6.10) -δ T y * ≤ 0, (6.11)

λδ T y * = 0, (6.12)

-z 0,t + µ t z * t = 0, ∀1 ≤ t ≤ n 1 , (6.13) 1 2 (z * t ) 2 -y * t ≤ 0, ∀1 ≤ t ≤ n 1 , (6.14) µ t ( 1 2 (z * t ) 2 -y * t ) = 0, ∀1 ≤ t ≤ n 1 . (6.15)
Since z 0,t = 0, we have µ t = 0, z * t = 1 µt z 0,t , ∀1 ≤ t ≤ n 1 by (6.13). Therefore, by (6.15),

y * t = 1 2 (z * t ) 2 , ∀1 ≤ t ≤ n 1
. By (6.10), the values of µ t = 1 -λδ t can be deduced from that of λ. Since µ > 0, we obtain that λ < 1 δ 1 . By (6.13),

z * t = z 0,t 1-λδt , ∀1 ≤ t ≤ n 1 . This shows that x * = U T z * = U T (I -λD) -1 z * 0 = (I -λS) -1
x 0 is an optimal solution for (6.7).

Theorem 1 shows that with a well chosen λ, the constraint in (6.7) can be replaced by a penalization.

There are in fact two cases. Either x 0 verifies the constraint, in which case λ = 0, and x * = x 0 is the solution. Otherwise, λ > 0. We replug the values of

y * t = 1 2 (z * t ) 2 = z 2 0,t 2(1 -λδ t ) 2
back into (6.12), and obtain that λ verifies

T t=1 δ t z 2 0,t 2(1 -λδ t ) 2 = 0.
When F r and F c are fixed, the subproblem of (6.4) on V can be separated into N constrained problems of the form, min

x x -x 0 2 -λx T ∆ ρn x, s.t. A n x = c n ,
x ≥ 0, (6.16) where x 0 is the n-th column of F r F T c , c n is the observations on the n-th column, and A n is a matrix which encodes the measurement operator over that column. The following theorem shows how to solve this problem.

Theorem 2. Suppose that S is a symmetric matrix with eigenvalues δ, and λ

1 > 0. Suppose A ∈ R m,l a full-rank matrix with m ≤ l, x 0 ∈ R l , c ∈ R m , λ ≥ 0. Define Q ≡ (I -λS) -1 A T (A(I - λS) -1 A T ) -1 . If λ < 1 δ ρ,1 , then Qc + (I -QA)(I -λS) -1 x 0 is a minimizer of min x x -x 0 2 -λx T Sx, s.t. Ax = c, (6.17)

CHAPTER 6. NONNEGATIVE MATRIX FACTORIZATION FOR TIME SERIES RECOVERY FROM A FEW TEMPORAL AGGREGATES

Proof. Let l be the dimension of c. Define I C as the indicator function for the constraint of (6.17), that is

I C (x) = 0, if Ax = c,
and

I C (x) = +∞, if Ax = c.
Problem (6.17) is then equivalent to

min x F (x) ≡ 1 2 x -x 0 2 - 1 2 λx T Sx + I C (x). (6.18)
The subgradient of (6.18

) is ∂F (x) = {x -x 0 -λSx -A T | ∈ R l }. When λ < 1 δ 1 , (6.18
) is convex. Therefore, x * is a minimizer if and only if 0 ∈ ∂F (x), and

Ax * = c. That is, ∃ ∈ R l , (I -λS)x * -x 0 -A T = 0, Ax * = c.

The vector thereby verifies

A(I -λS) -1 (x 0 + A T ) = c.
The l-by-l matrix A(I -λS) -1 A T is invertible, because l is smaller than m, and A is of full rank (because each measurement covers disjoint periods). Therefore,

= (A(I -λS) -1 A T ) -1 (c -A(I -λS) -1 x 0 ), x * = (I -λS) -1 (x 0 + A T ) = Qc + (I -QA)(I -λS) -1 x 0 .
In our particular problem, the eigenvalues of ∆ ρn are

δ ρn,t = 2 cos( t n 2 + 1 π) -2ρ n ,
with t taking every value from 1 to n 2 . This means that for most of the autocorrelation threshold that we could need to impose (-1 ≤ ρ n ≤ 1), ∆ ρn has both strictly positive and strictly negative eigenvalues, allowing the above theorems to apply.

Both I -λ∆ ρn and A n (I -λ∆ ρn ) -1 A T n are invertible with λ < δ ρn,1 . The matrix inversion only needs to be done once for each individual. After computing

Q n ≡ (I -λ∆ ρn ) -1 A T n (A n (I - λ∆ ρn ) -1 A T n ) -1 , Q n c n and (I -Q n A n )(I -λ∆ ρn ) -1
for each n, we use Algorithm 2 to solve (6.4).

Choosing λ An optimal value of λ could be calculated. Substituting the values of y * in (6.12), shows that the optimal λ is a root of the polynomial n 2 t=1 δ ρ,t z 2 0,t 2(1-λδρ,t) 2 . The root-finding is too expensive to perform at every iteration. However, the optimal λ verifies 0 < λ < 1 δ ρ, 1 ,

where δ ρ,1 = 2 cos( 1 n 2 +1 π) -2ρ is the biggest eigenvalue of ∆ ρ . This gives us a good enough idea about how large a λ to use. In the numerical experiments, we chose λ = min(1, 1 2 maxn δ ρn,1 ) in the penalization when the constraint in (6.7) is active, and λ = 0 (no penalization) when the constraint is verified by x 0 .

EXPERIMENTAL RESULTS

Algorithm 2 Block coordinate descent for NMF from temporal aggregates and autocorrelation penalty

input ρ n , A n , Q n , Q n c n , ∀1 ≤ n ≤ N, and 1 ≤ K ≤ min{T, N } Initialize F 0 r , F 0 c ≥ 0, V 0 = P A (F 0 r F 0 c ), i = 0 while Stopping criterion is not satisfied do F i+1 r = Update(F i r , F i c , V i ) F i+1 c = Update(F i+1 r , F i c , V i ) for all 1 ≤ n ≤ N do v i+1 n = (Q n c n + (I -Q n A n )(I -λ∆ ρn ) -1 F i+1 r h i+1 n ) + end for i = i + 1 end while output V i ∈ A, F i r ∈ R n 1 ×k + , F i c ∈ R n 2 ×k + 6.
3 Experimental results

Datasets

We use one synthetic dataset and three real-world electricity consumption datasets to evaluate the proposed algorithms. In each dataset, the individual autocorrelation is calculated on historical data from the corresponding datasets, not used for evaluation.

• Synthetic data: 20 independent Gaussian processes with Matern covariance function (shifted to be nonnegative) are sampled over 150 periods to form the factor matrix F r . A 20-by-120 weight matrix F c is generated by sampling from a standard normal distribution truncated at 0, independently for each entry. The data matrix is obtained as F r × F c (n 1 = 150, n 2 = 120). This matrix is exactly of rank 20.

• French electricity consumption (proprietary dataset): daily consumption of 636 mediumvoltage feeders gathering each around 1,500 consumers based near Lyon in France during 2012 (n 1 = 365, n 2 = 636).

• Portuguese electricity consumption [28] daily consumption of 370 Portuguese clients during 2014 (n 1 = 365, n 2 = 370).

• Electricity consumption of small Irish companies [29,[START_REF]Results of Electricity Coast-Benefit Analysis, Customer Behaviour Trials and Technology Trials Commission for Energy Regulation[END_REF] daily consumption of 426 small Irish companies during 200 days in 2010 (n 2 = 200, n 2 = 426).

Evaluation procedure

For each individual in a dataset, we generate observations by selecting a number of observation periods. The temporal aggregates are the sum of the time series between two consecutive observation periods. The observation periods are chosen in two possible ways: periodically (at regular intervals with the first observation period sampled at random), or uniformly at random. The regular intervals for periodic observations are p ∈ {2, 3, 5, 7, 10, 15, 30}. This is motivated by the real application where meter readings are recorded regularly. With random observations, we use sampling rates that are equivalent to the regular intervals. That is, the number of observations D verifies D n 1 n 2 = 1 p ∈ {0.5, 0.33, 0.2, 0.14, 0.1, 0.07, 0.03}. We apply the following methods to recover the data matrix from each set of sampled observations:

• interpolation Temporal aggregates are distributed equally over the covered periods.

• softImpute As an alternative method, we apply a state-of-art matrix completion algorithm to complete the cumulative matrix. The observed entries are the cumulative values of the column from the first period to the observation dates. We use a nuclear-norm minimization algorithm, implemented in the R package, softImpute [111], to complete the cumulative consumption matrix, before differentiating each column to the obtain recovered matrix.

To choose the thresholding parameter, we use the warm start procedure documented in softImpute.

• HALS, and NeNMF These are the proposed matrix recovery algorithms using two classical F r and F c update implementations: HALS, and NeNMF. When autocorrelation penalization is used, we choose λ = min(1, 1 2 maxn δ ρn,1 ), as explained in the previous section. The rank used in proposed algorithms is chosen by a 5-fold cross validation procedure: we split the observations randomly into 5 folds, and apply the algorithm to 4 of the 5 folds with ranks 2 ≤ k ≤ 30. We then calculate the 2 -distance between the temporal aggregates on the recovered matrix with the 1-fold holdout. Repeating this procedure onto the 5 folds separately, we choose the rank which minimizes the average 2 -distance, to perform the algorithm on all observations. The cross validation procedure is carried out on the unpenalized version of the algorithms, and used both in the validation of unpenalized and penalized algorithms.

With each recovered matrix V obtained in an algorithm run, we compute the relative rootmean-squared error (RRMSE):

RRMSE(V, V * ) = V -V * F V * F .
Each experiment (dataset, sampling scheme, sampling rate, recovery method, unpenalized or penalized) is run three times to calculate the average RRMSE. Qualitatively similar results have been obtained using other error metrics.

Results

In Figure 6.1, one example of each data is shown. In the experiment depicted by the examples, 20% random observations are sampled. That is, observations dates are chosen uniformly at random, with one temporal aggregate every five days in average. In each dataset, the ground truth is drawn in black. The unpenalized NeNMF recovered the times series in green. We can see that generally trends are well estimated. However, unpenalized recovery produces an estimation which is more wiggly than the ground truth. With the autocorrelation-penalized version of NeNMF (orange lines), this aspect is improved.

The average RRMSE of experiments is reported in Figure 6.2. The figure is zoomed to show the RRMSE of the proposed algorithms. Much higher error rates for reference methods are sometimes not shown.

On sample sets with random observation periods (lower panel), proposed methods (HALS and NeNMF, blue and purple lines), whether unpenalized (solid lines) or penalized (dashed lines), out-performs the interpolation benchmark (red solid lines) by large in all datasets. This is especially the case when the sampling rate is small, i.e. when the task is more difficult. On the Irish dataset (lower panel, furthest to the right), penalized HALS and NeNMF (dashed blue and purple lines) are an improvement to unpenalized HALS and NeNMF when the sampling rate is low. With periodic observations (upper panel), the RRMSE is higher for every method. Proposed unpenalized methods, HALS and NeNMF (blue and purple solid lines) are equivalent to interpolation benchmark (red solid lines) for synthetic data, but sometimes worse for real datasets. Real electricity consumption has significant weekly periodicity, which is poorly captured by observations at similar periods. However, this shortcoming of the unpenalized method is more than compensated for by the penalization (dashed blue and purple lines).

We notice that penalized HALS and NeNMF consistently outperform interpolation with both observation schemes. This makes penalized methods particularly useful for the application of electricity consumption reconstitution, where it may be costly to install a random observation scheme, or to change the current periodic observation scheme to a random one.

It is also interesting to note that the rank chosen by the cross validation procedure is higher in higher sampling rate scenarios (Figure 6.3). This shows that the cross validation procedure is able to relax the rank constraint when more information is available in the data.

The traditional matrix completion method seems to fail in this application: softImpute (green solid lines) only has comparable results to interpolation or proposed methods in two of the four datasets, with 50% sampling rate in the random sampling scheme, which is the easiest case. In most cases, softImpute has an RRMSE much larger than 100%, and thus is not shown in the graphic. This indicates that the cumulative matrix considered in this application does not verify assumptions which guarantee matrix completion success.

Introduction

In recent years, a large number of methods have been developed to improve matrix completion methods using side information [START_REF] Jain | Provable Inductive Matrix Completion[END_REF][START_REF] Rao | Collaborative Filtering with Graph Information: Consistency and Scalable Methods[END_REF][START_REF] Si | Goal-Directed Inductive Matrix Completion[END_REF]. By including features linked to the row and/or columns of a matrix, these methods have a better performance at estimating the missing entries.

This chapter is a joint work with Yohann De Castro, Yannig Goude, Georges Hébrail and Jean-Marc Azaïs. It is based on the submitted article [3].

INTRODUCTION

Moreover, we suppose that the matrix of interest, V * , stems from a generative low-rank nonnegative model, in the following sense:

1. The matrix V * is of nonnegative rank k, with k n 1 , n 2 . This means, k is the smallest number so that we can find two nonnegative matrices

F r ∈ R n 1 ×k + and F c ∈ R n 2 ×k + satisfying V * = F r F T c .
Note that this implies that V * is of rank at most k, and therefore is of low rank.

2. There are d 1 row features X r ∈ R n 1 ×d 1 and d 2 column features X c ∈ R n 2 ×d 2 connected to each row and column of V * . We note by x i r the i-th row of X r , and by x i c the i-th row of X c . There are two link functions f r : R

d 1 → R k and f c : R d 2 → R k , so that F r = (f r (X r )) + , F c = (f c (X c )) + ,
where f r (X r ) ∈ R n 1 ×k is the matrix obtained by stacking row vectors f r (x i r ), for 1

≤ i ≤ n 1 (idem for f c (X c ) ∈ R n 2 ×k
), and (•) + is the ramp function which corresponds to thresholding operation at 0 for any matrix or vector.

In this general setting, the features X r and X c , the measurement operator A, and the measurements α are observed. The objective is to estimate the true matrix V * as well as the factor matrices F r and F c , by estimating the link functions f r and f c .

To obtain a solution to this matrix recovery problem, we minimize the quadratic error of the matrix factorization. In Section 7.3, we will propose an algorithm for the following optimization problem: min

V,fr∈F k r ,fc∈F k c V -(f r (X r )) + (f c (X c )) T + 2 F s.t. A(V) = b, V ≥ 0, (7.3)
where F r ⊆ (R) R d 1 and F c ⊆ (R) R d 2 are some functional spaces in which the row and column link functions are to be searched.

By specializing X r , X c , and A, or restricting the search space of f r and f c , this general model includes a number of interesting applications, old and new.

The masks A 1 , ..., A N

As mentioned in Chapter 5, varying the masks allows us to consider various applications in the matrix recovery context.

• Complete observation:

N = n 1 n 2 , A i 1 ,i 2 = e i 1 e T i 2
, where e i is the i-th canonical vector. This means every entry of V * is observed.

• Matrix completion: the set of masks is a subset of complete observation masks, with N < n 1 n 2 .

• Matrix sensing: the design matrices A i are random matrices, sampled from a certain probability distribution. Typically, the probability distribution needs to verify certain conditions, so that with a large probability, A verifies the Restricted Isometry Property [78].

INTRODUCTION

in the sense that for any solution (E r , E c ) to (7.4), ((E r ) + , (E c ) + ) is a solution to (7.5).

A more immediate generalization of (7.5) to include exogenous variables would be in the form min

V,fr∈F k r ,fc∈F k c V -f r (X r )(f c (X c )) T 2 F s.t. A(V) = b, V ≥ 0, f r (X r ) ≥ 0, f c (X c ) ≥ 0. (7.6)
Solving (7.6) would involve identifying the subset of F r and F c that only produce nonnegative value on the row and column features, which could be difficult.

Table 7.1: Classification of matrix factorization with side information by the mask, the link function, and the features included as side information, with some problems previously addressed in the literature.

Link function Linear Other regression methods

Features Identity General numeric features

Kernel features General numeric features

Mask

Identity Matrix factorization

Reducedregression rank [START_REF] Velu | Multivariate Reduced-Rank Regression: Theory and Applications[END_REF][START_REF] Bunea | Joint Variable and Rank Selection for Parsimonious Estimation of High-Dimensional Matrices[END_REF][START_REF] Chen | Sparse Reduced-Rank Regression for Simultaneous Dimension Reduction and Variable Selection[END_REF] Multiple kernel learning [START_REF] Kekatos | Electricity Market Forecasting via Low-Rank Multi-Kernel Learning[END_REF] Nonparametric RRR [START_REF] Foygel | Nonparametric Reduced Rank Regression[END_REF] Matrix completion Matrix completion [86] IMC [START_REF] Jain | Provable Inductive Matrix Completion[END_REF][START_REF] Agarwal | Regression-Based Latent Factor Models[END_REF][START_REF] Xu | Speedup Matrix Completion with Side Information: Application to Multi-Label Learning[END_REF][START_REF] Chiang | Matrix Completion with Noisy Side Information[END_REF] GRMF [START_REF] Rao | Collaborative Filtering with Graph Information: Consistency and Scalable Methods[END_REF][START_REF] Si | Goal-Directed Inductive Matrix Completion[END_REF][START_REF] Abernethy | A New Approach to Collaborative Filtering: Operator Estimation with Spectral Regularization[END_REF] Rank-one projections [84,85] Temporal aggregates [2] General masks Matrix recovery [78,105] By using (7.3), we actually shifted the search space F r and F c to (F r ) + and (F c ) + which consists of composing all functions of F r and F c with the ramp function (thresholding at 0). In a word, (7.3)

is equivalent to min V,fr∈(Fr) k + ,fc∈(Fc) k + V -f r (X r )(f c (X c )) T 2 F s.t. A(V) = b, V ≥ 0. (7.7)
Problem (7.7) also helps us to reason on the identifiability of (7.3). In a sense, (7.3) is not well-identified: two distinct elements in F k c × F k r have the same evaluation value of the objective function, if they only differ on their negative parts. In fact, this does not affect the interpretation of the model, because these distinct elements correspond to the same element in (F r ) k + × (F c ) k + . Since we are only going to use the positive parts of the function both in recovery and prediction, this becomes a parameterization choice which has no consequence on the applications.

As a comparison, we also propose an algorithm for the following optimization problem in Section 7.3: min

fr∈F k r ,fc∈F k c b -A((f r (X r )) + (f c (X c )) T + ) 2 2 , ( 7.8)
which corresponds to 5.7 mentioned in Chapter 5. Instead of minimizing the low-rank approximation error for a matrix that matches the data as a linear matrix equation in (7.3), (7.8) minimizes the sampling error of an exactly low-rank matrix. Both have been studied in the literature. For example, objective functions similar to (7.3) have been considered in [78], and ones similar to (7.8) in [79]. We will see in both Sections 7.3 and 7.4 that (7.3) has a better performance than (7.8).

where D n is a n-by-n diagonal matrix with only strictly positive coefficients on the diagonal and zeros everywhere else, and the (m -n)-by-n matrix M 0 is a collection of the other m -n rows of M.

Definition 2 (Strongly Boundary Closeness). A nonnegative matrix M ∈ R m×n + is said to be strongly boundary close if the following conditions are satisfied.

1. M is boundary close: for all i, j ∈ {1, ..., n}, i = j, there is a row m in M which satisfies m i = 0, m j > 0;

2. There is a permutation of {1, ..., n} such that for all i ∈ {1, ..., n -1}, there are n -i rows m 1 , ..., m n-i in M which satisfy (a)

m j i = 0, n s=i+1 m j i > 0 for all j ∈ {1, ..., n -i}; (b) the square matrix (m j s ) 1≤j≤n-i,i+1≤s≤n is of full rank (n -i).
Strongly boundary closeness demands, modulo a permutation in {1, ..., n}, that for each 1 ≤ i ≤ n -1, there are n -i rows m 1 , ..., m n-i of M that have the following form,

   m 1 : m n-i    T =          : . . . : 0 • • • 0 m 1 i+1 • • • m n-i i+1 : . . . : m 1 n • • • m n-i n          (i -1) first rows i-th row is all zero      (n -i) × (n -i) full rank matrix (7.9)
These row vectors, m 1 , ..., m n-i , all have 0 on the i-th element, and its lower square matrix of is of full rank. There are therefore enough linearly independent points on each n -1-dimensional facet R n + , which shows that cone(M T ) is somewhat maximal in R n + . The following was proved in [12]: Theorem 4. [12] If F r is strongly boundary close, then the only simplicial cone with

k generators in R k + containing cone(F T r ) is R k + . Moreover, if F c is separable, then V * = F r F T c
is the unique NMF of V * up to permutation and scaling.

Identifiability with side information

The NMF with linear row features, V * = X r B r F T c , is said to be unique, if for all matrix pairs ( Br , Fc ) ∈ R d 1 ×k × R n 2 ×k that verifies X r Br ≥ 0, Fc ≥ 0, V * = X r Br Fc , we have B = B r , Fc = F c up to permutation of columns and scaling.

|The columns of M are strongly boundary close};

F (X) = {M ∈ R n 1 ×k + |rank(M) = k, span(M) ∈ span(X)}. Theorem 5. If E F (X r ) = ∅, and B r ∈ (X T r X r ) -1 X T r (E F (X r )), and
F c is separable, then the factorization V * = X r B r F T c is unique. Proof. Notice that for B r ∈ (X T r X r ) -1 X T r (E F (X r
)), the nonnegative matrix X r B r is strongly boundary close. The factorization (X r B r , F c ) is therefore unique. The model identifiability follows immediately, since X r is of full rank.

Example of X r that verifies E F (X r ) = ∅ For this theorem to have practical consequences, one needs to find appropriate row features so that E F (X r ) = ∅.

Here we provide a family of matrices X r so that E F (X r ) = ∅.

With a fixed k ≥ 2, suppose that X r has k(k -1)/2 columns, and at least k(k -1)/2 rows, with the first k(k -1)/2 + 1 rows defined as the following:

• the first row and column have 0 on the first entry and positive entries elsewhere;

• for 2 ≤ i ≤ k, X r has strictly positive entries on the first ((i -1)(i -2)/2 + 1) columns, from Row (i -1)(i -2)/2 + 3 to Row (i -1)(i -2)/2 + 1 + i, and zero entries everywhere else. These (k -1) rows are linearly independent.

Then we have E F (X r ) = ∅, because the following k(k -1)/2-by-k matrix B r is in this set:

• for 1 ≤ i ≤ k, B * K has i consecutive strictly positive entries on the i-th column, between Row i(i -1)/2 + 1 and Row i(i -1)/2 + i.

The following matrices instantiate the case of k = 4:

B r =             1 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 1             , X r =              
0 5 14 7 9 15 13 10 0 0 0 0 0 0 4 5 0 0 0 0 0 12 4 0 0 0 0 0 10 7 10 7 0 0 0 13 10 12 9 0 0 0 12 10 16 8 0 0 0 : : : : : : :

              ,F r =              
0 5 21 37 10 0 0 0 4 5 0 0 12 4 0 0 10 7 17 0 13 10 21 0 12 10 24 0 : : : :

              . 7.3. HALSX ALGORITHM component on X if for every x ∈ X and y i ∈ X i , we have g(x 1 , x 2 , ...tx i + (1 -t)y i , ..., x m ) ≤ max{g(x), g(x 1 , x 2 , ..., y i , ...x m )} for all t ∈ [0, 1]
. g is said to be strictly quasi-convex with respect to the i-th component, if with the additional assumption that y i = x i , we have

g(x 1 , x 2 , ...tx i + (1 -t)y i , ..., x m ) < max{g(x), g(x 1 , x 2 , ..., y i , ...x m )} for all t ∈]0, 1[.
It has been shown that if g is strictly quasi-convex with respect to the first m -2 blocks of components on X, then a limiting point produced by a Gauss-Seidel algorithm is a critical point [110]. This result is not directly applicable for the HALS algorithm. Typically, if f c,i , the i -th column of F c , is identically zero, the loss function is completely flat respect to f c,i , the i-th column of F r . Therefore the loss function is not strictly quasi-convex. In order to avoid this scenario, [18] suggests thresholding at a small positive number instead of at 0, when updating each column of the factor matrices.

In fact the convexity assumption of [110] can be slightly relaxed to directly apply to HALS, as demonstrated by the following proposition. Theorem 6. Suppose that the function g is quasi-convex with respect to x i on X, for i = 1, ..., m -2. Suppose that some limit points x of the sequence {x t } (t∈N) verify that g is strictly quasi-convex with respect to x i on the product set {x 1 } × {x 2 } × ... × X i ×, ... × {x m }, for i = 1, ..., m -2. Then every such limiting point is a critical point of Problem (3).

Compared to the result of [110], this shows that the strict convexity with respect to one block does not have to hold universally for feasible regions of other blocks. It only needs to hold at the limiting point. This theorem can be established following the proof of Proposition 5 of [110], using the following lemma. Lemma 1. Suppose that the function g is quasi-convex with respect to x i on X, for some i ∈ {1, ..., m}. Suppose that some limit points ȳ of {y t } verify that g is strictly quasi-convex with respect to x i on {ȳ 1 } × {ȳ 2 } × ... × X i ×, ... × {ȳ m }. Let {v t } be a sequence of vectors defined as follows:

v t j = y t j if j = i, arg min z i ∈X i g(y t 1 , ..., z i , ..., y t m ) if j = i.
Then, if lim t→+∞ g(y t )-g(v t ) = 0, we have

lim t→+∞ ||v t i -y t i || = 0. That is lim t→+∞ ||v t -y t || = 0.
Proof. (The proof of the lemma is based on [START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF].) Suppose on the contrary that ||v t i -y t i || does not converge to 0. Define τ k = ||v t i -y t i ||. Restricting to a subsequence, we can obtain that τ k ≥ τ 0 > 0. Define s t = v t -y t τ k . Notice that {s t } is of unit norm, and v t = y t + τ k s t . Since {s t } is on the unit sphere, it has a converging subsequence. By restricting to a subsequence again, we could suppose that {s t } converges to s.

Since

θ * ∈ arg min θ∈R d R -g(θ)(f c ) T 2 F ,
is a minimizer of a smooth convex problem,

∂ ∂θ R -g(θ)(f c ) T 2 F (θ * ) = ∇g θ * (R -g(θ * )(f c ) T )f c = 0. This means (R -g(θ * )(f c ) T )f c = 0, because ∇g θ * is of full rank. Consequently g(θ * ) = 1 f c 2 2 Rf c .
It has been shown in NMF literature (for example [18,Theorem 2]) that, (g(θ * )) + = arg min

f ∈R n 1 + R -f (f c ) T 2 F .
This is equivalent to

g(θ * ) ∈ arg min f ∈R n 1 L(f ), or 0 ∈ ∂L g(θ * ) .
We therefore conclude with 0 ∈ ∇g θ * ∂L g(θ * ) .

In many regression methods, even when a non-linear transformation is applied to the data, the regression function is linear in its parameters. A non-exhaustive list of methods include linear regression (g(θ) = X r θ), spline regression (g(θ) = φ(X r )θ), or support vector regression (SVR) (g(θ) = K(X r , X r )θ). In this case, g has a constant Jacobian matrix. In the case of linear and spline regression, the Jacobian matrix is of rank n 1 if there are no less features than examples. For SVR, this is true for any positive definite kernels. This allows us to apply the previous lemma to each column update step of Algorithm 4. By calculating

f t+1 r,i = arg min f ∈Fr R t - f (X r )(f t c,i ) T 2 F at
Step t for Column i in F r , we actually have

f t+1 r,i = arg min f ∈Fr R t -(f (X r )) + (f t c,i ) T 2 F .
This shows that at each iteration, we solve the subproblems of (7.3).

In these cases, by rewriting the functional space F r and F c in a parametric form, the search space is actually R r 1 and R r 2 , for some r 1 and r 2 .

Theorem 7. If n 1 ≤ r 1 , n 2 ≤ r 2 ,

every full-rank factorization produced by HALSX (Algorithm 4) is a critical point of Problem (7.3).

Remark 1. In order to obtain a similar result for Problem (7.6), one would need to ensure the obtained functions have non-negative values on the features. This could be done by alternating projection.

Designs and HALSX

At each iteration of Algorithm 4, we need to project the working matrix F t r (F t c ) T into the convex polytope defined by the measurements and nonnegativity:

V t = arg min V|A(V)=b,V≥0 V -F t r (F t c ) T 2 F . (7.12)
In general, the polytope projection can be obtained by alternating projection. Namely, we can alternate between:

• V = V + A † (b -A(V)); • v i,j = max(0, v i,j ),
where A † is the right pseudo-inverse of A, viewed as an N -by-n 1 n 2 matrix.

For some measurement operators, there are efficient ways to solve (7.12).

• Matrix completion mask:

v i,j = α l , if ∃1 ≤ l ≤ N, A l = e i e T j ; max(0, v i,j ), if not.
• Temporal aggregate mask: simplex projection (see [2] for details).

Linear HALSX

In this section, we consider HALSX with numeric row features and linear row link functions. That is, given X r and α = A(V * ), we need to solve min

V∈R n 1 ×n 2 ,Br∈R d 1 ×k ,Fc∈R n 2 ×k V -(X r B r ) + (F c ) T + 2 F s.t. A(V) = b, V ≥ 0, (7.13)
Following Algorithm 4, we need to update the columns of B r at each iteration. At the t-th step, for 1 ≤ i ≤ k, we solve the subproblem

arg min b r,i R t i -X r b r,i (f t c,i ) T 2 F ,
where

R t i = V t -k j=1,j =i X r b r,l f T c,l .
This minimization problem has a closed-form solution:

b t+1 r,i = 1 f t c,i 2 2 (X T r X r ) -1 X T r R t i f t c,i .
In order to accelerate the numerical algorithm, a QR decomposition of X r = QR is done before the iterations, where Q is an orthogonal matrix, and R is a square upper triangular matrix. When X r is of full rank, X T r X r is invertible. We compute one time (X T r X r ) -1 X r = R -1 Q T before the iterations, and use the result at each iteration.

HALSX ALGORITHM

Noting by X j = {a j l (x i )} i,l the design matrix, for h j = L j l=1 β j l a j l ∈ H(a j , L j ), an element of the functional space, we have

h j (X) = X j β j .
The whole model of g, can then be represented linearly:

g(µ) = Xθ + (X 1 , X 2 , ...)    β 1 β 2 :    = Xθ + j X j β j .
The dimension of H(a j , L j ), L j , controls the smoothness of the functions to be estimated. As little information is available on the degree of smoothness of the functions, we use a rather large L j , and add a penalty on the wiggliness, (h j ) 2 dx, as in [START_REF] Wood | Generalized Additive Models: An Introduction with R[END_REF]. The least squares estimator of this model is therefore

arg min θ,β 1 ,β 2 ,... g(µ) -Xθ - j X j β j 2 + j λ j (β j ) T S j β j ,
where λ j is the penalization parameter of the j-th non-linear component, and S j is a positive definite matrix depending on X and a j . The penalization parameter, λ j , is chosen by a generalized cross validation criterion.

HALSX-GAM

At each iteration of the algorithm, for i = 1, ..., k, we re-estimate the link function f r,i of the i-th column of F r as a GAM.

The subproblem for i is the following arg min

θ,β 1 ,β 2 ,... R t i -(X r θ + j=1 X j β j )(f t c,i ) T 2 F + j λ j (β j ) T S j β j . ( 7.14)
With fixed penalization parameters λ 1 , λ 2 , ..., the optimization above can be solved by

     θ β 1 β 2 :      t+1 = 1 f t c,i 2 ×     X T r X r X T r X 1 X T r X 2 • • • (X 1 ) T X r (X 1 ) T X 1 + λ 1 f t c,i 2 S j (X 1 ) T X 2 • • • : : . . . • • •     -1 ×      X T r (X 1 ) T (X 2 ) T :      R t i f t c,i .
In practice, we use the GAM estimation routines implemented in the R package mgcv [START_REF] Wood | Generalized Additive Models: An Introduction with R[END_REF] to choose the penalization parameter and estimate the model at the same time.

(A i f c )(A i f c ) T ) -1 N i=1 b i A i f c
This computation generalizes to linear exogenous variables, with the optimality condition:

β = ( N i=1 (XA i f c )(XA i f c ) T ) -1 N i=1 b i A i f c .
When the matrices to be inversed in these equations are not invertible, we will use the generalized inverse instead.

Using these elementary steps, we propose Algorithm 5 to solve Problem (7.8). Compared to Algorithm 4, Algorithm 5

• has one less block (the slack variable V is not present);

• checks the deviation with data more frequently;

• each subproblem is more costly because of the presence of A in the subproblem. As we will see in the detailed development of the computation, when N , the sample size (the dimension of image of A) is large, each update involves rather costly computations.

Algorithm 5 Hierarchical Alternating Least Squares with eXogeneous variables for NMF (HALSX2)

Require: Measurement operator A, measurements b, rank 1 ≤ k ≤ min{n 1 , n 2 }, features X r and X c , functional spaces F r and F c in which to search the link functions.

Initialize F 0 r , F 0 c ≥ 0, t = 0 while Stopping criterion is not satisfied do R t = b -A(F t r (F t c ) T ) for i = 1, 2, ..., k do R t = R t + A(f t r,i (f t c,i ) T ) f t+1 r,i = arg min f ∈Fr R t -A(f (X r )(f t c,i ) T ) 2 2 f t+1 r,i = max(0, f t+1 r,i (X r )) R t = R t -A(f t+1 r,i (f t c,i ) T ) end for for i = 1, 2, ..., k do R t = R t + A(f t+1 r,i (f t c,i ) T ) f t+1 c,i = arg min f ∈Fc R t -A(f t+1 r,i f (X c ) T ) 2 2 f t+1 c,i = max(0, f t+1 c,i (X c )) R t = R t -A(f t+1 r,i (f t+1 c,i ) T ) end for t = t + 1 end while output F t r ∈ R n 1 ×k + , f t r,1 , ..., f t r,k ∈ F r , F t c ∈ R n 2 ×k + , f t c,1 , ..., f t c,k ∈ F c .
• HALS, NeNMF, softImpute Matrix recovery/completion methods without side information. See the experiments of the previous chapter for more details.

The following methods use side information, and can be used for time series prediction of new columns and/or rows from incomplete data:

• individual_gam Estimating separate GAMs on each column or row, on the matrix obtained from empty_model or on the whole data matrix when it is available.

• factor_gam Estimating GAMs on the factors obtained by HALS or NeNMF.

• rrr [START_REF] Addy | Reduced-Rank Regression [R Package Rrr Version 1.0.0[END_REF] Applying reduced-rank regression on the matrix obtained from empty_model or on the whole data matrix when it is available.

• grmf [START_REF] Rao | Collaborative Filtering with Graph Information: Consistency and Scalable Methods[END_REF] A matrix completion algorithm using graph-based side information to enhance collaborative filtering performance.

• trmf [26] A matrix completion algorithm tailored to time series, by adding three penalization terms to the matrix factorization quadratic error. When only temporal aggregate measurements are available, we apply this method on the matrix obtained from empty_model.

• HALSX_model Algorithm 4.

For all matrix factorization methods (HALS, NeNMF, rrr, factor_gam, grmf, trmf, HALSX_model), we use the method with several ranks, then choose the best rank (k ∈ {2, 3, ..., 20} for synthetic and Portuguese data, k ∈ {2, 3, ..., 10} for French and MovieLens data). For trmf, we do a grid search on the three penalization parameters, and choose the best combination.

As explained in the previous section, a regression method needs to be specified in order to use HALSX_model. This regression method specifies the fonctional spaces F r and F c in which the link functions are to be looked for. In our experiments, we use four different regression methods with HALSX: linear model, GAM, support vector regression with linear kernel, and Gaussian process regression with radial basis kernel (lm, gam, svmLinear, gaussprRadial). We use the implementation of lm in standard R, and svmLinear and gaussprRadial of the R package kernlab [START_REF] Karatzoglou | Kernlab-an S4 Package for Kernel Methods in R[END_REF], through the caret API. As of gam, we use the mgcv implementation directly.

For each data matrix, we apply a linear measurement operator on an upper-left submatrix to obtain the measures (of dimension 100-by-130 for synthetic data, 730-by-270 for French data, 731-by-369 for Portuguese data, and 666-by-1189 for MovieLens data). A random shuffle between columns is done before dividing the matrix into submatrices. As in Chapter 6, on time series datasets (the first three), both periodic and random temporal aggregates are generated as observations. On the MovieLens dataset, observations are random entries, as in matrix completion. Using the measurements of the upper-left submatrix, we use each method to estimate a matrix factorization model, with or without side information.

We then report the matrix recovery error on this sampled submatrix for all methods, and the prediction errors on the rest of the data matrix for the methods can produce predictions for new columns and/or rows. We distinguish row prediction error, column prediction error, and rowcolumn prediction error, as the error calculated on the lower-left, upper-right and lower-right submatrix (see Figure 7.1).

We report the relative root-mean-squared error as error metric for time series datasets: show the prediction error on the three datasets. We can see that trmf [26] and rrr [START_REF] Addy | Reduced-Rank Regression [R Package Rrr Version 1.0.0[END_REF] are not well adapted to temporal aggregate measurements. When they are applicable (trmf is only applicable to row prediction, rrr only applicable to row or column predictions, not RC predictions), they have much worse performance than the other rate cases, the NMF methods without side information (HALS and NeNMF) work the best, with HALSX_model with lm, gam or svmLinear are second to best. HALSX_model with gaussprRadial does not work very in this case, as is the case with the other comparison methods (grmf, softImpute and empty_model) Figure 7.10 shows the prediction error on MovieLens data for new users and/or new movies. The order of the variants of the HALSX_models is conserved: gam and lm are the best, svmLinear is not very good for 10%, but better with higher sampling rates, and gaussprRadial does not work well in this problem. Otherwise, the factor_gam method is slightly better for column predictions (new movies) in higher sampling rate cases, but worse in other cases. Both individual_gam and rrr are much worse than the proposed methods.

• for electricity datasets, RRMSE(V, V * ) = V-V * F V * F ,

Conclusion

Motivated by electricity consumption estimation, we proposed a general approach for including side information on the columns and row in nonnegative matrix factorization methods, with general linear measurements. Based on a generative model, the framework we propose generalizes many prior works in multivariate regression and matrix factorization.

In order to explore the identifiability of the model, we established a sufficient condition on the features for the factorization to be unique. We deduced HALSX, a general algorithm to estimate the generative model, and showed that the algorithm converges to a critical point in rather mild conditions.

The proposed algorithm is tested on synthetic and real datasets, both in electricity consumption and recommendation systems. In various sampling scenarios, HALSX produced better or equivalent performance both in matrix recovery and in prediction, compared to a number of reference methods.

Consider the simple case where n 1 = 1, and the measurements are simply the observations at the first n 2 -1 locations. Since the matrix to be estimated is of dimension 1-by-n 2 , the rank k can only be 1. We set the column factor f c (X c ) to be 1, since it is just a real number. We are then in the exact same setting as univariate kriging considered in Chapter 3. That is, knowing the value of V * at locations x c,1 , x c,2 , ..., x c,n 2 -1 : v * 1 , v * 2 ..., v * n 2 -1 , we'd like to estimate the its value at position x c,n 2 : v * n 2 . In simple kriging, the idea is to suppose that the underlying Gaussian process is second-order stationary. By minimizing the quadratic risk, we obtain that the optimal prediction is

v n 2 = (C(x c,1 -x c,n 2 ), ..., C(x c,n 2 -1 -x c,n 2 ))(C(x c,i -x c,j ) 1≤i,j≤n 2 -1 ) -1 (v * 1 , v * 2 ..., v * n 2 -1 ) T ,
where C is the covariance function of the Gaussian process.

Another way to view the kriging result is the following: we transform the spatial coordinates with a feature map φ, so that

φ(x) = (C(x c,1 -x), C(x c,2 -x), ..., C(x c,n 2 -1 -x)) T ,
and we restrict the link function f c to be a linear function of the transformed features. This way, the matrix factorization problem (8.1) becomes min

w∈R n 2 -1 ,V∈R 1×n 2 V -(φ(X c )w) T 2 F , s.t. v 1 = v * 1 , ..., v n 2 -1 = v * n 2 -1 .
We can easily obtain that the optimal value of w is

(C(x c,i -x c,j ) 1≤i,j≤n 2 -1 ) -1 (v * 1 , v * 2 ..., v * n 2 -1
) T . This yields exactly the same optimal prediction value for v n 2 .

Matrix factorization and sum-product covariance functions

Here we try to recover a matrix V * with more than one row and column: both n 1 and n 2 are potentially larger than 2. The covariance structure of spatio-temporal processes is an object that is more complex, as discussed in Chapter 3. A spatio-temporal covariance function C can be evaluated on any spatio-temporal increment (s, τ ). Typically, we would evaluate C on (s, τ ) where s is the (element-wise) difference between two rows in X c , and τ the difference between two rows in X r . We start by drawing the analogy between a separable spatio-temporal process with a rank-1 matrix factorization model. We use the subscript index r and c to designate temporal and spatial components, as rows correspond to periods and columns correspond to spatially located individuals in our notation.

A separable spatio-temporal covariance function is one that verifies C(s, τ ) = C r (τ )C c (s). One spatio-temporal process that verifies such a covariance function is the product of a spatial Gaussian process Z c and a temporal Gaussian process Z s . In this case, it is clear that the matrix that we wish to recover, V * , has a rank-1 factorization f r f T c , where f r and f r are simply the observations of the temporal and spatial processes at the periods and the locations in the dataset. Obviously, a rank-1 matrix has proportional rows and columns. Therefore, if one single row and one column are entirely observed, the matrix completion or kriging problem becomes trivial in this case.

As mentioned in Chapter 3, a way to obtain slightly more general covariance structure is to use sum-product covariance. Here, consider a sum-product covariance with k independent separable components:

C(s, τ ) = k i=1 C i r (τ )C i c (s).

MATRIX FACTORIZATION AND SOCIO-DEMOGRAPHIC CLUSTERING

In a similar fashion as in the univariate case, we use the spatial and temporal covariance functions to create feature maps φ i r and φ i c , where

φ i c (x c ) = (C c (x c,1 -x c ), C c (x c,2 -x c ), ..., C c (x c,n 2 -x c )) T ,
for all spatial coordinates x c , and

φ i r (x r ) = (C r (x r,1 -x r ), C r (x r,2 -x r ), ..., C r (x r,n 1 -x r )) T ,
for all temporal coordinates x r . This leads to the following matrix factorization problem:

min

Wr∈R n 1 ×k ,Wc∈R n 2 ×k ,V∈R n 1 ×n 2 V - k i=1 φ i r (X r )w i r (φ i c (X c )w i c ) T 2 F , s.t. A(V) = α.
In practice, the idea would be not to specify the feature maps (hence the covariance components) by hand, but let Algorithm 4 find appropriate ones. This can be done by specifying for example a class of kernel functions, and let the regression procedure find the window sizes by cross validation for each rank.

Matrix factorization and socio-demographic clustering

In Chapter 4, we proposed a hybrid method which has two components:

• a clustering component, which groups individuals into clusters, based on socio-demographic characteristics;

• a regression component, which proposes an estimation of the electricity consumption of all individuals of a cluster, with temporal covariates such as the temperature, the hour of the day, and the time of year.

These two components can easily be transformed into the nonnegative matrix factorization framework V * = f r (X r )f c (X c ), where • X r is the temporal covariates, f r is the regression models of consumption, in the same way as in Chapter 7;

• X c is the socio-demographic characteristics, and f c is the classifier: that is, for an individual with socio-demographic features x c , f c (x c ) = e i if this individual is in the i-th cluster, where e i is the i-th canonical vector in R k . In the framework of Algorithm 4, this can be achieved by estimating f c as a multinomial logistic regression.

The method proposed in Chapter 4 is thus a special case of Algorithm 4, where we only do one iteration.

Matrix factorization as feature generation

It is interesting to note that the analogy between matrix factorization and some of the more recent feature generation methods in machine learning, in particular in the neural network literature.

Matrix factorization is an extremely simple neural network: a depth-two auto-encoder with linear activation. One of the reasons why people are interested in convergence results of matrix factorization is because these analyses can shed light on the theoretical understanding of deeper neural networks [90].

Recent trends in machine learning showed that automatic feature generation methods, such as word2vec in natural language processing (which is also a depth-two neural network) and variational auto-encoders in image processing, can have very good empirical performance in applications. Similarly, in this chapter, we saw that the two more hand-tuned methods proposed in Chapters 3 and 4 can actually be written as special cases of the general approach proposed in Chapter 7. Implementation of nonnegative matrix factorization algorithms During the thesis, an R package called meterModels has been developped to implement the algorithms proposed in the previous two chapters. In this chapter, we describe the main functionalities of meterModels. When it is necessary, we will also provide some simple code to illustrate how to use the package.

The meterModels package includes several matrix factorization algorithms, as well as many helper functions to perform experiments comparing different methods. Most of the package is written in R. Heavy computations are translated in C++ to increase the performance of the package, through the RCpp API. When available, openMP can be used to perform distributed computation in model estimation. On a multi-core machine, parallel estimation of multiple models can be used through the parallel package in R.

On a single machine, the package can be used on matrices with several hundred columns and 100,000 rows in less than ten minutes.

Data representation

Most matrix factorization software known to us can only deal with observations in the form of matrix entries, thus only adapted to matrix completion. Given the applications considered in this thesis, we developed meterModels to handle more general measurements. To achieve this, meterModels provides several data representation formats, for temporal aggregate data (Section 6) and general linear measurements. We describe these formats in this section.

This chapter is based on the vignette of meterModels.

CHAPTER 9. IMPLEMENTATION OF NONNEGATIVE MATRIX FACTORIZATION ALGORITHMS

incremental_data: A fully-observed data matrix

As indicated by the name, the object class incremental_data is particularly motivated by electricity metering. Basically, incremental_data is a full data matrix (in the field incremental_data$data), with a string in incremental_data$data_id as the name of the dataset, and optionally exogenous variables provided through incremental_data$exo_variable. # create an incremental_data object incre_d = incremental_data(x = example_matrix, data_id = "exponential")

index_data: Temporal aggregates

index_data is to be used to provide temporal aggregate data to a model estimation algorithm. The class index_data is thus named because a meter reading is sometimes called an "index".

To create an object of class index_data, there are two ways:

• We can provide an incremental_data object, and an object of the class sampling_scheme to tell the program how to sample the data matrix. For example, in the following code, we initiated a sampling scheme, sampling_incre_d_at which randomly takes 50% of the temporal aggregates:

seed = 42
# we will sample incre_d with a random sampling_scheme # with a 50% sampling rate and a fixed seed sample_incre_d_at = sampling_scheme( sampling_rate = 0.5, type = "random", seed = seed)

# create an index_data object index_d = index_data(x = sample_incre_d_at, incremental_data = incre_d) lapply(index_d,head) • Or we can provide the dimension of the matrix to be estimated, the values of readings, and the position of the readings in the aggregated version of the corresponding incremental_data object. In the following code, we create the same sample as above, using directly the samples in index_d. index_d2 = index_data(x = index_d$sample_index, sample_data = index_d$sample_data, m = index_d$dim_data [1], n = index_d$dim_data [2], data_id = index_d$data_id) lapply(index_d2,head)

## $sample_data ## [
## $sample_data ## [
# ## $data_id ## [1] "exponential" ## ## $dim_data ## [1] 10 15
The created index_data object has four fields:

• data_id for the ID of the data as in incremental_data,

• dim_data for the dimension of the underlying data matrix, • and sample_index and sample_data that indicate the index of the measurements and the readings in the cumulative matrix, obtained by calculating the cumulative sums of each column of the incremental_data matrix. This is the primary data format used in Chapter 6. When this data format is used, most of the new algorithms in meterModels fully take advantage of it to estimate models efficiently.

linear_measurement_data: General linear measurements

A more general data format, linear_measurement_data, is available, for cases where the observations are linear measurements, but not necessarily temporal aggregates. The masks are provided through a sparseMatrix object (from a widely used R package called Matrix). To create a linear_measurement_data object, we provide the sparse matrix which tells the program how to measure the data, and an incremental_data object which is the matrix to measure.

exp_data

exp_data is a class to be used when running many algorithm runs, in an experiment (see below). It is created from an object of the class incremental_data, with some additional information as to which part of the data matrix to be used as training data and test data.

For example, the following code wraps incre_d into an object on which we will estimate models using the whole matrix as train data (the parameters m_train and n_train control this behavior).

# define the experiment data exp_d = exp_data(x = incre_d, m_train = m, n_train = n) 9.2 matrix_model: Estimate a single model meterModels has a number of matrix factorization algorithms included. Here we will use one of the proposed algorithms, NeNMF, which is an extension of [20].

To estimate a model, call the function matrix_model with the data object and the name of the algorithm, and a fixed rank. To see the estimated matrix, use predict on the estimated model. nenmf_model = matrix_model(algorithm = "NeNMF", data = index_d, r = 2) nenmf_pred = predict(nenmf_model, completion = TRUE)

We can plot the original and estimated model to see how well the matrix is recovered. In Figure 9.1, we can see that most of the original matrix is well estimated. There is some problem on the first column, probably due to the fact that only two temporal aggregates are observed on this column.

plot_matrices(list(

original = example_matrix, nenmf_pred = nenmf_pred))

We can also call matrix_model with incremental_data objects, so that only to factorize the matrix. In this case, the estimated matrix is exactly the same as the original (Figure 9.2). In order to estimate an NMF with exogenous variables as side information (see Chapter 7), we use HALS_model as algorithm, and provide a formula_params object to specify which regression models to estimate: regression_model_full_data = matrix_model( algorithm = "HALS_model", data = incre_d, r = 2, exo_variable = list(X = data.frame(variable = rnorm(dim(incre_d) [1]))), formula = list(X = formula_params( "s(variable)", # this is the right part of the formula HALS_model_params = list(reg_method = "gam")))) Reference method Table 9.1 lists the algorithms included in meterModels. Among the reference methods, softImpute, grmf, trmf, and reduced_rank_regression are included by directly using the open-source programs of the respective authors. The other reference methods are in fact developed by us, although they are not the main propositions of this thesis. All of the proposed methods can be used jointly with the autocorrelation constraint, described in Chapter 6, even the ones proposed in Chapter 7.

Model evaluation

To evaluate obtained matrix factorization model, we can calculate the relative difference in a matrix norm (for example Frobenius norm), by using relative_diff_norm.

relative_diff_norm(matrix1 = incre_d$data, matrix2 = nenmf_pred, type = "F") ## [1] We can also check that the recovered matrix has almost the same measurements on the mask, by using evaluate_on_mask.

evaluate_on_mask(increData = incre_d$data, increEst = nenmf_pred, known = index_d$sample_index)

## [1] 3.230599e-17

Run an experiment to estimate several models simultaneously

Many parameters are available for training a matrix_model. Therefore, we often need to try a number of combinations of the parameters, to find the best one. This is why we have the function experiment. To use experiment, we need to first define an exp_data object, and an experiment specification object (exp_spec).

# define the experiment specs simple_spec = exp_spec( # try several ranks r = 2:5, # use NeNMF, HALS, and empty_model for constant interpolation algorithm = c("NeNMF", "HALS", "empty_model"), # save the models in a temporary directory save_model = "./exp_temp/", # use 3 parallel processes parallel_opts = list(mc.cores = 3), # we will only do nonnegative matrix factorization general_factorization = FALSE)

The estimated models are saved on disk at the path provided through save_model parameter in simple_spec. Here we are using 3 cores to run model training in parallel (parallel_ops). We are only going to do nonnegative matrix factorization (general_factorization = FALSE).

We launch the experiment by calling experiment, with the experimental data, the sampling scheme, and the experimentation specification.

# then we are good to go simple_exp = experiment(exp_data = exp_d, sampling_scheme = sample_incre_d_at, exp_spec = simple_spec) The experiment function returns a table with the evaluation metrics of the estimated models (Table 9.2). It looks like the HALS algorithm with rank fixed at 2 obtains the best model (smallest completion error). This is not surprising, since the data is a simulated matrix exactly of rank 2.

Using Table 9.2, we can generate plots similar to those reported in the last two chapters (Figure 9.3).

Distributed computation in model estimation

The NMF algorithms introduced in this thesis are intrinsically iterative1 . Given the format of the observations (either temporal aggregates or general linear measurements), distributed stochastic gradient descent algorithms proposed in matrix completion and recommendation systems are not easily adaptable to the problem [START_REF] Recht | Parallel Stochastic Gradient Algorithms for Large-Scale Matrix Completion[END_REF][START_REF] Yun | NOMAD: Non-Locking, stOchastic Multi-Machine Algorithm for Asynchronous and Decentralized Matrix Completion[END_REF][START_REF] Chin | A Fast Parallel Stochastic Gradient Method for Matrix Factorization in Shared Memory Systems[END_REF][START_REF] Chin | LIBMF: A Library for Parallel Matrix Factorization in Shared-Memory Systems[END_REF].

The parallel solution that we finally developped is based on the CCD++ algorithm [START_REF] Yu | Scalable Coordinate Descent Approaches to Parallel Matrix Factorization for Recommender Systems[END_REF]. The CCD++ algorithm is mostly a parallel HALS, which is discussed in previous chapters. When updating the columns of F r and F c , this algorithm updates the coefficients in these columns in parallel. In order to combine this algorithm with the simplex projection and autocorrelation penalization discussed in Chapter 6, this algorithm was re-implemented in Rcpp. Using openMP, both the update of factors and the projection into the simplex satisfying the linear measurement In meterModels, this parallelization can be used agnostically (not available with regression models). To choose to use this functionality, the user can use the opts parameter in matrix_model.

Part III

In this thesis, we studied methods to model multivariate electricity consumption time series from partial data. We considered the problem from several aspects, with either temporally aggregated or individually aggregated data, either to estimate past consumption or to predict future consumption. We studied solutions from three contexts: spatio-temporal statistics, clustering of individuals using socio-demographic information, and NMF-based methods. In NMF, we extended the literature to have solutions that are adapted to our problem, and proposed improvements to better achieve the estimation goals. Particularly, we proposed to take into account additional information such as temporal autocorrelation and exogenous variables in NMF, and deduced algorithms to solve these extended models. The proposed methods were tested on synthetic and real datasets, and showed promising results.

Perspective applications in electricity

In this thesis, empirical studies were carried out to evaluate the estimation and prediction performance of the matrix factorization methods. In these empirical studies, the methods are evaluated by examining the distance between the prediction or recovery results with the original data, under certain metrics, typically the Frobenius norm. It will be interesting to carry out studies that are more directly connected to application cases in electricity generation and consumption. That is, there can be statistics that are more important for the electric network, and thus it will be more insightful to examine the methods using metrics that are directly linked to these important statistics.

For example, one potential application is to use the matrix factorization method to estimate the daily peak in power in local transformers from residential smart meter data. The setting is the following: smart meters can provide daily electricity consumptions to the distribution system operator (DSO) with individual offsets. The DSO can then run the algorithms proposed in this thesis to estimate the individual electricity consumption at a much finer scale (10-minute interval for example). By aggregating the estimated fine-scale consumption time series at each transformer, the DSO has an estimation to the moment of the daily peak and the level of the peak. Preliminary studies showed this approach to be promising. Perspectives in this application include fine tuning the evaluation metrics so as to reflect the real concern of the DSO. Moreover, to better adapt the methods of this thesis to this application, it can be interesting to change the objective function in the optimization problem to directly estimate the peak.

Methodological perspectives

Estimation with both spatial and temporal aggregates The empirical studies in this thesis used either spatial (in Chapters 3 and 4) or temporal aggregates (in Chapters 6 and 7).

However, the algorithm proposed in Chapter 7 can in fact use general linear measurements. CHAPTER 10. PERSPECTIVES It would then be interesting to combine both partial observations together in one empirical study to evaluate the potential gain. In terms of electricity applications, we can potentially improve the individual consumption estimation, by using not only temporal aggregates but also temporally fine-scaled measures on the network, such as measured from the distribution network.

Given that the spatially and temporally aggregated data come from different sources, these two parts of data might be incoherent. It would be important to relax the linear equality constraint to deal with this situation.

Usage of social network data

It is also possible to generalize the approach to social network data. The socio-demographic information for individuals considered in this thesis is in the form of individual features, such as client type, age, and occupation. It can be interesting to use a form of social networks, which encode the links between individuals. Social networks are important in recommendation systems, as people who are connected over a social network tend to have similar tastes. In electricity consumption, a similar mechanism might take place.

Correlation between the data and the measurement operator

In most matrix sensing studies, it is supposed that the measurement operator is independent of the matrix. In matrix completion, the standard setting is that the observation pattern is uniformly random over all matrix entries. When the matrix to be completed is also considered to be a realization of a random distribution, it is often supposed to be independent from the observation.

However, this hypothesis is rarely verified. In a recommendation application, for example the MovieLens100k application considered in Chapter 7, it is often the case that users would give ratings to a number of correlated items. Moreover, users tend to try items on which they have an prior positive opinion. Therefore, the observation pattern should be correlated with the matrix to be completed. In the electricity consumption estimation problem, the correlation could be of following forms:

• Temporally aggregated data: particularly with periodic observations, if the observation period is a multiple of one of the intrinsic periods of the data, there might be interference.

Variations of a frequency similar to that of the observations can be masked. Empirically, we observed in Chapter 6 that periodic observations have lower accuracy. This result could stem from the correlation between the observation patterns and the data.

• Spatially aggregated data: correlated individuals can be clustered spatially. If the observations are spatially aggregated consumption over similar individuals, and the estimation does not take this correlation into account, the total variance can be overestimated by the model.

Deep learning with partial data

In the end of Chapter 8, we discussed briefly that nonnegative matrix factorization is a very shallow neural network with ReLU activation function. This fact thus raises the question whether the methods developed in this thesis can be generalized to more complex neural network architectures, or inversely, if deep learning can be modified to use partial data such as those considered in this thesis. There has been some interest in using deep learning jointly with compressive sensing techniques in the image processing community recently [START_REF] Kim | Accurate Image Super-Resolution Using Very Deep Convolutional Networks[END_REF][START_REF] Dong | Image Super-Resolution Using Deep Convolutional Networks[END_REF][START_REF] Bora | Compressed Sensing Using Generative Models[END_REF]. One important challenge in generalizing from matrix sensing using deep learning lies in finding a good representation for the observations which is compatible to the neural network architecture and computation methods.

Theoretical perspectives

We have proposed algorithms that are convergent to stationary points. As mentioned in Chapter 5, new results are emerging in global optimum guarantees for non-convex problems. In general, NMF is an NP-hard problem, and therefore, there can not be general guarantees. Recent results have shown that alternating direction algorithms can have global guarantees when one of the factors is initialized close to the ground truth [91,98]

Details

The entity names are supposed to be • either in the form "xxxSyyy" where "xxx" the higher level entity name, "S" is a separator character (parameter sep) and "yyy" is the lower level entity name.

• or in the form "xxxxxyyy" where the first five characters are the higher level entity name.

Value

an integer vector of the same length as entity_names specifying how to aggregate the entities.

This can be used either in aggregated evaluation functions such as increment_SPaggregate_evaluate or as an intrinsic property of a dataset so defined while creating an incremental_data object.

create_formula_string compare_two_matrix Calculate the deviation between the data matrix and the estimation matrix

Description

Calculate the deviation between the data matrix and the estimation matrix Usage compare_two_matrix(data, est, metrics) implemented_matrix_comparison_metrics()

Arguments data, est

The two matrices to be compared metrics A character vector specifying the metrics to be used. Should be a subset of implemented_matrix_comparison_metrics.

Value

A named vector of the distance/deviation between the two matrices.

Examples

compare_two_matrix(example_matrix, jitter(example_matrix), implemented_matrix_comparison_metrics())

create_formula_string

Create character formula from a vector of variable names

Description

This function creates a string which can be later transformed to the right part of an R formula. It is primarily used to create formula_params objects to be used in regression_models.

Usage

create_formula_string(var_names, type = "gam", k = 15, numeric) empty_formula

Arguments var_names

The vector of variables names to be put in the formula type Either "gam" for smooths or anything else for a lm kind of formula.

k

The dimension of the spline basis when using "gam". numeric A vector of the same length as var_names specifying if the variables are numeric. If not supplied, all variables are supposed to be numeric. If this is not the case, then there will be problems with smooth terms in estimation. empty_formula

See

Create an empty formula to be used as a placeholder

Description

This object is to be used when no regression componnent is to be estimated. It is basically an empty formula_params.

Usage

Arguments

x an object of the class index_data, incremental_data, or linear_measurement_data. rho either NULL or a vector of length n. When it is provided, it is the vector of temporal autocorrelation threshold to be imposed on each of the columns. lambda either a number which is 0 or positive, or "optimal". This is the penalization parameter to be used to impose autocorrelation constraint. ... additional parameters to pass on to the estimation routine.

Details

When the data are temporal aggregates (index_data), constant interpolation is returned.

If rho and lambda is provided, the projection of 0 matrix with autocorrelation penalization is returned.

When the data are general linear measurements, the orthogonal projection of 0 matrix into the linear subspace of matrices satisfying the data constraints is returned. For example, for matrix completion, this is equivalent to replacing every missing value with 0.

Value

Description

The object of the class exp_data is to be used in experiments with many parameters. Therefore, specific information concerning the data is to be passed through this function. The rows 1:m_train + m_rho of the data is to be used as training data.

Usage

The rows 1:m_test + m_rho+m_train is to be used as test data (when regression models are estimated). The same for n_train and n_test for the columns. formula a list of formula and additional parameters to pass on to specify regression models. Several pairs of formulae can be specified simultaneously and will result in separate model estimations by experiment. Each element of formula is of the format list(X = formula_params(""), Y = formula_params("")).

The fields X and Y are for row and column features respectively.

Value

an exp_data object.

experiment

• tol: tolerance for deciding whether to stop.

• maxit: the maximal number of iterations.

• minit: the minimal number of iterations.

• verbose: whether to provide more reporting in the computation or not.

• fixed_X, fixed_Y: whether X or Y should be fixed in the iterations.

• maxtime: maximal computing time allocated in seconds.

• ITER_MAX, ITER_MIN: the maximal and minimal number of iterations for NeNMF.

Value

In a typical experiment run, this function saves all estimated models and the model estimation result tables on the disk.

It also returns a data.frame of standard model evaluation metrics. Apart from information on the models that are estimated, the output table contains outputs of the internal function evaluate_output_list applied on the train data, the column test data, the row test data, and the column-row test data. The latter three are present whenever they are applicable: with m_test, n_test strictly larger than 0, in presence of exogenous variables, and with algorithms which estimate regression_models.

See Also

evaluate_output_list to know more about the output evaluation metrics. rho either NULL or a vector of length n. When it is provided, it is the vector of temporal autocorrelation threshold to be imposed on each of the columns. lambda either a number which is 0 or positive, or "optimal". This is the penalization parameter to be used to impose autocorrelation constraint.

Examples

ground_truth the ground truth matrix of incremental data. Value a factorization_model object. It has all fields of matrix_model, and also the following additional ones:

NNDSVD

• X,Y: the two factors.

• hist_obj, hist_prjg, hist_error, hist_time: reporting values throughout the iterations.

• iter: number of iterations used.

• niter: number of inner iterations using NeNMF algorithm.

• convergence: whether the algorithm arrived at convergence or the maximal number of iterations is attained.

• X_iterations, Y_iterations: the factors throughout the iterations. They are NULL if save_iterations = FALSE.

• other fields which are only applicable to regression_models. Arguments output_table an output_table returned by experiement, or obtained by merging several these outputs.

Details

find_optimal_rank finds the optimal rank for each combination of all other experiment parameters.

duplicate_complete_sampling_results duplicates the rows which are estimated using complete matrix observations in the table returned by experiment and changes the sampling_type value to "periodic" and "random". This is to help us plot error rates depending on the sampling rate.

transform_formula_to_factors Transform formulas in the output table into factors so that it would be easy to produce graphics.

Value

a post-processed output table. X the first factor matrix.

See

X_prime the second factor matrix.

metrics Which metrics to use in the evaluation? This character vector should be a subset of implemented_matrix_comparison_metrics(). output1, output2

the two outputs to compare known extra set of observation index, if we want to evaluate the two on mask with a different one than that included in output1. perm whether to output the found permutation or not.

in_between the pairwise distance between the models. rank the rank of the models these_models the list of models.

increData the incremental_data against with to evaluate model predictions.

aggreg whether to calculate on an aggregation scale.

verbose wheter to print the results. generate_exo_mat prediction_type Which predict methods to use in the evaluation? This character vector should be a subset of implemented_prediction_types(). newdata When evaluating a regression_model, provide new data through this variable to test the prediction of the regression_model. As in predict.regression_model, if supplied, this should be a named list with fields X and/or Y.

Details

For a list of model estimated with the same parameters, evaluate_output_list from different random initializations, this function outputs the average performance metrics for matrix recovery and prediction. evaluation_matrix_model uses this function to actually calculate the error rates.

When evaluating on train data, evaluate_output_list also calculates several metrics which evaluates the distance between the models estimated from different random initialization. This is done by calling compare_two_outputs and how_many_stationary_points, which compare both the matrix recovery and the factors, and estimate the number of distinct stationary points in the list. Even more internally, they call finds_permutation and compare_nmf_factors to permute the factors and actually compare them.

generate_exo_mat

Generate some random features

Description

This is mostly used in development and for creating simulations.

Usage

Details

As parameters, provide a subset of comma-separated character strings of the following keywords:

• "octave": imported methods from octave ("grmf" and "trmf")

• "factorization": matrix factorization methods without regression

• "autocorrelation": algorithms that can be used jointly with autocorrelation.

• "regression": regression-enhenced algorithms

• "benchmark": benchmarks incremental_data Arguments increData, obs either the incremental data matrix increData or the observations on the data matrix obs. One of cumData, increData, and obs has to be supplied for evaluate_on_mask. period interger. he aggregation scale. lissage a number between 0 and 1. The smoothing parameter. cumData, increData the matrix of cumulative or incremental data. One of these must be supplied.

If both supplied, will use increData. cumEst, increEst the matrix of cumulative or incremental estimation. One of these must be supplied. If both supplied, will use increEst. metrics a character vector specifying the comparison metrics to be calculated. linkage_vector integer vector specifying how to aggregate columns of estimation and data matrix. exo_variable exogeneous variables associated with the data matrix. Supposed to be a named list with (optional) fields X and Y. aggreg a numeric vector specifying how to aggregate the data matrix. Must be of length ncol{x}. This is optional.

Details

An index is a single-integer representation of position of an entry in a matrix (example:

M[index]).
A subscript is a pair-of-integer representation (example: M[i,j] where (i,j) is the subscript).

ind2sub translates indexes to subscripts.

sub2ind translates subscripts to indexes. ## S3 method for class sampling_scheme index_data(x, incremental_data, ...) ## S3 method for class numeric index_data(x, sample_data, m, n, data_id, ...)

Examples

Arguments

x either a sampling_scheme object or a numeric vector of the index of the sampled data. If x is a sampling_scheme, the matrix to be sampled should be supplied by incremental_data. If x is a numeric vector, sampled_data should be supplied. incremental_data an objecto of the incremental_data class to be sampled.

sample_data the observations of indexes (cumulative values). When supplied, should match the numeric vector x. m, n the dimension of the incremental data matrix to be estimated. In other words, they are the dimensions of diff(index_data), not of the sampled accumulated index matrix.

Value

a list with two fields:

• matrix_A: a sparse matrix representing the measurement operator

• measurements: a vector of length dim{matrix_A} [1]. data either an incremental_data or a linear_measurement_data object. The data is to be sampled with x.

Details

• If x is a numeric vector, matrix_A and dim_data should be supplied.

• If x is a sampling_scheme, data should be supplied, and it shoud be either incremental_data or linear_measurement_data.

• If x is a measurement_matrix, data should be of the class incremental_data.

• If x is a Matrix::sparseMatrix, it is supposed that it is a matrix with missing entries. The linear measurement operator is one which samples all observed entries in x. More specifically, x should not be a Matrix::nsparseMatrix (a binary pattern matrix), but one which actually contains the value of the observations. Otherwise there will be a warning.

Value

a linear_measurement_data object which includes the following fields:

• sample_data: the observations in the form of a numeric vector.

• measurement_matrix: a measurement_matrix of dimension length(sample_data)by-prod(dim_data)

• dim_data: the dimension of the matrix to be estimated.

• exo_variable: the exogenous variables. Named list with fields called X and Y. Should verify nrow(exo_variable$X) == dim_data [1] and nrow(exo_variable$Y) == dim_data [2].

• data_id: a character string. The name of the dataset.

See Also measurement_matrix normalize_factors • fixed_X, fixed_Y: whether X or Y should be fixed in the iterations.

• maxtime: maximal computing time allocated in seconds.

• ITER_MAX, ITER_MIN: the maximal and minimal number of iterations for NeNMF. • X0,Y0: initial values for X, Y • FuninR: If FuninR = "R" use the R version of subroutines. If FuninR = "cpp" use the C++ version. Default is "cpp". • ccdppitmax,ccdppthreads,ccdpplambda: the maximal number of iterations, the number of threads, and the penalization parameter for CCDPP algorithm.

rho either NULL or a vector of length n. When it is provided, it is the vector of temporal autocorrelation threshold to be imposed on each of the columns. lambda either a number which is 0 or positive, or "optimal". This is the penalization parameter to be used to impose autocorrelation constraint.

ground_truth the ground truth matrix of incremental data.

NNDSVD a logical to specify whether use SVD to initialize. save_iterations whether to save the factorizations during iterations. general_factorization

If TRUE does a general factorization without nonnegativity. If FALSE, which is the default value, does an NMF. IF is the string "random", randomly allow half of the iterations to explore negative values.

exo_variable A list of two fields named X and Y. Both are to be data.frames, of m or n rows, corresponding to the row individuals or the column individuals of the matrix to be recovered.

formula

A list of two fields named X and Y. Both are to be objects of the class formula_params.

... additional parameters to be passed to mgcv::gam for initializing the models when using HALS_models and gam. the sampling rate (0.2 for sampling 20% of data). For every sampling type except evolution_sampling, this should be a single number between 0 and 1. For evoluation_sampling, this should be a pair of sampling rates. sr [1] is the (periodic) sampling rate applied on regular users. sr [2] is the percentage of regular users. The rest of the users are observed completely (supposed to be users with Linky meters). seed the seed for sampling the phase.

Details

periodic_sampling proposes sampling entries which are at regular intervals (with random offset for each column).

random_sampling proposes sampling entries which are chosen uniformly at random amongst all entries in the matrix.

evolution_sampling proposes sampling entries at regular intervals for a part of the users, and complete observations for the rest of the users.

block_sampling supposes that a block of data is missing, and returns the entries that are not missing.

Value

The sampled (observed) entries' indexes. For mat, a matrix of dimension m-by-n, call mat[xxxxx_sampling(m, n, sr)] to retrieve the observed entries.

See Also

See factorization_model for examples. See plot_factors and plot_weights for additional parameters.

plot.regression_model Plot a regression_model

Description

This function plots the regression models (the models between factors and exogenous variables).

Usage plot.regression_model(x)

Arguments

x an object of class regression_model. x$algorithm should be either "HALSX1" or "factor_gam".

See Also

See regression_model for examples. The factor by which to order the individuals. For the graphic to be more informative, we reorder the individuals according to one of the factors.

See Also plot_factors

Examples plot_weights(example_matrix) regression_model rho either NULL or a vector of length n. When it is provided, it is the vector of temporal autocorrelation threshold to be imposed on each of the columns. lambda either a number which is 0 or positive, or "optimal". This is the penalization parameter to be used to impose autocorrelation constraint.

ground_truth the ground truth matrix of incremental data.

NNDSVD a logical to specify whether use SVD to initialize. save_iterations whether to save the factorizations during iterations. general_factorization

If TRUE does a general factorization without nonnegativity. If FALSE, which is the default value, does an NMF. IF is the string "random", randomly allow half of the iterations to explore negative values.

... additional parameters to passe on to NMF_increment.

Value

A regression_model object. Always has a field called models which consequently can be used in prediction. When the object returned is also a matrix_model or a factorization_model, corresponding fields are included.

See Also

formula_params to see details on how to provide parameters for the regression algorithms.

Description

Specify a sampling scheme with its type and parameters Usage sampling_scheme(sampling_rate, type = c("periodic", "random", "evolution", "block", "completion", "all"), seed = NULL) is.sampling_scheme(object)

3. 1

 1 Log-density of population in Lyon. . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Cloud of empirical semivariance log-density of population . . . . . . . . . . . . . 3.3 Empirical variogram log-density of population . . . . . . . . . . . . . . . . . . . . 3.4 Variogram models included in 'gstat' . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Fitted variogram models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3.6 Grid on the domain of interest. The log-density of population is assigned to the square containing the geometric center of each IRIS. . . . . . . . . . . . . . . . . 3.7 Kriging prediction value of log-density of population. . . . . . . . . . . . . . . . . 3.8 Kriging variance of log-density of population. . . . . . . . . . . . . . . . . . . . . 3.9 Block-kriging value of log-density of population. . . . . . . . . . . . . . . . . . . . 3.10 Block-kriging variance of log-density of population. . . . . . . . . . . . . . . . . . 3.11 Spatial average of the weekly consumption and an autoregressive model fitted on it. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Autocorrelation and partial autocorrelation function. . . . . . . . . . . . . . . . . 3.13 Spatial variograms fitted on the four weeks separately. . . . . . . . . . . . . . . . 3.14 A wireframe plot of the empirical and fitted spatio-temporal variogram. . . . . . 3.15 The empirical and fitted temporal variogram by temporal lag. . . . . . . . . . . . 3.16 Weekly electricity consumption of the MV feeders of the first four weeks. . . . . 3.17 Kriged weekly IRIS electricity consumption of the first four weeks . . . . . . . . LIST OF FIGURES 4.1 Overlapping of IRIS and area supplied by MV feeders in the 7th municipal district of Lyon. IRIS are delimited by black borders. Points of different colors represent transformers connected to different MV feeders -each color corresponds to one MV feeder. In this example, source zones are areas served by transformers of the same color, and target zones, the IRIS, are the areas delimited by the black lines.4.2 Electric load of an MV feeder during three years (left) and three weeks (right). . 4.3 Cluster model: GAM estimation of annual cycle for the consumption at 10 a.m. on weekdays. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Cluster model: GAM estimation of annual cycle for the consumption at 10 a.m. on weekend. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Cluster model: GAM estimation of the SR effect throughout the day. . . . . . . . 5.1 An illustration of nonnegative matrix factorizaiton for a matrix of daily electricity consumption series. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 An illustration of masks which correpond to the meter readings on individual meters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Left rank-14 factor estimated by NeNMF . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Right rank-14 factor estimated by NeNMF . . . . . . . . . . . . . . . . . . . . . . . 5.5 Left rank-14 factor estimated by HALS . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Right rank-14 factor estimated by HALS . . . . . . . . . . . . . . . . . . . . . . . 6.1 Recovery examples from each of the four datasets. The data matrix is sampled with the scheme random, with 20% sampling rate in this experiment. The black line is the ground truth. The green line the recovery made by unpenalized NeNMF. The red line is the recovery made by penalized NeNMF. . . . . . . . . . . . . . . .

7. 1

 1 Dividing the data matrix into four submatrices: one for generating observations to estimate the model, the other three for prediction. . . . . . . . . . . . . . . . . 7.2 Execution time per iteration of Algorithm 4 and Algorithm 5 without exogenous variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Total execution time of Algorithm 4 and Algorithm 5 without exogenous variables x d'électricité, les données de consommation peuvent être utilisées par les opérateurs de réseau de transport (RTE en France) et de distribution pour mesurer les consommations/productions des portefeuilles des différents responsables d'équilibre tels les producteurs d'électricité, les gros consommateurs, les agrégateurs d'effacement ou les traders.

. 1 .

 1 Chaque colonne représente un individu et chaque ligne une période. Les deux sources de données partielles sont des sommes d'entrées de matrice : le rectangle orange représente une agrégation temporelle, et le rectangle vert représente une agrégation d'individus (spatiale par exemple).

Figure 1 . 1

 11 Figure 1.1: Représentation matricielle de la variable d'intérêt. Dans cette thèse, une entrée de la matrice est souvent la consommation électrique d'un individu (une colonne), durant une période (une ligne). Les mesures de capteurs sur le réseau sont souvent des agrégations d'individus (somme d'entrées dans le rectangle vert). Les relevés de compteurs individuels sont souvent des agrégations temporelles (somme d'entrées dans le rectangle orange).

Figure 2 .

 2 Figure 2.1: The matrix representation of the variable of interest. Typically, an entry of the matrix is the electricity consumption of an individual (a column), during a period (a row). Measures taken by sensors on the network are often aggregations of individuals (sum of the entries in the green box). Measures read on user meters are often temporally aggregated (sum of the entries in the orange box).

  A more model-based approach to solve areal interpolation is by considering the observations as the average value of the variable of interest, Z(D 1 ), ..., Z(D n ), on a collection of subsets, D 1 , ..., D n , of the domain D. ([41, Chap. 6], [54, Chap. 9], [48, Chap. 4]).

Figure 3 .

 3 Figure 3.4 shows the values of exponential, Gaussian, and a 1.5-order Matérn variogram with the same range and sill. The graphics are generated by the gstat Package in R [62].

Figure 3 . 1 :Figure 3 . 2 :Figure 3 . 3 :Figure 3 . 4 :

 31323334 Figure 3.1: Log-density of population in Lyon.

Figure 3 .

 3 Figure 3.7 shows the kriging predictions for each little squares on the grid. The three estimated variograms of Figure 3.5 are used. In comparison, a linear interpolation based on triangulation of the data points is also shown. Note that in linear interpolation, only points within the convex hull of the data points are evaluated.

Figure 3 . 5 :Figure 3 . 6 : 4 Figure 3 . 7 :Figure 3 . 8 :

 353643738 Figure 3.5: Fitted variogram models

Figure 3 . 9 :

 39 Figure 3.9: Block-kriging value of log-density of population.

Figure 3 . 10 :

 310 Figure 3.10: Block-kriging variance of log-density of population.

21 -Figure 3 . 11 :

 21311 Figure 3.11: Spatial average of the weekly consumption and an autoregressive model fitted on it.

Figure 3 . 12 :

 312 Figure 3.12: Autocorrelation and partial autocorrelation function.

Figure 3 . 13 :Figure 3 . 14 :

 313314 Figure 3.13: Spatial variograms fitted on the four weeks separately.

Figure 3 . 16 :

 316 Figure 3.16: Weekly electricity consumption of the MV feeders of the first four weeks.

4. 3 .Figure 4 . 2 :

 342 Figure 4.2: Electric load of an MV feeder during three years (left) and three weeks (right).

4. 3 .

 3 APPLICATION ON REAL AND SYNTHETIC LOAD DATA q q

Figure 4 . 3 :

 43 Figure 4.3: Cluster model: GAM estimation of annual cycle for the consumption at 10 a.m. on weekdays.

Figure 4 . 4 :

 44 Figure 4.4: Cluster model: GAM estimation of annual cycle for the consumption at 10 a.m. on weekend.

Figure 4 . 5 :

 45 Figure 4.5: Cluster model: GAM estimation of the SR effect throughout the day.

Figure 5 . 1 :

 51 Figure 5.1: An illustration of nonnegative matrix factorizaiton for a matrix of daily electricity consumption series.

. 5 )

 5 The design matrices A 1 , ..., A N are called masks. The matrices on the right part of Figure 5.2 show a subsection of the masks correponding to the indices on individual electricity meters shown on the left of Figure 5.2.

Figure 5 . 2 :

 52 Figure 5.2: An illustration of masks which correpond to the meter readings on individual meters.

Figure 5 . 3 :Figure 5 . 4 :

 5354 Figure 5.3: Left rank-14 factor estimated by NeNMF

Figure 5 . 5 :Figure 5 . 6 :

 5556 Figure 5.5: Left rank-14 factor estimated by HALS

  is the index set over h(d) consecutive periods of Consumer n d , starting from Period t 0 (d) + 1. Each measurement covers a disjoint index set. All entries of X are not necessarily involved in the measurements.

Figure 6 . 1 :

 61 Figure 6.1: Recovery examples from each of the four datasets. The data matrix is sampled with the scheme random, with 20% sampling rate in this experiment. The black line is the ground truth. The green line the recovery made by unpenalized NeNMF. The red line is the recovery made by penalized NeNMF.

Figure 7 . 3 :

 73 Figure 7.3: Total execution time of Algorithm 4 and Algorithm 5 without exogenous variables

Figure 7 . 8 :

 78 Figure 7.8: Prediction performance on real Portuguese electricity data

Contents 9 . 1

 91 Data representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 9.1.1 incremental_data: A fully-observed data matrix . . . . . . . . . . . . 114 9.1.2 index_data: Temporal aggregates . . . . . . . . . . . . . . . . . . . . 114 9.1.3 linear_measurement_data: General linear measurements . . . . . . . 115 9.1.4 exp_data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 9.2 matrix_model: Estimate a single model . . . . . . . . . . . . . . . . . . . . . . 116 9.3 Model evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 9.4 Run an experiment to estimate several models simultaneously . . . . . . . . . 119 9.5 Distributed computation in model estimation . . . . . . . . . . . . . . . . . . . 120

#

  create a rank-2 nonnegative 10-by-15 matrix as example m = 10 n = 15 r = 2 example_matrix = matrix(rexp(m * r), ncol = r) %*% matrix(rexp(n * r), ncol = n)

#

  A measurement matrix that samples the first five entries of a matrix. matrix_A = measurement_matrix( "matrix", matrix_A = Matrix::sparseMatrix(i = 1:5, j = 1:5, dims = c(5, m*n))) lin_m_d = linear_measurement_data( x = matrix_A,

9. 2 .Figure 9 . 1 :Figure 9 . 2 :

 29192 Figure 9.1: The original matrix and the matrix estimated by NeNMF using 50% random temporal aggregates.

Figure 9 . 3 :

 93 Figure 9.3: Comparison of the models estimated by calling experiment

9. 5 .

 5 DISTRIBUTED COMPUTATION IN MODEL ESTIMATION constraint are distributed through the multiple cores.

  Also formula_params and regresion_model Examples create_formula_string(LETTERS[1:5]) create_formula_string(LETTERS[1:5], numeric = c(FALSE, FALSE, TRUE, TRUE, TRUE)) create_formula_string(LETTERS[1:5], type = "lm")

  ex_sp = exp_spec(r = 2:3, algorithm = c("NeNMF", "HALS", "HALS_model"), lambda = 0, random_start = 2, general_factorization = FALSE) ## Not run: example_exp = experiment(example_exp_data, example_sampling_scheme, ex_sp) example_exp = find_optimal_rank(example_exp) ## End(Not run) factorization_model Fit a factorization model Description Fit a factorization model with index data or incremental data factorization_model • FuninR: If FuninR = "R" use the R version of subroutines. If FuninR = "cpp" use the C++ version. Default is "cpp". • ccdppitmax,ccdppthreads,ccdpplambda: the maximal number of iterations, the number of threads, and the penalization parameter for CCDPP algorithm.

Examples

  fac_mod = factorization_model(example_index_data, r = 2, algorithm = "HALS") # predict with completion = TRUE or FALSE predict(fac_mod, completion = TRUE) predict(fac_mod, completion = FALSE) # plot the model plot(fac_mod) find_optimal_rank find_optimal_rank Post-processing of experiment outputs Description These functions help post process experiment outputs.

  Also experiment for examples. finds_permutation finds_permutation Evaluate a list of matrix_models Description Internal functions used in experiment to report estimated model performance. Usage finds_permutation(correlation_matrix) compare_nmf_factors(X, X_prime, metrics = implemented_matrix_comparison_metrics()[1]) compare_two_outputs(output1, output2, known = NULL, perm = FALSE, matrix_A = NULL) how_many_stationary_points(in_between, rank) evaluate_output_list(these_models, increData, aggreg, verbose = FALSE, prediction_type = implemented_prediction_types(), metrics = implemented_matrix_comparison_metrics(), newdata = NULL) Arguments correlation_matrix the correlation matrix between two factors matrices of the same size used to find the permutation of factors which minimizes their distance.

ValueA

  named list with • exo_matrix: the variables • features: the transformed variables • formula_alone: the formula used to transform variables • formula_caret: formula_alone without the splines get_data_test get_data_test Subsample or aggregate the data matrix to desired frequency and dimension Description This function subsamples or aggregates a data matrix whose columns are time series to a desired frequency. It can also calculates the cumulative sums of the columns and/or normalizes the columns to have mean 1. Usage get_data_test(data_matrix, base_freq = 1, row_number = NULL, normalize = TRUE, cumul = TRUE, fun = mean) Arguments data_matrix a numeric matrix. base_freq the frequency desired.row_number the intended row_number.

  normalizewhether the data is to be normalized cumul whether the cumulative index is to be returned fun the aggregation function to be used. Default is mean.DetailsNAs are removed in calculating the aggregation values. Value a matrix of dimension row_number-by-ncol(data_matrix). (rexp(m * r), ncol = r) %*% matrix(rexp(n * r), ncol = n) system.time(get_data_test(example_matrix, 100)) system.time(aggregateData(example_matrix, 100, mean)) ## End(Not run) get_formula_param get_formula_param Retrieves parameters from a formula_params Description Retrieves parameters from a formula_params Usage get_formula_param(x, what) Arguments x the formula_params from which to get parameters. what the name of the parameters to be retrieved. Value a vector of parameters which match what. NULL is returned if nothing is found. implemented_algorithms The implemented algorithms Description Lists the algorithms using keywords. Usage implemented_algorithms(...)

  known the index of the observations on the matrix. Used in evaluate_on_mask. cumEst the matrix of cumulative estimated consumption incremental_data Put the data matrix in form Description Put the data matrix in form Usage incremental_data(x, exo_variable = NULL, aggreg = NULL, data_id) ## S3 method for class incremental_data dim(x) Arguments x the numeraic data matrix.

  number of rows in the matrix ind vector of integers. The indexes to be translated.r, cthe integer vectors of row and column subscripts.

  subs = ind2sub(10, 1:20) inds = sub2ind(10, subs[,1], subs[,2]) index_data index_data Create an object of index data Description Create an object of index data, either with incremental data and a sampling scheme, or directly with index data. Usage index_data(x, ...) is.index_data(x) ## S3 method for class index_data dim(x)

  (rexp(m * r), ncol = r) %*% matrix(rexp(n * r), ncol = n) # create an incremental_data object incre_d = incremental_data(x = example_matrix, data_id = "exponential") seed = 42 # we will sample incre_d with a random sampling_scheme # with a 50% sampling rate and a fixed seed sample_incre_d_at = sampling_scheme( sampling_rate = 0.5, type = "random", seed = seed)# create an index_data object index_d = index_data(x = sample_incre_d_at, incremental_data = incre_d) lapply(index_d,head) index_d2 = index_data(x = index_d$sample_index, sample_data = index_d$sample_data, m = index_d$dim_data[1], n = index_d$dim_data[2], data_id = index_d$data_id) index_data_to_matrix_measurement Convert index_data to a measurement matrix and the observations Description Convert index_data to a measurement matrix and the observations Usage index_data_to_matrix_measurement(measurements_index_data) sparse_matrix_to_matrix_measurement(sparse_matrix) Arguments measurements_index_data an index_data object. sparse_matrix a sparse matrix where the zero entries are unobserved.

normalize_factorsA

  Normalize the factors in a factorization_model Description This function rescales output$X to have columns which have averages of one. output$Y is modified accordingly to keep the model equivalent. factorization_model with normalized factors. periodic_sampling Helper functions to subsample the matrices Description These functions are returned wrapped inside sampling_schemes. Usage periodic_sampling(m, n, sr, seed = NULL) random_sampling(m, n, sr, seed = NULL) evolution_sampling(m, n, sr, seed = NULL) block_sampling(m, n, sr, seed = NULL) Arguments m, n integers. The number of rows and columns in the matrix. sr

Description

  Plots the factors and the weights of a factorization_model, using plot_factors and plot_weights.Usageplot.factorization_model(x, ...)Argumentsx an object of class factorization_model, aka, has the fields X and Y. ... parameters passed to plot_factors and plot_weights.

#

  in a factorization model. This function takes a numeric matrix and plot each column as a curve. Underlying this function is the idea that factors that are obtained are time series. Usage plot_factors(W, start = NULL, interval = NULL, ...) Arguments W a numeric matrix (factors obtained in NMF typically). start, interval options to put time x-axis values for factors. start is a POSIXct. interval is a difftime. We can also specify the periods covered by the factors plot_factors(example_matrix, start = as.POSIXct(Sys.Date()), interval = as.difftime(1, units = "days")) plot_weights Plot weights Description Plot a numeric matrix by interpreting each column as weights of factors for one individual. Each column of H is renormalized to sum 1. Usage plot_weights(H, order_by_factor = 1, ...) Arguments H a numeric matrix (weights obtained in NMF). order_by_factor

#

  When not specified explicitly, # predict treats a regression_model that is also a factorization_model as # a regression_model. # Predict a regression_model without newdata to see the prediction # values of the regression models. predict(reg_mod) # Call predict.factorization_model explicitly on regression_model to see # prediction by the factorization_model predict.factorization_model(reg_mod, completion = TRUE) predict.factorization_model(reg_mod, completion = FALSE)# Predict using newdata predict(reg_mod, newdata = list(X = example_exogenous_variables$X[1:2,])) predict(reg_mod, newdata = list(Y = example_exogenous_variables$Y[1,,drop=FALSE])) predict(reg_mod, newdata = list(X = example_exogenous_variables$X[1:2,], Y = example_exogenous_variables$Y[1,,drop=FALSE]))# When not specified explicitly, it is the regression_model part that is # plotted.

  relative_diff_norm plot(reg_mod) # We can force it to plot factorization_model. plot.factorization_model(reg_mod) relative_diff_norm Relative difference in norm between two matrices Description Relative difference in norm between two matrices Usage relative_diff_norm(matrix1, matrix2, type = "F") matrices to compare. type currently unused. For future development can be used to inplement other metrix than Frobenius norm. sampling_scheme A wrapper class for the sampling schemes

  1.1 Représentation matricielle de la variable d'intérêt. Dans cette thèse, une entrée de la matrice est souvent la consommation électrique d'un individu (une colonne), durant une période (une ligne). Les mesures de capteurs sur le réseau sont souvent des agrégations d'individus (somme d'entrées dans le rectangle vert). Les relevés de compteurs individuels sont souvent des agrégations temporelles (somme d'entrées dans le rectangle orange). . . . . . . . . . . . . . . . . . . . . . . . . . .

  Les jeux de données réelles de consommation électrique sont les suivants : 1.4. CONTRIBUTIONS À LA FACTORISATION DE MATRICE NONNÉGATIVE • Consommation électrique française (Données fournies par Enedis) la consommation journalière des 636 départs HTA (haute tension A) dans la région autour de Lyon. Chaque individu dans ce jeu de données correspond à la consommation d'un groupe d'environ 1500 clients. Les données brutes sont disponibles au pas de 10 minutes, pendant 3 ans. Pour les expériences des chapitres 6 et 7, la consommation est agrégée au pas journalier.

Prédiction pour les nouveaux clients et/ou les nouvelles périodes

Comme schématisé dans la figure

1

.1, deux applications sont considérées : Reconstitution de consommation passée/estimation de score de films (chapitre 6 et 7) Nous mesurons les agrégations temporelles sur un jeu de données de consommation en tirant au hasard un certain nombre de dates d'observation. Nous utilisons une méthode pour estimer la matrice de consommation, et mesurons la performance en comparant les estimations avec la consommation réelle. Pour le jeu de données de score de films, nous tirons au hasard un sous-ensemble de score comme observations, et complétons la matrice des scores. Ensuite nous comparons les estimations avec les scores qui n'ont pas été observés. (chapitre 7) Nous cachons toutes les observations sur certains individus et certaines dates (ou certains films dans le cas de MovieLens), et appliquons la NMF avec variables exogènes sur les observations restantes. Nous produisons ensuite des prédictions pour les dates et les individus cachés. Nous les comparons avec la réalité pour mesurer la performance de la méthode.

  Contents 2.1 Motivation for electricity consumption estimation and prediction . . . . . . . . 11 2.2 Problem formalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3 Classical and novel solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.3.1 Spatio-temporal kriging . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.3.2 Prediction of multiple electricity consumption series . . . . . . . . . . 14 2.3.3 Matrix-factorization based methods . . . . . . . . . . . . . . . . . . . . 15

  Variogram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Kriging in univariate change of support problem . . . . . . . . . . . . . . . . 3.2.1 Kriging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Block kriging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3 Areal interpolation beyond kriging . . . . . . . . . . . . . . . . . . . . 3.3 Multivariate kriging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Cross-covariance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Cokriging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Change of support in spatio-temporal model . . . . . . . . . . . . . . . . . . . Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	Chapter
	3	Spatio-temporal change of sup-
		port
	Contents
		3.1 Introduction to spatial statistics . . . . . . . . . . . . . . . . . . . . . . . . . .
		3.1.1 Basics in geostatistics . . . . . . . . . . . . . . . . . . . . . . . . . . .
		3.1.2 3.4.1 Spatio-temporal covariance and kriging . . . . . . . . . . . . . . . . . .
		3.5 Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
		3.5.1 Variogram of the density of population in Lyon . . . . . . . . . . . . .
		3.5.2 Kriging of a univariate variable . . . . . . . . . . . . . . . . . . . . . .
		3.5.3 Block kriging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
		3.5.4 Kriging of spatio-temporal data . . . . . . . . . . . . . . . . . . . . . .
		3.6

Table 3 . 1 :

 31 Estimated parameters of variogram models

		Nugget (b) Sill (c)	Scale (a)
	Exponential	0.015	4.469 236729.503
	Gaussian	0.019	0.075	3238.211
	Matern	0.020	0.215	3731.366

Table 3 . 2 :

 32 Estimated parameters of Matern variogram models for the first four weeks of consumption data

	Week Nugget (b) Sill (c) Scale (a)
	1	0.244	0.081	633.972
	2	0.252	0.143	631.750
	3	0.283	0.111	455.967
	4	0.258	0.137	440.876

Table 3 .

 3 The empirical and fitted temporal variogram by temporal lag.

	##	min	max
	## x 402208.4 412869.6
	## y 4879405.6 4889993.4

3: Estimated parameters of spatio-temporal variogram models

Model component Nugget (b) Sill (c) Scale (a)

Table 3 . 4 :

 34 Error rates of the kriging estimation compared to the simulation consumption of IRIS (ground truth)

		Separable sumProduct spatialKrige
	MAPE	1.614	1.643	1.513
	RRMSE	4.479	4.549	4.229

Figure 3.17: Kriged weekly IRIS electricity consumption of the first four weeks

Table 4 . 1 :

 41 Rules of association of client type and distribution transformers. The + (or -) sign for either variable signifies that the value of the variable is larger (or smaller) than the median of all transformers.

		Apartments/CP -Apartments/CP +
	Warehouses/CP -	rural homes	apartment block
	Warehouses/CP +	industrial client	professional

MV feeder distribution zones with insufficient historical data and IRIS (Ilots Regroupés pour l'Information Statistique), an administration zone of approximately 2,000 habitants, used by the French National Institute of Statistics for releasing official socio-economic data. Note that feeder distribution zones and IRIS overlap.

Table 4 . 2 :

 42 Mean prediction error rate on the real dataset.

	Specification Model type	IRIS	CP
			clus-	clus-
			ters	ters
	Parametric Observed level	18.62	17.90
		(source)
	Parametric Observed level	19.13	18.54
		(target)
	Parametric Estimated	38.20	37.76
		level (target)
	GAM	Observed level	18.90	17.59
		(source)
	GAM	Observed level	19.38	18.28
		(target)
	GAM	Estimated	38.41	37.53
		level (target)

Table 4 . 3 :

 43 Mean prediction error rate on simulated consumption.

		IRIS clusters CP clusters
	MV feeders (source)	13.20	11.22
	MV feeders (target)	13.87	11.45
	IRIS estimation (target)	16.93	14.09

  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.2 Multivariate time series estimation and prediction as nonnegative matrix recovery . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 5.2.1 Nonnegative matrix factorization . . . . . . . . . . . . . . . . . . . . . 56 5.2.2 General model definition . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5.2.3 Applications of matrix recovery . . . . . . . . . . . . . . . . . . . . . . 58 5.3 Convergence of nonnegative matrix factorization algorithms . . . . . . . . . . 59 5.3.1 Global convergence of alternating minimization algorithms in non-

	to nonnegative
	matrix factorization
	Contents
	5.1

convex problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 5.3.2 What would global convergence mean? . . . . . . . . . . . . . . . . . . 60 5.3.3 Link with general matrix factorization . . . . . . . . . . . . . . . . . . 60 5.3.4 Empirical convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 5.3.5 Local optimality conditions . . . . . . . . . . . . . . . . . . . . . . . . 63 5.3.6 Projections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.3.7 Applying second order condition . . . . . . . . . . . . . . . . . . . . . 65

Table 5 . 1 :

 51 Relative error (RRMSE) in Frobenius norm of the factors produced by the NeNMF algorithm.

	Rank Objective function Error on F_r Error on F_c
	4	5e-04	0.0111	0.0145
	5	6e-04	0.0304	0.0420
	6	6e-04	0.0243	0.0369
	7	6e-04	0.0345	0.0406
	8	5e-04	0.0215	0.0325
	9	5e-04	0.0375	0.0509
	10	5e-04	0.0308	0.0410
	11	4e-04	0.0316	0.0444
	12	4e-04	0.0216	0.0270
	13	5e-04	0.0332	0.0508
	14	4e-04	0.0487	0.0577

Table 5 . 2 :

 52 Relative error (RRMSE) in Frobenius norm of the factors produced by the HALS algorithm.

	Rank Objective function Error on F_r Error on F_c
	4	1e-04	0.0038	0.0045
	5	8e-04	0.0245	0.0317
	6	8e-04	0.0169	0.0307
	7	5e-04	0.0176	0.0235
	8	7e-04	0.0231	0.0370
	9	7e-04	0.0228	0.0320
	10	5e-04	0.0196	0.0265
	11	6e-04	0.0197	0.0291
	12	6e-04	0.0304	0.0369
	13	8e-04	0.0206	0.0295
	14	1e-03	0.0350	0.0524

Table 9 . 1 :

 91 Included algorithms in meterModels

	Algorithm	Type
	HALS	Proposed method
	NeNMF	Proposed method
	HALSX1	Proposed method
	HALS_model	Proposed method
	softImpute	Reference method
	HALSX2	Reference method
	factor_gam	Reference method
	grmf	Reference method
	trmf	Reference method
	individual_gam	Reference method
	reduced_rank_regression Reference method
	empty_model

  0.258141By using evaluate_matrix_model, we can evaluate the model using several metrics. The last two fields in the returned list below, regression and prediction, are reserved for regression models (Chapter 7), hence are NA here.

		9.4. RUN AN EXPERIMENT TO ESTIMATE SEVERAL MODELS SIMULTANEOUSLY
	## $regression
	## Frob	ND NRMSE RMSE MAPE
	##	NA	NA	NA	NA	NA
	##
	## $prediction
	## Frob	ND NRMSE RMSE MAPE
	##	NA	NA	NA	NA	NA
		evaluate_matrix_model(nenmf_model, incre_d)
		## $completion
		##	Frob		ND	NRMSE	RMSE	MAPE
		## 0.25814097 0.07393886 0.39200128 0.77130529 0.04642855
		##
		## $factorization
		##	Frob		ND	NRMSE	RMSE	MAPE
		## 0.25814137 0.07410952 0.39200188 0.77130647 0.04655324
		##

Table 9 . 2 :

 92 Part of the evaluation table generated by experiment.

	r algorithm	completion.Frob factorization.Frob time_spent
	2 NeNMF	0.2499	0.2499	0.6408
	3 NeNMF	0.3899	0.3899	0.8724
	4 NeNMF	0.4723	0.4723	0.5726
	5 NeNMF	0.5118	0.5118	0.5602
	2 HALS	0.1961	0.1961	0.6340
	3 HALS	0.3137	0.3137	0.8288
	4 HALS	0.4015	0.4015	0.4019
	5 HALS	0.4470	0.4470	0.3711
	NA empty_model	0.6415	NA	0.0089
		Recovery error rate
		Execution time

  an empty_model which is a matrix_model with minimal fields.

		example_matrix exp_data
	reg_mod = regression_model(x = example_index_data,
	exp_data	exo_variable = example_exogenous_variables, Create a set of data to be used in an experiment
		formula = example_formula,
		r = 2,
		algorithm = "HALS_model")
	# For a regression_model, is no newdata is supplied,
	# the prediction error is calculated on original data, hence equal to
	# regression error.
	evaluate_matrix_model(reg_mod, example_incremental_data)
	evaluate_matrix_model(reg_mod, example_incremental_data,
		newdata = example_exogenous_variables)
	example_matrix	Example data structures
	Description
	Example data structures
	Usage
	example_matrix
	example_sampling_scheme
	example_index_data
	example_incremental_data
	example_linear_measurement_data
	example_exogenous_variables
	Examples
	model = empty_model(x = example_index_data) example_formula
	model$M model = empty_model(x = example_linear_measurement_data) example_exp_data
	model$M
	pred_empty_model = predict(model, completion = TRUE) Format
	# empty_model cannot predict with completion=FALSE, since it s not An object of class matrix with 13 rows and 11 columns.
	# a factorization_model
	is.null(predict(model, completion = FALSE))

  • X0,Y0: initial values for X, Y • FuninR: If FuninR = "R" use the R version of subroutines. If FuninR = "cpp" use the C++ version. Default is "cpp". • ccdppitmax,ccdppthreads,ccdpplambda: the maximal number of iterations, the number of threads, and the penalization parameter for CCDPP algorithm.

	Arguments
	exp_data	An exp_data object.
	sampling_scheme
		A sampling_scheme object.
	exp_spec	An exp_spec object.
	verbose	Whether messages reporting the start and the end of each model estimation
		should be printed in stdout.
	save_model	a character string specifying where to save the estimated models.
	eval_only	Should we estimate (potentially) new models or only evaluate old models
		saved in save_model.
	See Also
	experiment to run experiments.
	Examples
	ex_sp = exp_spec(r = 2:10)
	experiment	The main experiment function
	Description
	This function is the main function to use to run a series of experiments.
	Usage
	experiment(exp_data, sampling_scheme, exp_spec, verbose = TRUE)
	evaluate_past_experiement(exp_data, sampling_scheme, exp_spec,
	verbose = FALSE)

Value

An object of class exp_spec which wraps the provided parameters AND an exp_grid in which each row is a model to be estimated. Incompatible parameter combinations are removed by an internal function simplify_exp_grid.

  a logical to specify whether use SVD to initialize. save_iterations whether to save the factorizations during iterations. general_factorization If TRUE does a general factorization without nonnegativity. If FALSE, which is the default value, does an NMF. IF is the string "random", randomly allow half of the iterations to explore negative values.

	...	additional parameters to pass on to the estimation routine.

formula

A list of two fields named X and Y. Both are to be objects of the class formula_params.

  Arguments m number of rows in the exogenous variable table. d number of variables in the exogenous variable table. k dimension of the spline basis to be used to generate the features.

	generate_exo_mat(m, d, k = 3)

  , a sampling_scheme, a measurement_matrix, or a Matrix::sparseMatrix. matrix_A the measurement matrix used to generate the measures. Should be an object of class Matrix::sparseMatrix of dimension length(x)-by-prod(dim_data). dim_data a length-2 vector with the dimension of the measured matrix exo_variable optional exogenous variables. data_id the ID of the dataset.

			linear_measurement_data
	Examples linear_measurement_data(x, exo_variable = NULL,
		measurements_m = index_data_to_matrix_measurement(example_index_data) data_id = "")
	is.linear_measurement_data(x)
	## S3 method for class linear_measurement_data
	is.incremental_data dim(x)
			Reports whether x is an incremental_data object
	Arguments
	x	Description	a numeric vector
		Reports whether x is an incremental_data object
		Usage
		is.incremental_data(x)
		Arguments
		x	An object to test
		linear_measurement_data
			Create a data object with general linear measurements
		Description
		Create a data object with general linear measurements
		Usage
		linear_measurement_data(x, ...)
		## S3 method for class numeric
		linear_measurement_data(x, matrix_A, dim_data,
		exo_variable = NULL, data_id)
		## S3 method for class sampling_scheme
		linear_measurement_data(x, data,
		exo_variable = NULL, ...)
		## S3 method for class measurement_matrix
		linear_measurement_data(x, data, ...)
		## S3 method for class sparseMatrix

  FuninR: If FuninR = "R" use the R version of subroutines. If FuninR = "cpp" use the C++ version. Default is "cpp".• ccdppitmax,ccdppthreads,ccdpplambda: the maximal number of iterations, the number of threads, and the penalization parameter for CCDPP algorithm.[30] Commission for Energy Regulation, Dublin. "Results of Electricity Coast-Benefit Analysis, Customer Behaviour Trials and Technology Trials Commission for Energy Regulation." In: (2011)

	Cited on pages 9,
	75.

set_missing_parameters_to_default •

Pour assurer la sécurité et la stabilité des réseaux électriques, l'électricité injectée dans les réseaux doit correspondre à la demande à tout moment. C'est la responsabilité des entreprises d'électricité comme Électricité de France (EDF) de gérer cette équilibre offre-demande. La demande d'électricité dépend essentiellement du portefeuille de clients et des conditions extérieures, comme par exemple la météorologie, sur lesquelles les entreprises d'énergie ont peu de marge de manoeuvre. Or, l'électricité est une forme d'énergie qui peut difficilement se stocker. Ainsi les entreprises d'électricité gèrent l'équilibre offre-demande en ajustant la production d'électricité et la configuration des réseaux, selon les estimations et prédictions de demande en temps réel.On peut donc comprendre pourquoi les données de consommation électrique sont importantes pour une entreprise comme EDF. Pour la planification de production nucléaire et la planification des réseaux, il est important d'avoir des prévisions de demande de moyen ou long terme (à partir de 6 mois). Pour la planification des opérations réseaux, il est important d'avoir des prévisions de moyen ou court terme, par exemple des prévisions horaires ou quotidiennes avec un horizon entre un jour et quelques mois. Pour profiter de la production décentralisée d'énergie en plein développement (par exemple photovoltaïque), il est intéressant d'avoir des données ou prévisions sur une échelle spatiale assez fine.Outre ces problèmes dans le domaine du management d'énergie, les données de consommation électrique sont utiles pour beaucoup d'autres applications. Dans le contexte de marché ouvert

On entend par "un individu", l'unité basique considérée dans l'application. Ceci peut être une personne physique, mais aussi un compteur électrique, ou encore un capteur sur un niveau du réseau électrique.

Historically, this is mostly done by hydroelectric energy storage, namely pumping water into a reservoir in higher locations during low-demand periods, and using that water to generate electricity during peaks. Nowadays, battery storage is extensively researched, but remains expensive and experimental for the moment.

Jérémy Lesuffleur, personnal communication, 2017.

We will call a basic unit an "individual", if it is the smallest unit considered in a given context, whether it is actually a physical person, a meter, a household or a group of clients.

SIRENE, Systeme Informatique pour le Repertoire des Entreprises et de leurs Etablissements, the INSEE database of all companies, public and private organization.

The scales are relative on Figs.4.3 and 4.4, since the estimation is applied to normalized time series.

Results for other hours are qualitatively similar.

This intuition is confirmed by the comparison of Algorithm 1 and our implementation of the algorithm proposed in[79].

The distributed computation in model estimation is implemented by Gustavo Castro during his internship at EDF Lab Paris-Saclay.

Chapter

4

Spatial estimation of electricity consumption using sociodemographic information

Contents SOCIO-DEMOGRAPHIC INFORMATION

In our tests, hourly consumption data is used. To account for intra-day variations, we model the 24 hours separately, resulting in 24 consumption models.

Computationally, the parametric specification is fitted using classical least square square method (function lm in the statistical programming language R). The estimation of a GAM model is handled in the mgcv package in R [START_REF] Wood | Generalized additive models: an introduction with R[END_REF].

Cluster layer

We use the standard k-means algorithm in data mining to form clusters ([77, Chap 14]) of geographical zones. Given the fixed number of clusters, and characteristic variables for individuals, this algorithm alternatively iterate between:

• calculating the centroids of the variables of given clusters;

• placing individuals into clusters with the closest centroid.

Starting from a random partition, k-means results in a partition of the individuals into clusters of similar characteristics.

We form clusters of geographical zones. Therefore, the "individuals" are geographical zones. Because source zones and target zones have fundamentally different status in our methodology, we use a two-step procedure: the k-means algorithm is run on the source zones to form clusters; the target zones are then distributed into the cluster with the closest centroid. This assures that in each cluster, there are always some source zones. We consider the following two types of characteristics to form clusters:

• socio-demographic variables (IRIS): the percentages of population in each age group, the number of office buildings, shops, and factories in the zone, obtained from public census data;

• client information (CP): the percentages of residential, business, corporate, or industrial clients in contractual power (the subscribed maximal power), obtained from grid operators' knowledge of their clients.

The two groups of variables result in two clustering, called IRIS and CP clustering respectively.

For each cluster, we estimate one common consumption model using data of all included source zones. This model is then transported to the target zones of the same cluster: we use the regression model to estimate their load, by plugging in temperature and calendar information of the periods of interest.

Since consumption models are sensitive to the relative scaling between source zones of one cluster, the estimation procedure must be fed with normalized load time series. When calendar information and temperature is plugged into the model for prediction and estimation, the level of the consumption has to be supplied separately. In the application, we either predict the level by a linear regression model with CP variables as explanatory variables, or simply use the observed level of target zones.

Application on real and synthetic load data

In the application, we consider a region around Lyon in France. The source zones are geographical zones served by medium-voltage (MV) feeders. Two kinds of target zones are considered: On the samples with random observation periods, proposed methods (HALS and NeNMF, blue and purple lines, both penalized and unpenalized) out-performs the interpolation benchmark (solid red line). On the samples with regular observation periods, unpenalized HALS and NeNMF (solid blue and purple lines) are similar to the interpolation benchmark, whiled penalized HALS and NeNMF (dashed blue and purple lines) are an important improvement. The softImpute method (solid green line) only has comparable performance in two of the datasets, in the easiest task (50% sampling rate at random periods). In most cases, RRMSE of softImpute is larger than 100%.

Part II

Nonnegative matrix factorization with general linear measurements and its applications in time series recovery and prediction

Perspectives

Motivated by a new industrial application, we extended NMF to use temporal aggregates as input data, by adding a projection step into NMF algorithms. With appropriate projection algorithms, this approach could be further generalized to other types of data, such as disaggregating spatially aggregated data, or general linear measures. When such information is available, we introduce a penalization on individual autocorrelation, which improves the recovery performance of the base algorithm. This component can be generalized to larger lags (with a matrix ∆ with 1's further off the diagonal), or multiple lags (by adding several lag matrices together). It is also possible to generalize this approach to other types of expert knowledge, such as exogenous variables, which is discussed in the next chapter. The rank chosen by the cross validation procedure generally increases with the sampling rate, for the four datasets. This shows that the procedure is able to relax the rank constraint when more information is available in the data. INFORMATION In this chapter, we generalize this idea to nonnegative matrix factorization. Although more difficult to identify, NMF often has a better empirical performance, and provides factors that are more interpretable in an appropriate context. We propose an NMF method that takes into account side information. Given some observations of a matrix, this method jointly estimates the nonnegative factors of the matrix, and regression models of these factors on the side information. This allows us to improve the matrix recovery performance of NMF. Moreover, using the regression models, we can predict the value of interest for new rows and columns that are previously unseen. We develop this method in the general matrix recovery context, where linear measurements are observed instead of matrix entries. This choice is again motivated by applications in electricity consumption. We are interested in estimating and predicting the electricity load from temporal aggregates. In the context of load balancing in the power market, electric transmission system operators (TSO) of the electricity network are typically legally bound to estimate the electricity consumption and production at a small temporal scale (half-hourly or hourly), for market participants within their perimeter, i.e. utility providers, traders, large consumers, groups of consumers, etc. [2]. Most traditional electricity meters do not provide information at such a fine temporal scale. Although smart meters can record consumption locally every minute, the usage of such data can be extremely constrained for TSOs, because of the high cost of data transmission and processing and/or privacy issues. Nowadays, TSOs often use regulator-approved proportional consumption profiles to estimate market participants' load. In Chapter 6, we proposed to solve the estimation problem by NMF using temporal aggregate data.

Using the method developed in this chapter, we put in parallel temporal aggregate data with features that are known to have a correlation with electricity consumption, such as the temperature, the day of the week, or the type of client. This not only improves the performance of load estimation, but also allows us to predict future load for users in the dataset, and estimate and predict the consumption of new users previously unseen. In electrical power networks, load prediction for new periods is useful for balancing offer-demand on the network, and prediction for new individuals is useful for network planning. Moreover, by examining the relationship between external features, the factors produced by NMF can be much more interpretable.

In the rest of this section, we introduce the general framework of this method, and the related literature. In Section 7.2, we deduce a sufficient condition on the side information for the NMF to be unique. In Section 7.3, we present HALSX, an algorithm which solves the NMF with side information problem, that we prove to converge to stationary points. In Section 6.3, we present experimental results, applying HALSX on simulated and real datasets, both for the electricity consumption application, and for a standard task in collaborative filtering.

General model definition

In this chapter, we use general linear measurements. We do not restrict us to temporally aggregated observations. As mentioned in Chapter 5, we are interested in reconstructing a nonnegative matrix

where A : R n 1 ×n 2 → R N is a linear operator. Formally, A can be represented by A 1 , ..., A N , N design matrices of dimension n 1 × n 2 , and each linear measurement can be represented by

2)

The design matrices A 1 , ..., A N are called masks. INFORMATION

• Rank-one projections [84,85]: the design matrices are random rank-one matrices, that is A i = α i β T i , where α i and β i are respectively random vectors of dimension n 1 and n 2 . The main advantage to this setting is that much less memory is needed to store the masks, since we can store the vectors α i and

. In [84,85], theoretical properties are proved for the case where α i and β i are vectors with independent Gaussian entries and/or drawn uniformly from the vectors of the canonical basis.

• Temporal aggregate measurements: in this case, the matrix is composed of n 1 time series concerning n 2 individuals, and each measure is a temporal aggregate of the time series of an individual. The design matrices are defined as

t=t 0 (i)+1 e t e T s i , where s i is the individual concerned by the i-th measure, t 0 (i) + 1 the first period covered by the measure, and h(i) the number of periods covered by the measure.

The features X r and X c

• Individual features: X r = I n 1 , X c = I n 2 . Basically, no side information is available.

The row individuals and column individuals are each different.

• General numeric features: X r ∈ R n 1 ×d 1 and X c ∈ R n 2 ×d 2 . This includes all numeric features.

• Features generated from a kernel: certain information about the row and column individuals may not be in the form of a numeric vector. For example, if the row individuals are vertices of a graph, their connection to each other is interesting information for the problem, but it is difficult to encode as real vectors. In this case, features can be generated through a transformation, or by defining a kernel function.

The link functions f r and f c

• Linear: F r = f r (X r ) = X r B r , and

In this case, we need to estimate B r and B c to fit the model. With identity matrices as row and column features, this case is reduced to the traditional matrix factorization model with

When the features are generated from a kernel function, even a linear link function permits non-linear relationship between the features and the factor matrices.

• General regression models: when the relationship between the features and the variable of interest is not linear, any off-the-shelf regression methods can be plugged in to search for a non-linear link function.

The choice of the optimization problem (7.3)

Notice that with individual row and column features, linear link functions and complete observations, (7.3) becomes min

This is equivalent to the classical NMF problem, min

Prior works

Table 7.1 shows a taxonomy of matrix factorization models with side information, by the mask, the link function and the features used as side information.

There is an abundant literature that studies the general matrix factorization problem under various measurement operators, when no additional information is provided (see [78,85,86,104] for various masks considered, and [90] for a recent global convergence result with RIP measurements). The NMF with general linear measurements is studied in various applications [2,79,83].

On the other hand, with complete observations, the multiple regression problem taking into account the low-rank structure of the (multi-dimensional) variable of interest is known as reduced-rank regression. This approach was first developed very early (see [START_REF] Velu | Multivariate Reduced-Rank Regression: Theory and Applications[END_REF] for a review). Recent developments on rank selection [START_REF] Bunea | Joint Variable and Rank Selection for Parsimonious Estimation of High-Dimensional Matrices[END_REF], adaptive estimation procedures [START_REF] Chen | Reduced Rank Regression via Adaptive Nuclear Norm Penalization[END_REF], using nonparametric link function [START_REF] Foygel | Nonparametric Reduced Rank Regression[END_REF], often draw the parallel between reduced-rank regression and the matrix completion problem. However, measurement operators other than complete observations or matrix completion are rarely considered in this community.

Building on theoretical boundaries on matrix completion, the authors of [START_REF] Jain | Provable Inductive Matrix Completion[END_REF][START_REF] Xu | Speedup Matrix Completion with Side Information: Application to Multi-Label Learning[END_REF][START_REF] Chiang | Matrix Completion with Noisy Side Information[END_REF] showed that by providing side information (the matrix X), the number of measurements needed for exact matrix completion can be reduced. Moreover, the number of measurements necessary for successful matrix completion can be quantified by measuring the quality of the side information [START_REF] Chiang | Matrix Completion with Noisy Side Information[END_REF].

Collaborative filtering with side information has received much attention from practitioners and academic researchers alike, for its huge impact in e-commerce and web applications [START_REF] Jain | Provable Inductive Matrix Completion[END_REF][START_REF] Rao | Collaborative Filtering with Graph Information: Consistency and Scalable Methods[END_REF][START_REF] Si | Goal-Directed Inductive Matrix Completion[END_REF][START_REF] Agarwal | Regression-Based Latent Factor Models[END_REF][START_REF] Xu | Speedup Matrix Completion with Side Information: Application to Multi-Label Learning[END_REF][START_REF] Chiang | Matrix Completion with Noisy Side Information[END_REF]. One of the first methods for including side information in collaborative filtering systems (matrix completion masks) was proposed by [START_REF] Abernethy | A New Approach to Collaborative Filtering: Operator Estimation with Spectral Regularization[END_REF]. The authors generalized collaborative filtering into a operator estimation problem. This method allows more general feature spaces than a numerical matrix, by applying a kernel function to side information. [START_REF] Si | Goal-Directed Inductive Matrix Completion[END_REF] proposed choosing the kernel function based on the goal of the application. [START_REF] Kekatos | Electricity Market Forecasting via Low-Rank Multi-Kernel Learning[END_REF] applied the kernelbased collaborative filtering framework to electricity price forecasting. Their kernel choice is determined by multi-kernel learning methods.

To the best of our knowledge, matrix factorization (nonnegative or not) with side information, from general linear measurements has rarely been considered, nor is general non-linear functions other than with features obtained from kernels. This chapter aims at proposing a general approach which fills this gap.

Identifiability of nonnegative matrix factorization with side information

Matrix factorization is not a well-identified problem: for one pair of factors (F r , F c ), with V * = F r F T c , any invertible matrix R produces another pair of factors, since (F r R)(F c (R -1 ) T ) T is also equal to V * . In order to address this identifiability problem, one has to introduce extra constraints on the factors.

When the nonnegativity constraint is imposed on F r and F c , however, it has been shown that sometimes the only invertible matrices that verify F r R ≥ 0 and R -1 F c ≥ 0 are the composition of a permutation matrix and a diagonal matrix with strictly positive diagonal elements. A nonnegative matrix factorization is said to be "identified" if the factors are unique up to permutation and scaling. The identifiability conditions for NMF are a hard problem, INFORMATION because it turns out that NMF identifiability is equivalent to conditions that are computationally difficult to check. In this section, we review some known necessary and sufficient conditions for NMF identifiability in the literature, and develop a sufficient condition for NMF identifiability in the context of linear numerical features.

The reason for studying identifaibility of NMF with side information is two-fold. First of all, it is possible that using side information could restrict the solution space of the NMF problem, and makes the factorization problem more identifiable. Secondly, if we wish to find interpretation of the factors obtained in NMF, it is desirable if the factors are "unique".

In order to simplify our theoretical analysis, we focus on the complete observation case in this section (every entry in V * is observed). Without loss of generality, we derive the sufficient condition for row features. That is, we will derive conditions on V * ∈ R n 1 ×n 2 + and X r ∈ R n 1 ×d 1 , so that the nonnegative matrix factorization V * = X r B r F T c , with X r B r ≥ 0, F c ≥ 0, is unique. A generalization to column features can be easily obtained. In this section, we assume that in addition to be of nonnegative rank k, matrix V * is also exactly of rank k.

Identifiability of NMF

The authors of [5] and [12] proposed two necessary and sufficient conditions for the factorization to be unique. Both conditions use the following geometric interpretation of NMF introduced by [5].

Since V * = F r F T c , the columns of V * are conical combinations of the columns in F r . Formally, cone(F r ), the conical hull of the columns of F r , is a polyhedral cone contained in the first orthant of R n 1 . As V * is of rank k, the rank of F r is also k. This implies that the extreme rays (also called generators) of cone(F r ) are exactly the columns of F r , which are linearly independent. cone(F r ) is therefore

• containing all columns of V * .

Inversely, if we take any cone F ⊆ R n 1 verifying these three conditions, and define a matrix F whose columns are the k generators of F, there will be a nonnegative matrix G, so that V * = FG. The uniqueness of NMF is therefore equivalent to the uniqueness of simplicial cones of k generators contained in the first orthant of R n 1 and containing all columns V * .

In [12], an equivalent geometric interpretation in R k is given in the following theorem:

Despite the apparent simplicity of the theorem, the necessary and sufficient conditions are very difficult to check. Based on the theorem above, several sufficient conditions have been proposed. The most widely used condition is called the separability condition. Before introducing this condition (in its least restrictive version presented by [12]), we need the following two definitions. INFORMATION Definition 1 (Separability). Suppose m ≤ n. A nonnegative matrix M ∈ R m×n + is said to be separable if there is a m-by-m permutation matrix Π which verifies

For a given full-rank matrix X ∈ R n 1 ×d 1 , consider the following two sets of matrices:

As a conclusion to this section, we notice that since the sufficient condition proposed by Theorem 5 is based on a known uniqueness condition of classical NMF, this does not make the problem more identifiable. Rather, we derived a sufficient condition on the feature matrix and the coefficient matrix (X r and B r ) for the NMF with side information to be unique. This is a generalization of identifiability analysis of NMF to the case with side information.

HALSX algorithm

In this section, we propose HALSX, or Hierarchical Alternating Least Squares with eXogeneous variables, a general algorithm to estimate the nonnegative matrix factorization problem with side information, from linear measurement, by solving (7.3). It is an extension to a popular NMF algorithm: Hierarchical Alternating Least Squares (HALS) (see [18,19]).

Before discussing HALSX, we will first present a result on the local convergence of Gauss-Seidel algorithms. This result guarantees that any legitimate limiting points generated by HALSX are stationary points of (7.3).

While presenting specific methods to estimate link functions, we will only discuss row features, as a generalization to column features is immediate.

Relaxation of convexity assumption for the convergence of Gauss-Seidel algorithm

To show that all legitimate limiting points of HALSX are stationary points, we first extend a classical result concerning block nonlinear Gauss-Seidel algorithm [110,Proposition 4].

Consider the minimization problem, min g(x)

where g is a continuously differentiable real-valued function, and the feasible set X is the Cartesian product of closed, nonempty and convex subsets

Suppose that the global minimum is reached at a point in X. The m-block Gauss-Seidel algorithm is defined as Algorithm 3.

Define formally the notion of component-wise quasi-convexity.

This segment has strictly positive dimension in the subspace corresponding to X i .

By the definition of {v t }, g(v t ) ≤ g(y t 1 , ..., z i , ..., y t m ), for all t, and for all z i ∈ X i . In particular,

By quasi-convexity of g on X, g(y t + τ 0 s t ) ≤ max{g(y t ), g(v t )} = g(y t ).

Taking the limit when t converges to +∞ on both equalities, we obtain

In other words, g(ȳ

HALSX algorithm

To solve (7.3), we propose HALSX (Algorithm 4). When complete observations are available, the feature matrices are identity matrices, and when only linear functions are allowed as link functions, Algorithm 4 is equivalent to HALS [19].

From Theorem 6, one deduces that every full-rank limiting point produced by the popular HALS algorithm is a critical point.

In this algorithm, at each elementary update step, we first look for a link function which minimizes the quadratic error, without concerning ourselves with its nonnegativity. The obtained evaluation of the minimizer function is then thresholded at 0 to update the factors.

To obtain that the limiting points of HALSX are stationary points, we need to ensure that for some functional spaces F r and F c , such an update solves a corresponding subproblem of (7.3). To do this, we will use the following proposition:

+ are not identically equal to zero, and

If ∇g θ * , the Jacobian matrix of g at θ * , is of rank n 1 , then θ * is also a solution to

+ not identically equal to zero. We will note by L the loss function, so that

The function L is convex. Problem (7.11), which can be rewritten as min

Algorithm 4 Hierarchical Alternating Least Squares with eXogeneous variables for NMF (HALSX)

Require: Measurement operator A, measurements b, features X r and X c , functional spaces F r and F c in which to search the link functions, and 1

Stopping criterion

As in classical NMF algorithms, we will use the Karush-Kuhn-Tucker conditions (KKT) to provide a stopping criterion. The KKT conditions of (7.13) are,

where A • B is the entry-wise product (Hadamard product) for A, B of the same dimension, and ∇ x f (x 0 ) is the subgradient of the function f at point x 0 , with respect to the variable x. Note that V ≥ 0 and A(V) = b are always satisfied at the end of an iteration.

) are respectively in the subgradient with respect to F c and B r , we will stop the algorithm when the norm of following vector

is smaller than times its initial value, with a small . For the algorithms presented in the next sections, this stopping criterion is generalized quite easily.

HALSX with smoothing splines

The computation considered above can estimate an NMF with linear features fairly efficiently. However, in real applications, linear link functions are too restrictive. In the following, we will estimate non-linear link functions that are Generalized Additive Models (GAM, [START_REF] Wood | Generalized Additive Models: An Introduction with R[END_REF]).

A Generalized Additive Model is a generalization to Generalized Linear Model (GLM) which includes additive non-linear components. Consider n observations x i , y i , for 1 ≤ i ≤ n, where x i is the vector of features, and y i is an observation of a random variable Y i . Suppose that Y i = µ i + i , where i are independent identically distributed zero-mean Gaussian variables, and µ i = E(Y i ) has the following relationship to the features:

where θ is the vector of parametric model components, g is a known, monotonic, twice-differentiable function, h 1 , h 2 , h 3 , ..., are the non-linear functions to be estimated.

We note by X the matrix grouping the features of all observations. We use penalized regression spline to fit the GAMs. For j = 1, 2, 3, ..., define a spline basis a j = (a j 1 , a j 2 , ...) in which h j , the j-th component of the GAM, is to be estimated. Practically, we search for h j is the L-dimensional vector space

HALSX with other regression models

We can replicate the strategy above to work with other regression models. As mentioned before, the convergence to critical point is guaranteed, as long as the regression model estimation can be re-parameterized to verify the conditions of Proposition 1. Using this strategy, many off-theshelf algorithms for regression model training can be plugged in. In the experiments described in the next section, we use the predictive model API provided in the R package caret [START_REF] Kuhn | Building Predictive Models in R Using the Caret Package[END_REF].

Meta-parameters in regression models

As in HALSX-GAM, we build the estimation of meta-parameters using cross validation as a part of the link function estimation step, can treat them indifferently as regular parameters.

An HALS-like algorithm for (7.8)

Before detailing Algorithm 5 which aims to solve (7.8), we will first develop the elemental HALS iteration in the context of (7.8) where no supplemental information is supplied for the factorization model, namely X r = I n 1 , X c = I n 2 . Indeed, when updating one column of F r , the sub-problem becomes: how to solve arg min

and for all b ∈ R N , A * , the transpose of A is defined by

The left-hand side of the equation can be written as

which leads to the following symmetric n 1 -by-n 1 system on f :

Complexity

At each sub-iteration Algorithm 5, we need to calculate N n 1 -by-n 1 or n 2 -by-n 2 matrices ((A i f c )(A i f c ) T ), then inverse the sum of the these N matrices. While the computation of the sum is map-reducible, on a single-threaded machine, this can be computationally expensive when N is large. Each iteration of Algorithm 5 has a multiplicative complexity of O(kN (n 2 1 + n 2 2 )), while each iteration of Algorithm 4 has a complexity of O(k(n 2 1 +n 2 2 )+N n 1 n 2 ) with general linear measurement operator. This difference in complexity can be very important when N or k is large.

Experiments

Datasets

We use one synthetic dataset and three real datasets to evaluate the proposed methods. • Portuguese electricity consumption [28] daily consumption of 370 Portuguese clients from 2010 to 2014 (T = 1461, N = 369). The daily temperature in Portugal and calendar variables are used as side information.

• MovieLens 100k [START_REF]K Dataset[END_REF] an anonymized public dataset with 100,000 movie scores for 1682 movies from 943 users (T = 943, N = 1682). This is a standard public dataset for matrix completion. Note that the data matrix is not complete. Error rates are calculated on the vector of available scores. The genres of the movies, and gender, age and profession of the users are used as side information.

The subsets of the two real-world electricity consumption datasets are also used in Chapter 6.

Validation procedure

Several methods that can be used in time series recovery and prediction are compared to the HALSX algorithm proposed in this chapter.

Among these methods, the following methods (introduced in Chapter 6) are used to compare the matrix recovery/completion performance.

• empty_model For temporal aggregates only: to recover the target matrix, temporal aggregates are distributed equally over the covered periods. INFORMATION • for MovieLens, we calculate the l 2 -norm version on the vector of all available movie scores.

In preliminary tests, other error metrics were also evaluated, and were not qualitatively different from the chosen metric.

Execution time and precision between HALS and HALS2

To show-case the computational complexity difference discussed in Section 7.3.7, we run HALS (Algorithm 4) and HALS2 (Algorithm 5) on random temporal aggregate masks, with no side information. In HALS, we treated observations as general linear measurements, without using the efficient simplex projection specific to temporal aggregate masks discussed in Chapter 6.

The per iteration execution time is shown in Figure 7.2. The difference in complexity discussed in Section 7.3.7 is fairly clear here. In HALS2, the execution time per iteration increases both with the rank and the number of samples in data, at least for sampling rate from 14.3% on. For low sampling rates, HALS2 often diverges, where as HALS is always stable.

As discussed in Section 7.3.7, although the execution time per iteration is greater, HALS2 could be more efficient if it does much less iterations than HALS. In Figure 7.3, we can see that this is not the case: for problems with high sampling rates, HALS2 indeed does less iterations to converge. However, the total execution time is still larger than HALS. For lower sampling rates, HALS2 has troubles converging, and only stops when the maximal execution time allowed (300 seconds) is reached. This is also confirmed in Figure 7. 4, where HALS2 has much worse recovery error than HALS.

Given these results, we will only use the HALSX (Algorithm 4), the version of HALS with side information, in the following tests.

Performance on temporal aggregate measurements

On the synthetic and the two real electricity consumption datasets, we use Algorithm 4 to perform matrix recovery and prediction on temporal aggregate measurements. Every method INFORMATION We use two types of temporal aggregate measurements: periodic and random. In periodic measurements, each scalar measure covers a fixed number of periods of one individual. In random measurements, the number of periods covered by a measure is random (see Chapter 6 for more details). In electricity consumption data, periodic measurements are closer to the actual meter reading schedules of utility companies, while matrix recovery with random measurements is an easier problem. For both sampling types, we sample 10%, 20%, ..., 50% of the data to show-case the matrix recovery performance of the proposed method.

For the synthetic data, we use the true row and column features used to produce the simulations. For the French electricity data, the row features are variables known to have an influence on electricity consumption: the temperature, the day type (weekday, weekend, or holiday), the position of the year. The column features are the percentage of residential, professional, or industrial usages in the group of users for each column. For the Portuguese electricity data, as no individual features are available, we only use the same row features as for the French dataset (temperature, day type, position of the year). Figure 7.5 shows the matrix recovery error. For most of the scenarios, HALSX_models (red lines with symbols) are comparable or better than the other methods without side information. The only case where an HALSX_model is a little worse is when compared with HALS and NeNMF (two NMF methods) in synthetic data with random measurements, which is the least close to the real application. As is seen in Chapter 6, the softImpute [111] method is not well adapted to temporal aggregate measurements, and has much higher error (higher than the maximal value in these graphics. When comparing the four regression methods used in HALSX, we see that GAM and Gaussian process regression (gam and gaussprRadial) are the best for synthetic data, linear model (lm) is the best for Portuguese data, and Gaussian process regression (gaussprRadial) is the best for French data. methods, except in Synthetic data with complete observations. In most cases, HALSX_models are comparable to or better than factor_gam and individual_gam, which shows that using side information while estimating the factorization model produces factors more adapted for prediction. When comparing the regression methods used with HALSX_models in prediction, gam is consistently the best for synthetic and Portuguese data, gaussprRadial best for the French dataset.

It is interesting to note that the performance HALSX_models is the least sensitive to sampling rates: it is mostly constant from 30% of data. This shows that using side information is supplementary to observing more data.

Performance on matrix completion mask

On the MovieLens dataset, we use Algorithm 4 to perform matrix recovery and prediction with uniformly sampled matrix entries. The sampling rates are {10%, 20%, 30%, 40%, 50%, 90%}. As is the case for time series datasets, we use samples from the upper-left submatrix to estimate the model, evaluate matrix recovery on that submatrix, and evaluate row and/or column prediction errors. Every method except for trmf is used in this setting. As side information, we use the gender, the age, and the occupation of the users and the genre (a dimension-19 binary variable) of the movies.

For grmf, we produce a graph where each individual is connected to its ten nearest neighbors with euclidean distance with the features, as suggested in the reference [START_REF] Rao | Collaborative Filtering with Graph Information: Consistency and Scalable Methods[END_REF]. As the parameters estimated in grmf is per user/movie, it can not be used to predict new individuals, even though it uses side information. At this point of the thesis, we have seen both several methods that are derived from different contexts: spatial statistics (Chapter 3), socio-demographic clustering (Chapter 4), and NMF (Chapters 5, 6, 7). In this chapter, we discuss the link between the main proposition of this thesis, using nonnegative matrix factorization to recover and predict partially observed time series, with the two state-of-art methods investigated in Chapters 3 and 4.

In 8.1, we apply the proposed general matrix factorization framework to spatial and spatiotemporal kriging problems. In Section 8.2, we argue that the method proposed in Chapter 4 can be seen as a special case of matrix factorization with exogenous variables.

Kriging and matrix factorization

Consider the following matrix factorization model: min

are the functional spaces in which we try to find the link functions, X r and X c are known row and column features, A a linear measurement operator, and α is the measurements that we obtained by applying A on the ground truth V * . For the moment, we do not constrain the factors to be nonnegative.

In particular, we will consider the column features X c to be spatial coordinates of the column individuals, and the row features X r to be temporal coordinates (timestamps).

Kriging as matrix factorization with a feature map

The link between Gaussian process regression and kernel regression has been extensively studied (see [44,Chap. 2]). In this section, we first examine a toy model, to see that univariate kriging is a special case of (8.1).

Appendices

Appendix

A

Reference manual of package 'meterModels' Value a named list. Each element of the list is a named vector of the evaluation metrics for one of the prediction_type.

See Also

implemented_prediction_types and implemented_matrix_comparison_metrics for details of the different types of evaluation value possible.

Examples

fac_mod = factorization_model(example_index_data, r = 2, algorithm = "HALS") # For a factorization_model, regression and prediction errors are NA evaluate_matrix_model(fac_mod, example_incremental_data) exp_spec See Also formula_params for specification of the regression formulae, experiment for using an exp_data in action.

Examples

# create exp_data with incremental_data example_exp_data = exp_data(example_incremental_data, m_rho = 0, m_train = 2, n_train = 2, formula = example_formula)

# create exp_data with linear_measurement_data example_exp_data = exp_data(example_linear_measurement_data, m_rho = 0, m_train = 2, n_train = 2, formula = example_formula) exp_spec

Create an object specifying parameters of an experiment

Description

Create an object specifying parameters of an experiment Usage exp_spec(r, algorithm = implemented_algorithms("all"), random_start = 3, lambda = c(0, 1, "optimal"), general_factorization = c(TRUE, FALSE), parallel_opts = list(mc.cores = 1, mc.preschedule = FALSE), opts = list(), save_model = file.path(getwd(), "exp_tmp"), eval_only = FALSE) is.exp_spec(x)

Arguments r, algorithm, lambda, general_factorization

The same parameters as in factorization_model, except they can be vectors so that a combination of all parameters are tested.

random_start an integer to specify the number of random runs to turn for each experiment spec. When using parallel to parallelize, these random runs with one single experiment spec will be run on the same core. opts a named list to set the parameters for the estimation algorithm with the following fields:

• tol: tolerance for deciding whether to stop.

• maxit: the maximal number of iterations.

• minit: the minimal number of iterations.

• verbose: whether to provide more reporting in the computation or not.

• fixed_X, fixed_Y: whether X or Y should be fixed in the iterations.

• maxtime: maximal computing time allocated in seconds.

• ITER_MAX, ITER_MIN: the maximal and minimal number of iterations for NeNMF. • X0,Y0: initial values for X, Y implemented_prediction_types • "partial_data": algorithms that can be used with partial data.

• "completion": algorithms that can be used to do matrix completion.

• "factorized": all matrix factorization algorithms, including regression-enhenced ones.

• "linear_measurement_data": algorithms that can be used with objects of the class linear_measurement_data.

• "index_data": algorithms that can be used with objects of the class index_data.

• "proposed": algorithms proposed by the thesis.

• "reference": everything except "proposed"

• "all": everything.

When multiple parameters are put in, the intersection of the sets are returned. Default is "all".

Examples implemented_algorithms() implemented_algorithms("reference","factorized")

implemented_prediction_types

The prediction types

Description

The prediction types Usage implemented_prediction_types()

Details

For object of class matrix_model, whenever applicable, the four implemented types are

• "factorization": object$X %*% object$Y.

• "completion": object$X %*% object$Y projected into the subspace satisfying the data constraint.

• "regression": predict(object$models$X) %*% predict(object$models$Y).

• "factorization": same as "regression", but with new data.

Examples implemented_prediction_types() implemented_sampling_type

The implemented sampling types

Description

Lists the implemented sampling types for each type of data Usage implemented_sampling_type(type = "index_data") Arguments type either "index_data" or "linear_measurement_data".

Value

a vector of the sampling types.

Examples implemented_sampling_type("index_data")

increment_evaluate Evaluates a matrix estimation against the data matrix.

Description

Evaluates a matrix estimation against the data matrix. theta the smoothing parameter. Should be between 0 and 1.

Value

A numeric vector of the same length as vect.

match_aggregate_max References the package R39Toolbox developped by R39.

Examples

lissageV(1:10, 0.1) match_aggregate_max

Evaluation functions concerning aggregated (daily) peaks

Description

These functions compare the periodic peaks (for example daily or hourly or weekly) of the individual aggregated ground truth with that of an estimation. This is to evalute if the NMF algorithms can be applied to estimate periodic peaks from temporally aggregated individual data.

Usage

diff_aggregate_max_positive(matrix1, matrix2, timestep, aggregate_vector) diff_aggregate_max_negative(matrix1, matrix2, timestep, aggregate_vector)

Arguments matrix1, matrix2

The ground truth and the estimation. timestep

The temporal aggregation step at which to calculate the peaks. If the each row of matrix1 and matrix2 are half hourly datat, and we are interested in daily peaks, timestep should be 48. aggregate_vector

The aggregation vector to aggregate individuals.

Details

match_aggregate_max calculate the percentage of peak moments that are correctly estimated.

diff_puissance_moment_pmax calculates the average deviation between the estimation and ground truth at the peak moments ( both estimated and real). matrix_model diff_aggregate_max calculates the average deviation between the estimaed peaks and the real peaks.

diff_aggregate_max_positive calculates the average positive deviation between the estimaed peaks and the real peaks.

diff_aggregate_max_negative calculates the average negative deviation between the estimaed peaks and the real peaks.

Value

a single value representing an evaluation metric.

matrix_model

General matrix model

Description

Create a general matrxi model. An object of matrix_model class has a predict method.

There are two types of matrix models, factorization_model and regression_model. The difference between them is that predict.regression_model allows new data. This is where the dispatching between factorization and regression models happens based on the data and the algorithm. ... other parameters to be passed on to factorization_model or regression_model.

Value

an object of class matrix_model which has at least the following fields:

• M: the recovered matrix.

• algorithm: the algorithm.

• obs: a named list of the observations provided for model estimation.

• args_list: a named list of the other arguments.

If algorithm is not compatible with the data and other parameters that are provided, this function will return an empty_model, which is a naive model. parameters to be used to pass one representation of the measurement operator. See Details for more explanation.

Details

The parameter type can be of the following values:

• "completion" a binary matrix of the same dimension as the matrices to be measured.

The "pattern" matrix should be supplied under the parameter name pattern_mat. • "matrix" a matrix which is the measurement operator in the matrix form. The parameter matrix_A should be supplied. It should be a sparseMatrix (class from the Matrix package). • "list_mat" a list of matrices, each of the same dimension as the matrices to be measured. Each matrix is a "mask". The parameter list_mat should be supplied. • "list_sub", "list_ind", "vec_ind". All three are used for completion data. "list_ind" is a list of length of nrow of the matrix to be measured, in which each element is a vector of the row indexes of the observations. "list_sub" is a list of vectors of the "sub" (row column pair) of the observations (the parameter name is list_subs). "vec_ind" is a vector of the indexes (by seeing the matrix as a very long vector) of the observations. opts a named list to set the parameters for the estimation algorithm with the following fields:

• tol: tolerance for deciding whether to stop.

• maxit: the maximal number of iterations.

• minit: the minimal number of iterations.

• verbose: whether to provide more reporting in the computation or not.

Details

When not specified explicitly, predict treats a regression_model that is also a factoriza-tion_model as a regression_model.

Value

If newdata is supplied and object is a regression_model with appropriate $algorithm, a matrix of dimension nrow(newdata$X)-by-nrow(newdata$Y). If only one of the newdata data.frames is provided, the missing dimension is replaced by that of the data used to estimate object.

If no newdata is supplied, the matrix recovery result is returned, either projected or not, depending on the value of completion.

If the parameters are not compatible, NULL is returned.

print.formula_params

See Also

For examples, see matrix_model, regression_model, factorization_model, and empty_model. A named list of parameters for the "HALS_model" algorithm. Should include • reg_model: a character string specifying which regression method to use. Supported values are "gam" or regression methods supported by caret. If is "gam", the mgcv package will be used. If is anything else, the caret package will be used for estimation of the regression models. • Other parameters are passed in the estimation iterations to either mgcv::gam or caret::train.

HALSX1_params

A named list of parameters for the "HALSX1" algorithm. Can use this to pass on a parameter called sp to specify penalization parameter in spline regression. trmf_params A named list of parameters for the "trmf" algorithm. The parameters are optional:

• lag_idx: The number of periods considered in the autoregressive models of the factors. • maxit: The maximal number of iterations. • lambdas: A vector of three numbers: penalization parameters on X, Y and M (the matrix to be recovered).

Details

A formula_params object has a formula_alone which is a string that looks like the righthand side of a typical R formula. Other fields are used to store algorithm-specific parameters, so HALS_model_params is for the "HALS_model" algorithm etc..

Value

a formula_params object. Basically a wrapper of the parameters.

https://topepo.github.io/caret/available-models.html for the list of available regression methods in caret.

Examples

formula_params(create_formula_string(colnames( example_exogenous_variables$X), type = "gam"), HALS_model_params = list(reg_method = "gam")) formula_params(create_formula_string(colnames( example_exogenous_variables$X), type = "lm"), HALS_model_params = list(reg_method = "xgbTree", # All tuning parameters to pass to # caret is to be supplied here subsample = 0. algorithm a character string with the name of the algorithm to be used. See implemented_algorithms for the complete list. opts a named list to set the parameters for the estimation algorithm with the following fields:

• tol: tolerance for deciding whether to stop.

• maxit: the maximal number of iterations.

• minit: the minimal number of iterations.

• verbose: whether to provide more reporting in the computation or not.

• fixed_X, fixed_Y: whether X or Y should be fixed in the iterations.

• maxtime: maximal computing time allocated in seconds.

• ITER_MAX, ITER_MIN: the maximal and minimal number of iterations for NeNMF. • X0,Y0: initial values for X, Y • FuninR: If FuninR = "R" use the R version of subroutines. If FuninR = "cpp" use the C++ version. Default is "cpp". • ccdppitmax,ccdppthreads,ccdpplambda: the maximal number of iterations, the number of threads, and the penalization parameter for CCDPP algorithm.

Arguments sampling_rate a number strictly between 0 and 1, or a numeric vector. If type == "evolution", it is two numbers between 0 and 1. If type == "all", this is ignored. • tol: tolerance for deciding whether to stop.

• maxit: the maximal number of iterations.

• minit: the minimal number of iterations.

• verbose: whether to provide more reporting in the computation or not.

• fixed_X, fixed_Y: whether X or Y should be fixed in the iterations.

Time series recovery and prediction with regression-enhanced nonnegative matrix factorization applied to electricity consumption

Keywords Spatial analysis, Time series, Electricity consumption, Source separation, Nonnegative matrix factorization Abstract We are interested in the recovery and prediction of multiple time series from partially observed and/or aggregate data. Motivated by applications in electricity network management, we investigate tools from multiple fields that are able to deal with such data issues. After examining kriging from spatio-temporal statistics and a hybrid method based on the clustering of individuals, we propose a general framework based on nonnegative matrix factorization. This framework takes advantage of the intrinsic correlation between the multivariate time series to greatly reduce the dimension of the parameter space. Once the estimation problem is formalized in the nonnegative matrix factorization framework, two extensions are proposed to improve the standard approach.

The first extension takes into account the individual temporal autocorrelation of each of the time series. This increases the precision of the time series recovery. The second extension adds a regression layer into nonnegative matrix factorization. This allows exogenous variables that are known to be linked with electricity consumption to be used in estimation, hence makes the factors obtained by the method to be more interpretable, and also increases the recovery precision. Moreover, this method makes the method applicable to prediction. We produce a theoretical analysis on the framework which concerns the identifiability of the model and the convergence of the algorithms that are proposed. The performance of proposed methods to recover and forecast time series is tested on several multivariate electricity consumption datasets at different aggregation level.