diff --git a/01_intro.R b/01_intro.R new file mode 100644 index 0000000..87a0e28 --- /dev/null +++ b/01_intro.R @@ -0,0 +1,50 @@ +# ML 1 01 Intro + +# Exemple de fonction non-linéaire pour générer un jeu de données artificiel. +f <- function(x) {exp(x) * cos(2*pi*sin(pi*x))} + +# Génération d'un jeu de données synthétique qui est l'image par f d'un +# échantillon uniforme de l'intervalle [0,1] auquel nous ajoutons un bruit +# gaussien de moyenne nulle et d'écart type sd. +gendat <- function(N, sd) { + # N: nombre d'observations à générer + # sd: écart type d'un bruit gaussien de moyenne nulle + X = runif(N) + Y = f(X) + rnorm(N, mean=0, sd=sd) + dim(X) <- c(N,1) # en général chaque observation est décrite par plusieurs + # variables et X est une matrice avec autant de lignes que + # d'observations et autant de colonnes que de variables. Sur + # notre exemple, chaque observation n'est décrite que par une + # seule variable. + list(X = X, Y = Y) +} + +# Affichage simultanée du jeu de données bruité et de la courbe de la fonction +# utilisée pour le générer. +plt <- function(data, f, ...) { + xs = seq(0,1,length.out=100) + plot(xs, f(xs), type="l", ...) + points(data$X, data$Y) +} + +# Résolution du système linéaire correspondant à la matrice de Vandermonde. +# Autrement dit, découverte d'un polynôme qui passe par chaque point. +polyreg1 <- function(data) { + xs <- c(data$X) # on transforme la matrice X, de dimension Nx1 sur notre + # exemple, en un vecteur + vandermonde = outer(xs, 0:(length(xs)-1), "^") + solve(vandermonde, data$Y) +} + +# Evaluation d'un polynôme en un point. +polyeval <- function(coef,x) { + powers = 0:(length(coef)-1) + f = function(y) { sum(coef * y^powers) } + sapply(x,f) +} + +# Affichage de la courbe d'un polynôme défini par ses coefficients. +pltpoly <- function(coef) { + xs = seq(0,1,length.out=100) + lines(xs, polyeval(coef,xs), lty="dotted") +} \ No newline at end of file diff --git a/01_intro.Rmd b/01_intro.Rmd new file mode 100644 index 0000000..a727a75 --- /dev/null +++ b/01_intro.Rmd @@ -0,0 +1,118 @@ +--- +title: "ML 01 Introduction" +output: + bookdown::pdf_document2: + number_section: yes +toc: false +classoption: fleqn +--- + +```{r, include=FALSE} +source("ML1_01_intro.R", local = knitr::knit_global()) +``` + +# Contexte + +A l’occasion des trois premières révolutions industrielles, des tâches, auparavant réservées au travail manuel de l’Homme, ont été automatisées. Il semble envisageable d’associer au tournant du 21ème siècle une quatrième révolution portée par l’automatisation de la capacité à prédire, essentielle pour le processus de prise de décision dans les secteurs de l’industrie, du commerce et des services. + +Cette transformation se fonde sur des évolutions scientifiques et techniques majeures. Elle est ainsi associée à une discipline, le machine learning ou apprentissage automatique de modèles prédictifs par extrapolation à partir de données générées par des processus physiques, numériques ou biologiques. Ces développements algorithmiques, en particulier la redécouverte des réseaux de neurones profonds, ont révélé sous un nouveau jour leur potentiel autour des années 2010 grâce, d’une part, à la création de jeux de données volumineux dans des domaines variés comme la reconnaissance de la parole, la vision par ordinateur, les données multimedia, le traitement de la langue naturelle, la robotique, les véhicules autonomes… et, grâce d’autre part, à une croissance rapide des capacités de calcul et de stockage aux coûts toujours plus abordables. + +Par ailleurs, cette automatisation des prédictions s’accompagne d’un renouveau des formes de jugement dans les processus de prise de décision avec un couplage de plus en plus fin entre d’un côté, des experts humains et de l’autre, des chaînes de traitement automatique des données qui aboutissent sur la mise en production d’algorithmes prédictifs. Pour la réussite de cette intégration, les compétences de l’ingénieur informaticien sont essentielles. Ce dernier, puisqu’il comprend en profondeur le fonctionnement et les limites des algorithmes qu’il déploie, est capable d’en mesurer les risques et les biais pour éclairer le jugement de ceux qui utiliseront ses réalisations logicielles pour prendre des décisions. + +Ainsi, ce cours introductif à l’apprentissage automatique a pour objectif d’offrir des connaissances fondamentales et des compétences pratiques qui aideront l'ingénieur à tenir ce rôle essentiel. + +# Objectifs + +La discipline de l’apprentissage automatique, ou machine learning, élabore des algorithmes et des méthodes pour découvrir des régularités dans des données multidimensionnelles afin, entre autres, d’automatiser la prédiction. Elle peut se subdiviser en trois catégories. D’abord, l’apprentissage non (ou semi) supervisé qui s’attache à découvrir des structures dans les données non étiquetées à travers des approches comme le clustering, la réduction dimensionnelle, les modèles génératifs… Ensuite, l’apprentissage par renforcement, dans le cadre duquel un agent interagit avec son environnement en adaptant son comportement pour maximiser une fonction de récompense. Enfin, l’apprentissage supervisé, qui fait l’objet de ce module, a quant à lui pour objectif d’apprendre à prédire l’association entre un objet décrit selon plusieurs dimensions et une étiquette. + +Par exemple, il peut s’agir d’associer aux quartiers d’une ville le prix médian d’un logement. Dans ce cas, un quartier peut être décrit par la proportion de zones résidentielles, le taux de criminalité, le nombre moyen de pièces par habitat, etc. Ici, nous faisons face à un problème dit de « régression » où la valeur à prédire, autrement l’étiquette associée à chaque observation, est continue. Lorsque la variable à prédire est discrète, il s’agit d’un problème dit de « classification », comme détecter un objet dans une image ou décider si une transaction bancaire risque d’être une fraude. Nous considérerons ces deux catégories de problèmes. + +Ainsi, ce cours a pour objectif d’introduire quelques concepts fondamentaux de l’apprentissage supervisé et de montrer leurs interconnexions variées dans le cadre de développements algorithmiques qui permettent d’analyser des jeux de données dans une visée avant tout prédictive. Ainsi, les propositions théoriques mèneront à l’écritures de programmes qui implémentent ou utilisent quelques modèles essentiels de l’apprentissage supervisé. + +Pour faciliter l’acquisition des connaissances, le cours est accompagné de notebooks manipulables, rédigés dans le langage de programmation R. Ces mises en pratique systématiques doivent permettre de faire le lien entre des concepts fondamentaux et leur application dans des projets d’analyse de données. + +# Ajustement de courbe + +$\mathbf{x}^{(1)} \dots \mathbf{x}^{(P)}$ sont des vecteurs de $\mathbb{R}^N$ associés aux valeurs, aussi appelées étiquettes, $y^{(1)} \dots y^{(P)}$ de $\mathbb{R}$. Nous cherchons une fonction $f(\mathbf{x}) : \mathbb{R}^N \rightarrow \mathbb{R}$ qui modélise la relation entre les observations $\mathbf{x}$ et les étiquettes $y$. + +La fonction $f$ peut avoir une forme paramétrique, comme par exemple : +\[ +f(\mathbf{x}) = a_0 + a_1x_1 + a_2x_2 + \dots + a_Nx_N +\] +Si $P=N+1$, les paramètres $a_0, a_1, \dots a_N$ sont solution d'un système linéaire : +\[ +\begin{cases} +y^{(1)} &= a_0 + a_1 x_1^{(1)} + a_2 x_2^{(1)} + \dots + a_N x_N^{(1)} \\ +y^{(2)} &= a_0 + a_1 x_1^{(2)} + a_2 x_2^{(2)} + \dots + a_N x_N^{(2)} \\ +\dots &= \dots \\ +y^{(P)} &= a_0 + a_1 x_1^{(P)} + a_2 x_2^{(P)} + \dots + a_N x_N^{(P)} \\ +\end{cases} +\] +Ce système s'écrit également sous forme matricielle : +\[ +\left( \begin{array}{cccc} +1 & x^{(1)}_1 & \dots & x^{(1)}_N \\ +1 & x^{(2)}_1 & \dots & x^{(2)}_N \\ +\dots & \dots & \dots & \dots \\ +1 & x^{(P)}_1 & \dots & x^{(P)}_N +\end{array} \right) +\left( \begin{array}{c} +a_0 \\ a_1 \\ \dots \\ a_N +\end{array} \right) += +\left( \begin{array}{c} +y^{(1)} \\ y^{(2)} \\ \dots \\ y^{(P)} +\end{array} \right) +\] +Chaque ligne $i$ de la matrice du terme de gauche de l'égalité ci-dessus est le vecteur ligne $\mathbf{x}^{(i)T}$ avec l'addition d'un premier terme constant qui correspond au paramètre $a_0$. En nommant cette matrice $\mathbf{X}^T$, le système linéaire ci-dessus s'écrit : +\[ +\mathbf{X}^T \mathbf{a} = \mathbf{y} +\] +Soit le cas particulier où $x$ est un scalaire et $f$ est un polynome de degré $N$ : +\[ +f(x) = a_0 + a_1x + a_2x^2 + \dots + a_Nx^N +\] +Avec $P=N+1$ observations et les étiquettes associées $\left( x^{(k)}, y^{(k)} \right)$, les coefficients de ce polynôme sont solution d'un système linéaire : + +\[ +\left( \begin{array}{ccccc} +1 & x^{(1)} & (x^{(1)})^2 & \dots & (x^{(1)})^N \\ +1 & x^{(2)} & (x^{(2)})^2 & \dots & (x^{(2)})^N \\ +\dots & \dots & \dots & \dots \\ +1 & x^{(P)} & (x^{(P)})^2 & \dots & (x^{(P)})^N +\end{array} \right) +\left( \begin{array}{c} +a_0 \\ a_1 \\ \dots \\ a_N +\end{array} \right) += +\left( \begin{array}{c} +y^{(1)} \\ y^{(2)} \\ \dots \\ y^{(P)} +\end{array} \right) +\] + +La matrice du terme de gauche de l'égalité ci-dessus est traditionnellement appelée "matrice de Vandermonde". + +# Interpolation polynomiale sur un jeu de données synthétique + +Soit un exemple de fonction non-linéaire, $f(x) = e^x \times cos(2 \pi \times sin(\pi x))$, utilisée pour générer un jeu de données synthétique. + +```{r} +set.seed(1123) +# Image par f d'un échantillon uniforme sur l'intervalle [0,1], avec ajout d'un +# bruit gaussien de moyenne nulle et d'écart type 0.2 +data = gendat(10,0.2) +plt(data,f) +``` + +Sur cet exemple, nous résolvons le système linéaire de Vandermonde pour découvrir un polynôme qui passe par chaque point du jeu de données. + +```{r} +# Résolution du système linéaire correspondant à la matrice de Vandermonde. +# coef contient les coefficients d'un polynôme qui passe par chaque point du jeu +# de données. +coef = polyreg1(data) +plt(data,f) +pltpoly(coef) +``` + +Il est improbable que ce polynôme, passant exactement par chaque observation, puisse offrir de bonnes capacités prédictives. Vérifier par exemple que, sur notre exemple synthétique, pour cinq points générés à partir de la fonction $f$ et avec l'ajout d'un bruit gaussien (par exemple d'écart type $0.2$), le polynôme découvert, de degré quatre, peut être très éloigné de la fonction génératrice. C'est un exemple du phénomène de sur-apprentissage. Pour limiter ce problème, nous cherchons à découvrir un polynôme de degré plus faible. Il ne passera pas exactement par toutes les observations, mais il prédira probablement mieux les étiquettes associées à de nouvelles observations. \ No newline at end of file diff --git a/02_moindres_carres.R b/02_moindres_carres.R new file mode 100644 index 0000000..5a16461 --- /dev/null +++ b/02_moindres_carres.R @@ -0,0 +1,8 @@ +# Résolution du système linéaire correspondant à la matrice de Gram pour un +# polynôme de degré fixé. +polyreg2 <- function(data, degre) { + xs <- c(data$X) + A = outer(xs, 0:degre, "^") + gram = t(A) %*% A + solve(gram, as.vector(t(A) %*% data$Y)) +} \ No newline at end of file diff --git a/02_moindres_carres.Rmd b/02_moindres_carres.Rmd new file mode 100644 index 0000000..ce2a8f2 --- /dev/null +++ b/02_moindres_carres.Rmd @@ -0,0 +1,86 @@ +--- +title: "ML 02 Méthode des moindres carrés" +output: + bookdown::pdf_document2: + number_section: yes +toc: false +classoption: fleqn +--- + +```{r, include=FALSE} +source("ML1_01_intro.R", local = knitr::knit_global()) +source("ML1_02_moindres_carres.R", local = knitr::knit_global()) +``` + +# Espace de fonctions + +Soit un espace vectoriel composé de fonctions. Une base de cet espace est un ensemble de fonctions ($f_1, f_2, \dots f_N$) tel que toute fonction de l'espace s'exprime comme combinaison linéaire des fonctions de base. +\[ +f(\mathbf{x}) = a_1 f_1(\mathbf{x}) + a_2 f_2(\mathbf{x}) + \dots + a_N f_N(\mathbf{x}) +\] +Pour un jeu de données $\left\{ \mathbf{x^{(k)}},y^{(k)}\right\}_{k=1}^{N}$ de taille $N$, les coefficients $a_i$ sont solution d'un système linéaire. +\[ +\left( \begin{array}{ccccc} +f_1(\mathbf{x^{(1)}}) & f_2(\mathbf{x^{(1)}}) & \dots & f_N(\mathbf{x^{(1)}}) \\ +f_1(\mathbf{x^{(2)}}) & f_2(\mathbf{x^{(2)}}) & \dots & f_N(\mathbf{x^{(2)}}) \\ +\dots & \dots & \dots & \dots \\ +f_1(\mathbf{x^{(N)}}) & f_2(\mathbf{x^{(N)}}) & \dots & f_N(\mathbf{x^{(N)}}) +\end{array} \right) +\left( \begin{array}{c} +a_1 \\ a_2 \\ \dots \\ a_N +\end{array} \right) += +\left( \begin{array}{c} +y^{(1)} \\ y^{(2)} \\ \dots \\ y^{(N)} +\end{array} \right) +\] +Nous notons ce système linéaire $\mathbf{Ax} = \mathbf{b}$. + +# Expression matricielle + +Le système linéaire $\mathbf{Ax} = \mathbf{b}$ avec $\mathbf{A} \in \mathbb{R}^{M \times N}$ n'a pas de solution quand le nombre d'observations dépasse le nombre de fonctions de base (c'est-à-dire, $M>N$). Une approche possible est alors de chercher une approximation $\mathbf{Ax} \approx \mathbf{b}$ qui minimise la somme des carrés des erreurs : $\|\mathbf{Ax}-\mathbf{b}\|^2_2$. + +\begin{align*} + & \|\mathbf{A}\mathbf{x}-\mathbf{b}\|^2_2 \\ += \{ & \|\mathbf{x}\|_2 = \sqrt{\mathbf{x}\cdot\mathbf{x}} \} \\ + & \left(\mathbf{A}\mathbf{x}-\mathbf{b}\right) \cdot \left(\mathbf{A}\mathbf{x}-\mathbf{b}\right) \\ += \{ & \text{Par définition du produit scalaire euclidien} \} \\ + & \left(\mathbf{A}\mathbf{x}-\mathbf{b}\right)^T \left(\mathbf{A}\mathbf{x}-\mathbf{b}\right) \\ += \{ & \text{propriété de la transposition} \} \\ + & \left(\mathbf{x}^T\mathbf{A}^T - \mathbf{b}^T \right) \left(\mathbf{A}\mathbf{x}-\mathbf{b}\right) \\ += \{ & \text{multiplication} \} \\ + & \mathbf{x}^T\mathbf{A}^T\mathbf{A}\mathbf{x} - \mathbf{x}^T\mathbf{A}^T\mathbf{b} - \mathbf{b}^T\mathbf{A}\mathbf{x} + \mathbf{b}^T\mathbf{b} \\ += \{ & \mathbf{b}^T\mathbf{A}\mathbf{x} \text{ étant une valeur scalaire, } \mathbf{b}^T\mathbf{A}\mathbf{x} = \left(\mathbf{b}^T\mathbf{A}\mathbf{x}\right)^T = \mathbf{x}^T\mathbf{A}^T\mathbf{b} \} \\ + & \mathbf{x}^T\mathbf{A}^T\mathbf{A}\mathbf{x} - 2\mathbf{x}^T\mathbf{A}^T\mathbf{b} + \mathbf{b}^T\mathbf{b} +\end{align*} + +Cette dernière expression quadratique en $\mathbf{x}$ correspond à une surface convexe. Donc son minimum peut être calculé en annulant sa dérivée (penser à une courbe $y = a+bx+cx^2$ dont l'unique extremum est atteint lorsque la pente est nulle). + +\begin{align*} + & \mathbf{0} = 2\mathbf{A}^T\mathbf{A}\mathbf{x} - 2\mathbf{A}^T\mathbf{b} \\ +=& \\ + & \mathbf{A}^T\mathbf{A}\mathbf{x} = \mathbf{A}^T\mathbf{b} +\end{align*} + +Ainsi, quand $M>N$, la solution approximée $\mathbf{x}$, telle que $\mathbf{Ax} \approx \mathbf{b}$ par minimisation de la somme des carrés des erreurs, est la solution du système linéaire suivant où $\mathbf{A}^T\mathbf{A}$ est appelée la matrice de Gram. +\[ +\mathbf{A}^T\mathbf{A} \mathbf{x} = \mathbf{A}^T\mathbf{b} +\] + +# Méthode des moindres carrés appliquée à la régression polynomiale + +Pour un polynôme de degré $N-1$, les fonctions de bases mentionnées ci-dessus sont : $f_1(x)=1$, $f_2(x)=x$, $f_3(x)=x^2$,..., $f_N(x)=x^{N-1}$. Elles permettent de définir la matrice des données $\mathbf{A}$ et la matrice de Gram $\mathbf{A}^T\mathbf{A}$. + +Nous reprenons l'exemple synthétique du précédent chapitre et nous résolvons le système linéaire correspondant à la matrice de Gram pour un polynôme de degré fixé. + +```{r} +set.seed(1123) +# Image par f d'un échantillon uniforme sur l'intervalle [0,1], avec ajout d'un +# bruit gaussien de moyenne nulle et d'écart type 0.2 +data = gendat(10,0.2) +coef = polyreg2(data,3) +plt(data,f) +pltpoly(coef) +``` + +Ce polynôme de degré trois modélise mieux la fonction génératrice inconnue que celui de degré quatre qui ne commettait aucune erreur sur les données observées. \ No newline at end of file diff --git a/03_tikhonov.R b/03_tikhonov.R new file mode 100644 index 0000000..2a47c5d --- /dev/null +++ b/03_tikhonov.R @@ -0,0 +1,10 @@ +# Résolution d'un système linéaire correspondant à la matrice de Gram pour +# un polynôme de degré fixé et avec l'ajout d'un facteur de régularisation en +# norme L2 dont l'importance est contrôlée par l'hyperparamètre alpha. +ridge <- function(alpha, data, degre) { + xs <- c(data$X) + A <- outer(xs, 0:degre, "^") + gram <- t(A) %*% A + diag(gram) <- diag(gram) + alpha + solve(gram, as.vector(t(A) %*% data$Y)) +} \ No newline at end of file diff --git a/03_tikhonov.Rmd b/03_tikhonov.Rmd new file mode 100644 index 0000000..e0fa4ba --- /dev/null +++ b/03_tikhonov.Rmd @@ -0,0 +1,106 @@ +--- +title: "ML 03 Régularisation de Tikhonov" +output: + bookdown::pdf_document2: + number_section: yes +toc: false +classoption: fleqn +--- + +```{r, include=FALSE} +source("ML1_01_intro.R", local = knitr::knit_global()) +source("ML1_03_tikhonov.R", local = knitr::knit_global()) +``` + +# Problèmes linéaires mal posés + +Avec moins d'observations que de fonctions de base ($M1]<-1; # domaine [0,1] de la matrice initiale peut ne plus être respecté + image(t(apply(Y,2,rev)), col=grey(seq(0,1,length=256)), axes=FALSE, ...) +} \ No newline at end of file diff --git a/05_presentation_svd.Rmd b/05_presentation_svd.Rmd new file mode 100644 index 0000000..1b701fc --- /dev/null +++ b/05_presentation_svd.Rmd @@ -0,0 +1,141 @@ +--- +title: "ML 05 Présentation de la décomposition en valeurs singulières (SVD)" +output: + bookdown::pdf_document2: + number_section: yes +toc: false +classoption: fleqn +--- + +```{r, include=FALSE} +source("ML1_05_presentation_svd.R", local = knitr::knit_global()) +``` + +# Contexte + +Nous présentons la décomposition en valeurs singulières, une méthode générique de compression d'un tableau rectangulaire de données. Elle possède de nombreuses applications. En particulier, elle permet de calculer en un seul entraînement les régressions ridge pour toutes les valeurs possibles de l'hyperparamètre $\alpha$. Plus généralement, elle permet souvent d'adopter une perspective géométrique pour mieux comprendre les problèmes et algorithmes de l'apprentissage automatique à partir de données. + +# Recherche d'un sous-espace qui minimise les carrés des écarts à l'espace d'origine + +Soit un tableau de données $\mathbf{X}$ composé de $N$ observations en lignes, chacune décrite par $P$ variables en colonnes. $x_{ij}$ est la valeur de la variable $j$ pour l'observation $i$. Les vecteurs lignes forment $N$ points de l'espace $\mathbb{R}^P$. Les vecteurs colonnes forment $P$ points de l'espace $\mathbb{R}^N$. Nous cherchons à projeter le nuage des points lignes sur un sous-espace $\mathcal{H} \subset \mathbb{R}^P$, tout en minimisant les déformations. + +Pour commencer, nous considérons le meilleur sous-espace $\mathcal{H}$ à une dimension, c'est-à-dire une droite définie par son vecteur directeur unitaire $\mathbf{v}$ ("unitaire" signifie ici que sa norme euclidienne est égale à $1$, soit $\sqrt{\mathbf{v}^T\mathbf{v}}=1$, autrement dit, $\mathbf{v}^T\mathbf{v}=1$). Soit $M_i$ un des $N$ points de $\mathbb{R}^P$. A ce point correspond le vecteur $\mathbf{OM_i}$ aussi noté $\mathbf{x_i}$ car ses coordonnées se lisent sur la i-ème ligne de $\mathbf{X}$. Soit $H_i$ la projection de $M_i$ sur la droite $\mathcal{H}$. + +![svd1](images/svd1.jpeg) + +La longueur $OH_i$ est le produit scalaire de $\mathbf{x_i}$ et de $\mathbf{v}$. +\[ +OH_i = \mathbf{x_i}^T\mathbf{v} = \sum_{j=1}^{P} x_{ij} v_j +\] +Nous obtenons un vecteur des $N$ projections $OH_i$ par le produit de la matrice $\mathbf{X}$ et du vecteur $\mathbf{v}$ : $\mathbf{X}\mathbf{v}$. Choisissons pour $\mathcal{H}$ le sous-espace qui minimise la somme des carrés des écarts entre chaque point et sa projection : $\sum_i M_i H_i^2$. + +Le théorème de Pythagore donne $\sum_i M_i H_i^2 = \sum_i OM_i^2 - \sum_i OH_i^2$. + +Puisque $\sum_i OM_i^2$ est indépendant de $\mathbf{v}$, nous devons maximiser +\begin{equation} +\sum_i OH_i^2 = (\mathbf{X}\mathbf{v})^T(\mathbf{X}\mathbf{v}) = \mathbf{v}^T \mathbf{X}^T \mathbf{X} \mathbf{v} +\label{eq:svd-quadratic} +\end{equation} +...sous la contrainte $\mathbf{v}^T\mathbf{v}=1$. + +Nous notons $\mathbf{v_1}$ ce vecteur directeur du meilleur sous-espace de dimension 1 pour le nuage des $N$ observations de $\mathcal{R}^P$. Le meilleur sous-espace de dimension $2$ contient $\mathbf{v_1}$ et il faut le compléter avec un second vecteur unitaire $\mathbf{v_2}$, orthogonal à $\mathbf{v_1}$ et qui maximise $\mathbf{v_2}^T \mathbf{X}^T \mathbf{X} \mathbf{v_2}$. Etc. + +Nous montrerons plus loin que $\mathbf{v_1}$ est le vecteur propre de $\mathbf{X}^T \mathbf{X}$ associé à la plus grande valeur propre $\lambda_1$, cette dernière étant justement le maximum de la forme quadratique de l'équation (\ref{eq:svd-quadratic}). De même, $\mathbf{v_2}$ est le vecteur propre de $\mathbf{X}^T \mathbf{X}$ associé à la seconde plus grande valeur propre $\lambda_2$. Etc. + +Nous venons de voir que les vecteurs orthogonaux $\mathbf{v_1}, \mathbf{v_2},\dots, \mathbf{v_k}$ forment une base d'un sous-espace $k$-dimensionnel sur lequel peuvent être projetés les $N$ vecteurs lignes de $\mathbf{X}$ pour minimiser les carrés des écarts à l'espace d'origine qui était de dimension (au plus) $max(N,P)$. + +Nous construirons de même des vecteurs orthogonaux $\mathbf{u_1}, \mathbf{u_2},\dots, \mathbf{u_k}$ qui forment une base d'un sous-espace $k$-dimensionnel sur lequel peuvent être projetés les $P$ vecteurs colonnes de $\mathbf{X}$ pour minimiser les carrés des écarts à l'espace d'origine qui était de dimension (au plus) $max(N,P)$. + +Nous trouverons de façon similaire que $\mathbf{u_1}$ est le vecteur propre de $\mathbf{X} \mathbf{X}^T$ associé à la plus grande valeur propre $\mu_1$. De même, $\mathbf{u_2}$ est le vecteur propre de $\mathbf{X} \mathbf{X}^T$ associé à la seconde plus grande valeur propre $\mu_2$. Etc. + +Les vecteurs $\mathbf{u_1}, \mathbf{u_2}, \dots$ et $\mathbf{v_1}, \mathbf{v_2}, \dots$ sont les axes principaux de la matrice des données $\mathbf{X}$. + +Nous pouvons écrire : + +\begin{equation} +\begin{cases} +\mathbf{X}^T\mathbf{X}\mathbf{v_\alpha} = \lambda_\alpha \mathbf{v_\alpha} \\ +\mathbf{X}\mathbf{X}^T\mathbf{u_\alpha} = \mu_\alpha \mathbf{u_\alpha} +\end{cases} +\label{eq:svd-meme-valeurs-propres} +\end{equation} + +En prémultipliant la première des deux équations du système (\ref{eq:svd-meme-valeurs-propres}) par $\mathbf{X}$, nous obtenons : + +\[ +\left(\mathbf{X}\mathbf{X}^T\right)\mathbf{X}\mathbf{v_\alpha} = \lambda_\alpha \left(\mathbf{X}\mathbf{v_\alpha}\right) +\] + +Cette dernière équation nous indique qu'au vecteur propre $\mathbf{v_\alpha}$ de $\mathbf{X}^T\mathbf{X}$ correspond un vecteur propre $\left(\mathbf{X}\mathbf{v_\alpha}\right)$ de $\mathbf{X}\mathbf{X}^T$ de même valeur propre $\lambda_\alpha$. Comme $\mu_1$ est la plus grande valeur propre de $\mathbf{X}\mathbf{X}^T$, nous avons montré que $\lambda_1 \leq \mu_1$ + +De même, en prémultipliant la seconde des deux équations du système (\ref{eq:svd-meme-valeurs-propres}) par $\mathbf{X}^T$, nous montrerions que $\mu_1 \leq \lambda_1$. Donc, $\lambda_1 = \mu_1$ et en général, $\lambda_\alpha = \mu_\alpha$. + +Calculons la norme euclidienne, aussi appelée norme L2, de $\mathbf{X}\mathbf{v_\alpha}$. + +\begin{align*} + & \text{Norme euclidienne de } \mathbf{X}\mathbf{v_\alpha} \\ += \{ & \text{ Par définition de la norme euclidienne.} \} \\ + & \sqrt{\left(\mathbf{X}\mathbf{v_\alpha}\right)^T\left(\mathbf{X}\mathbf{v_\alpha}\right)} \\ += \{ &\text{Propriété de la transposition} \} \\ + & \sqrt{\mathbf{v_\alpha}^T\mathbf{X}^T\mathbf{X}\mathbf{v_\alpha}} \\ += \{ & \text{Voir équation (\ref{eq:svd-meme-valeurs-propres})} \} \\ +& \sqrt{\lambda_\alpha\mathbf{v_\alpha}^T\mathbf{v_\alpha}} \\ += \{ & \text{Par construction, } \mathbf{v_\alpha} \text{est de norme 1.} \} \\ +& \sqrt{\lambda_\alpha} \\ +\end{align*} + +Nous avons montré plus haut que $\left(\mathbf{X}\mathbf{v_\alpha}\right)$ est un vecteur propre de $\mathbf{X}\mathbf{X}^T$ de valeur propre associée $\lambda_\alpha$. Or, $\mathbf{u_\alpha}$ est le vecteur propre unitaire (i.e., de norme euclidienne égale à 1) de $\left(\mathbf{X}\mathbf{v_\alpha}\right)$ associé à la valeur propre $\lambda_\alpha$. C'est pourquoi, nous pouvons écrire : + +\[ +\begin{cases} +\mathbf{u_\alpha} = \frac{1}{\sqrt{\lambda_\alpha}} \mathbf{Xv_\alpha} \\ +\mathbf{v_\alpha} = \frac{1}{\sqrt{\lambda_\alpha}} \mathbf{X}^T\mathbf{u_\alpha} +\end{cases} +\] + +A partir de l'égalité $\mathbf{X}\mathbf{v_\alpha}=\mathbf{u_\alpha}\sqrt{\lambda_\alpha}$, nous post-multiplions par $\mathbf{v_\alpha}^T$ et nous sommons sur l'ensemble des $P$ axes principaux (en supposant, sans perte de généralité, que $N>P$) : + +\[ +\mathbf{X}\left(\sum_{\alpha=1}^{P}\mathbf{v_\alpha}\mathbf{v_\alpha}^T\right) = \sum_{\alpha=1}^{P}\sqrt{\lambda_\alpha}\mathbf{u_\alpha}\mathbf{v_\alpha}^T +\] + +Soit $\mathbf{V}$ la matrice formée des $P$ vecteurs $\mathbf{v_\alpha}$ en colonnes. Comme les vecteurs $\mathbf{v_\alpha}$ sont orthogonaux deux à deux et de norme $1$, $\mathbf{V}\mathbf{V}^T$ est la matrice identité $\mathbf{I_P}$. Donc, $\sum_{\alpha=1}^{P}\mathbf{v_\alpha}\mathbf{v_\alpha}^T=\mathbf{V}\mathbf{V}^T=\mathbf{I_P}$. Nous obtenons finalement : + +\begin{equation} +\mathbf{X} = \sum_{\alpha=1}^{P}\sqrt{\lambda_\alpha}\mathbf{u_\alpha}\mathbf{v_\alpha}^T +\label{eq:svd} +\end{equation} + +Nous avons montré que la matrice $\mathbf{X}$ peut s'écrire comme une somme de matrices de rang 1 qui sont les produits des vecteurs $\mathbf{u}$ (appelés *vecteurs singuliers à gauche*) et $\mathbf{v}$ (appelés *vecteurs singuliers à droite*) pondérés par les *valeurs singulières* $\sqrt{\lambda_\alpha}$. L'équation (\ref{eq:svd}) est celle de la *décomposition en valeurs singulières* de $\mathbf{X}$. + +Nous pouvons écrire ce résultat sous forme matricielle en introduisant les matrices $\mathbf{U}$, $\mathbf{V}$ et $\mathbf{D}$. Sans perte de généralité, nous exprimons les résultats pour $N>P$. La matrice $\mathbf{U}$, de dimension $N \times P$, contient en colonnes les vecteurs singuliers à gauche $\mathbf{u_1}, \mathbf{u_2}, \dots, \mathbf{u_P}$. La matrice $\mathbf{V}$, de dimension $P \times P$ contient en colonnes les vecteurs singuliers à droite $\mathbf{v_1}, \mathbf{v_2}, \dots, \mathbf{v_P}$. La matrice diagonale $\mathbf{D}$ contient sur sa diagonale les valeurs singulières $\sqrt{\lambda_1}, \sqrt{\lambda_2},\dots,\sqrt{\lambda_P}$. La décomposition en valeurs singulière de toute matrice $\mathbf{X}$ s'écrit alors : + +\[ +\mathbf{X} = \mathbf{U} \mathbf{D} \mathbf{V}^T +\] + +Notons que le rang de la matrice $\mathbf{X}$, c'est-à-dire le nombre de colonnes indépendantes (ou le nombre de lignes indépendantes, c'est équivalent), est égal au nombre de valeurs singulières non nulles. + +Par ailleurs, nous obtenons une reconstitution approchée de rang $k$ de $\mathbf{X}$ en ne conservant dans l'équation (\ref{eq:svd}) que les $k$ plus grandes valeurs singulières $\sqrt{\lambda_1},\sqrt{\lambda_2}, \dots, \sqrt{\lambda_k}$ et en annulant toutes les autres. + +## Illustrations du SVD avec la compression d'une image + +Considérons l'image d'un hortensia en fleurs. + +![hortensia](images/hortensia.jpg)\ + +Nous convertissons cette image, originellement au format JPEG, en un format en niveaux de gris, simple à manipuler, appelé PGM (*Portable Grey Map*). Dans ce format non compressé, l'image est représentée par une matrice dont chaque valeur, comprise entre 0 et 1, représente le niveau de gris d'un pixel. Nous opérons cette transformation grâce au programme *imagemagick* avec la commande `convert hortensia.jpg hortensia.pgm`. + +Ensuite, nous appliquons la décomposition en valeurs singulières à la matrice de l'image en niveaux de gris. Nous utilisons le résultat du SVD pour reconstituer une version compressée de l'image en ne conservant qu'un certain nombre des plus grandes valeurs singulières. + +```{r, warning=FALSE} +X <- read_grey_img("images/hortensia.pgm") + +par(mfrow=c(1,2), oma=c(1,1,0,0)+0.1, mar=c(0,0,1,1)+0.1); +print_grey_img(X, main="image originale d=256", asp=1); +print_grey_img(compress_SVD(X,16), main="d=16", asp=1); +print_grey_img(compress_SVD(X,32), main="d=32", asp=1); +print_grey_img(compress_SVD(X,64), main="d=64", asp=1); +print_grey_img(compress_SVD(X,128), main="d=128", asp=1); +print_grey_img(compress_SVD(X,256), main="d=256", asp=1); +``` diff --git a/06_matrice_definie_non_negative.Rmd b/06_matrice_definie_non_negative.Rmd new file mode 100644 index 0000000..82a732e --- /dev/null +++ b/06_matrice_definie_non_negative.Rmd @@ -0,0 +1,71 @@ +--- +title: "ML 06 Matrice définie non négative" +output: + bookdown::pdf_document2: + number_section: yes +toc: false +classoption: fleqn +--- + +Avant de présenter le cœur de la preuve de l'existence de la décomposition en valeurs singulières, nous rappelons quelques propriétés des matrices symétriques définies non négatives. Une matrice définie non négative (souvent noté PSD pour "positive semi-definite") $\mathbf{M}$ représente une métrique car elle définit un produit scalaire et donc une géométrie, c'est-à-dire qu'elle donne un sens aux notions de distance et d'angle. + +Le produit scalaire défini par $\mathbf{M}$ entre deux vecteurs $\mathbf{x}$ et $\mathbf{y}$ pourra se noter $<\mathbf{x},\mathbf{y}>_{\mathbf{M}} = \mathbf{x}^T\mathbf{M}\mathbf{y}$. Pour que ce produit scalaire soit valide, $\mathbf{M}$ doit respecter la contrainte $\mathbf{x}^T\mathbf{M}\mathbf{y} \geq 0$ pour tout $\mathbf{x}$ et $\mathbf{y}$. + +Une fois donné un produit scalaire, les notions liées de norme et d'angle suivent~: $\|\mathbf{x}\|_{\mathbf{M}} = \sqrt{<\mathbf{x},\mathbf{x}>_{\mathbf{M}}}$ et $cos^2(\mathbf{x},\mathbf{y}) = <\mathbf{x},\mathbf{y}>_{\mathbf{M}}/\|\mathbf{x}\|_{\mathbf{M}}\|\mathbf{y}\|_{\mathbf{M}}$. + +La géométrie euclidienne "classique" (associée à la base canonique) correspond à l'emploi de la matrice identité pour métrique, $\mathbf{M}=\mathbf{I}$. + +Aussi, les valeurs propres d'une matrice PSD sont toutes positives~: + +\begin{align*} + & \mathbf{x}^T\mathbf{M}\mathbf{x} \geq 0 \\ +\Rightarrow \{ & \lambda \text{ est une valeur propre de } \mathbf{M} \} \\ + & \mathbf{x}^T\lambda\mathbf{x} \geq 0 \\ += \{ &\text{utilisation de la métrique identité} \} \\ + & \lambda \|x\|^2_{\mathbf{I}} \geq 0 \\ += \{ &\|x\|^2_{\mathbf{I}} \geq 0 \} \\ + & \lambda \geq 0 +\end{align*} + +Soit $\mathbf{V}$ la matrice des vecteurs propres de $\mathbf{M}$ et $\mathbf{L}$ la matrice diagonale de ses valeurs propres positives. + +\begin{align*} + & \mathbf{M}\mathbf{V} = \mathbf{V}\mathbf{L} \\ += \phantom{\{} & \\ + & \mathbf{M} = \mathbf{V}\mathbf{L}\mathbf{V}^{-1} \\ += \{ & \mathbf{S} = \sqrt{\mathbf{L}}\} \\ +& \mathbf{M} = \mathbf{V}\mathbf{S}\mathbf{S}\mathbf{V}^{-1} \\ += \{ & \quad \text{Les vecteurs propres distincts sont orthogonaux deux à deux : } \mathbf{V}^{-1} = \mathbf{V}^T. \\ +\phantom{=}\phantom{ }\phantom{\{} & \quad \mathbf{S}=\mathbf{S}^T. \text{On pose } \mathbf{T}=\mathbf{V}\mathbf{S}. \} \\ +& \mathbf{M} = \mathbf{T}\mathbf{T}^T \\ +\end{align*} + +Ainsi, toute matrice définie non négative $\mathbf{M}$ se décompose en $\mathbf{M}=\mathbf{T}\mathbf{T}^T$ où $\mathbf{T}$ est une transformation linéaire composée d'un changement d'échelle des axes du repère ($\mathbf{S}$) et d'une rotation ($\mathbf{V}$). + +Nous remarquons par exemple que l'équation d'un cercle selon la métrique $\mathbf{M}$ correspond à celle d'une ellipse selon la métrique identité $\mathbf{I}$ : + +\begin{align*} + & \|\mathbf{x}\|_{\mathbf{M}}^2 = c \\ += \{ &\text{C'est l'équation d'un cercle, c est une constante.}\} \\ + & < \mathbf{x}, \mathbf{x}>_{\mathbf{M}} = c \\ += \{ &\text{Par définition du produit scalaire.} \} \\ + & \mathbf{x}^T \mathbf{M} \mathbf{x} = c \\ += \{ &\mathbf{M} = \mathbf{T}^T\mathbf{T}\} \\ + & \mathbf{x}^T \mathbf{T}^T \mathbf{T} \mathbf{x} = c \\ += \{ &\text{Par définition du produit scalaire, on retrouve le cercle déformé en une ellipse par } \mathbf{T} \} \\ + & \|\mathbf{T}\mathbf{x}\|_{\mathbf{I}}^2 = c \\ +\end{align*} + +Ainsi, une matrice symétrique définie non négative projette le cercle unité sur une ellipse dont les axes orthonormaux sont les vecteurs propres et les longueurs des axes sont les valeurs propres. + +Enfin, l'angle formé entre $\mathbf{x}$ et $\mathbf{Mx}$ est inférieur ou égal à $\pi/2$ radians : + +\begin{align*} + & \mathbf{x}^T \mathbf{M} \mathbf{x} \geq 0 \\ += \{&\text{Géométrie du produit scalaire. On note } \theta \text{ l'angle entre les deux vecteurs.}\} \\ + & \|\mathbf{x}\|_\mathbf{I} \; \|\mathbf{Mx}\|_\mathbf{I} \; cos(\theta) \geq 0 \\ +\Rightarrow \phantom{\{}&\\ + & |\theta| \leq \pi/2 +\end{align*} + +Ainsi, le vecteur résultat $\mathbf{Mx}$ reste du même côté que $\mathbf{x}$ par rapport à l'hyperplan perpendiculaire à $\mathbf{x}$ qui sépare l'espace en deux. Nous comprenons que le concept de matrice symétrique définie non négative est une généralisation pour un espace multidimensionnel du concept de nombre positif sur la droite réelle. \ No newline at end of file diff --git a/07_multiplicateurs_de_lagrange.Rmd b/07_multiplicateurs_de_lagrange.Rmd new file mode 100644 index 0000000..4c02d3c --- /dev/null +++ b/07_multiplicateurs_de_lagrange.Rmd @@ -0,0 +1,76 @@ +--- +title: "ML 07 Optimisation sous contraintes et multiplicateurs de Lagrange" +output: + bookdown::pdf_document2: + number_section: yes +toc: false +classoption: fleqn +--- + +Avant de pouvoir présenter la preuve de l'existence de la décomposition en valeurs singulières, nous devons introduire une stratégie souvent très utile pour résoudre un problème d'optimisation sous contraintes : la méthode dite des multiplicateurs de Lagrange. + +Prenons un exemple simple à deux dimensions. Nous cherchons à maximiser $F(\mathbf{x}) = x_1^2 + x_2^2$ sous la contrainte $K(\mathbf{x})=b$, avec $K(\mathbf{x}) = a_1 x_1 + a_2 x_2$. C'est-à-dire que le point solution $\mathbf{x^*}$ doit être situé sur la ligne $K$. + +![lagrange](images/lagrange.jpg) +Sur cet exemple, les lignes de niveaux, ou contours, de $F$ sont des cercles centrés sur l'origine du repère cartésien. +La solution $\mathbf{x*}$ appartient au contour de $F$ tangent à la contrainte $K$. Il faut donc que leurs gradients soient alignés : $\nabla F = \lambda \nabla K$. + +\[ +\nabla F = +\left( \begin{array}{c} +\partial F / \partial x_1 \\ +\partial F / \partial x_2 +\end{array} \right) += +\left( \begin{array}{c} +2 x_1 \\ +2 x_2 +\end{array} \right) +\quad ; \quad +\nabla K = +\left( \begin{array}{c} +\partial K / \partial x_1 \\ +\partial K / \partial x_2 +\end{array} \right) += +\left( \begin{array}{c} +a_1 \\ +a_2 +\end{array} \right) +\] +Nous avons donc un système de trois équations à trois inconnues ($x_1$, $x_2$ et $\lambda$) : + +\[ +\begin{cases} +2 x_1 = \lambda a_1 \\ +2 x_2 = \lambda a_2 \\ +a_1 x_1 + a_2 x_2 = b +\end{cases} +\] + +En résolvant ce système par simples substitutions, nous trouvons la solution : + +\[ +\begin{cases} +x_1^* = \frac{a_1 b}{a_1^2 + a_2^2} \\ +x_2^* = \frac{a_2 b}{a_1^2 + a_2^2} \\ +\lambda^* = \frac{2b}{a_1^2 + a_2^2} +\end{cases} +\] +Ce même calcul s'écrit d'une façon plus systématique en rassemblant les équations $\nabla F = \lambda \nabla K$ et $K(\mathbf{x})=b$ par l'introduction du Lagrangien : + +\[ +\mathcal{L}(\mathbf{x},\lambda) = F(\mathbf{x}) - \lambda(K(\mathbf{x})-b) +\] + +En effet, en annulant les dérivées partielles de $\mathcal{L}$, nous retrouvons les équations $\nabla F = \lambda \nabla K$ et $K(\mathbf{x})=b$ : + +\[ +\begin{array}{lll} +\partial \mathcal{L} / \partial x_1 = 0 &\Leftrightarrow &\partial F / \partial x_1 = \lambda \partial K / \partial x_1\\ +\partial \mathcal{L} / \partial x_2 = 0 &\Leftrightarrow &\partial F / \partial x_2 = \lambda \partial K / \partial x_2\\ +\partial \mathcal{L} / \partial \lambda = 0 &\Leftrightarrow &K(\mathbf{x}) = b\\ +\end{array} +\] + +Ainsi, le problème d'optimisation sous contrainte exprimé en fonction de $F$ et $K$ peut s'écrire comme un problème d'optimisation non contraint, en fonction de $\mathcal{L}$, dans un espace de plus grande dimension que celui d'origine puisque s'ajoute aux dimensions de $\mathbf{x}$ celle de $\lambda$. \ No newline at end of file diff --git a/08_derivation_svd.Rmd b/08_derivation_svd.Rmd new file mode 100644 index 0000000..d716d96 --- /dev/null +++ b/08_derivation_svd.Rmd @@ -0,0 +1,81 @@ +--- +title: "ML 08 Dérivation du SVD" +output: + bookdown::pdf_document2: + number_section: yes +toc: false +classoption: fleqn +--- + +Après ces préambules sur la notion de métrique associée à une matrice définie non négative et sur la méthode d'optimisation dite des multiplicateurs de Lagrange, nous reprenons la dérivation de la preuve de l'existence de la décomposition en valeurs singulières pour toute matrice. + +Rappelons que nous avons noté $\mathbf{v_1}$ le vecteur directeur du meilleur sous-espace de dimension 1 pour le nuage des $N$ observations de $\mathcal{R}^P$. Nous prouvons maintenant que $\mathbf{v_1}$ est le vecteur propre de $\mathbf{X}^T\mathbf{X}$ associé à la plus grande valeur propre $\lambda_1$. + +Le résultat est sans difficulté plus général. Nous le prouvons pour toute matrice symétrique $\mathbf{A}$ (non seulement $\mathbf{X}^T\mathbf{X}$) et toute métrique de $\mathcal{R}^P$ représentée par une matrice symétrique définie non négative $\mathbf{M}$ (non seulement la matrice identité $\mathbf{I}$). + +Nous rappelons que $\mathbf{v}^T\mathbf{A}\mathbf{v} = \Sigma a_{i,j}v_iv_j$. Ainsi, comme $\mathbf{A}$ et $\mathbf{M}$ sont symétriques (i.e., $a_{i,j}=a_{j,i}$ et $m_{i,j}=m_{j,i}$), les dérivées partielles des formes quadratiques ont les formes suivantes : +\[ +\frac{\partial \mathbf{v}^T\mathbf{A}\mathbf{v}}{\partial \mathbf{v}}=2\mathbf{A}\mathbf{v} \quad \text{et} \quad \frac{\partial \mathbf{v}^T\mathbf{M}\mathbf{v}}{\partial \mathbf{v}}=2\mathbf{M}\mathbf{v} +\] +Pour trouver le maximum de $\mathbf{v}^T\mathbf{A}\mathbf{v}$ en intégrant la contrainte unitaire sur $\mathbf{v}$ (i.e., $\mathbf{v}^T\mathbf{M}\mathbf{v}=1$), nous annulons les dérivées du langrangien $\mathcal{L}$ : + +\[ +\mathcal{L} = \mathbf{v}^T\mathbf{A}\mathbf{v} - \lambda (\mathbf{v}^T\mathbf{M}\mathbf{v} - 1) +\] + +\begin{align*} + & \frac{\partial \mathcal{L}}{\partial \mathbf{v}} = 0 \\ += \{& \text{voir calcul des dérivées ci-dessus} \} \\ + & 2\mathbf{A}\mathbf{v} - 2\lambda\mathbf{M}\mathbf{v} = 0 \\ += \phantom{\{} & \\ + & \mathbf{A}\mathbf{v} = \lambda\mathbf{M}\mathbf{v} \\ += \{& \mathbf{v}^T\mathbf{M}\mathbf{v}=1 \} \\ + & \lambda = \mathbf{v}^T\mathbf{A}\mathbf{v} +\end{align*} + +Nous découvrons que la valeur du multiplicateur de Lagrange $\lambda$ est le maximum recherché. Par ailleurs, nous avons : + +\begin{align*} + & \mathbf{A}\mathbf{v} = \lambda\mathbf{M}\mathbf{v} \\ += \{& \mathbf{M} \text{ est définie non négative et donc inversible} \} \\ + & \mathbf{M}^{-1}\mathbf{A}\mathbf{v} = \lambda\mathbf{v} +\end{align*} + +Donc $\mathbf{v}$ est le vecteur propre de $\mathbf{M}^{-1}\mathbf{A}$ associé à la plus grande valeur propre $\lambda$. Nous notons cette solution $\mathbf{v_1}$ et la valeur propre correspondante $\lambda_1$. + +Nous cherchons ensuite $\mathbf{v_2}$, unitaire ($\mathbf{v_2}^T\mathbf{M}\mathbf{v_2}=1$), orthogonal à $\mathbf{v_1}$ ($\mathbf{v_1}^T\mathbf{M}\mathbf{v_2}=0$) et qui maximise $\mathbf{v_2}^T\mathbf{A}\mathbf{v_2}$. Pour ce faire, nous annulons les dérivées du Langrangien ci-après. + +\[ +\mathcal{L} = \mathbf{v_2}^T\mathbf{A}\mathbf{v_2} - \lambda_2 (\mathbf{v_2}^T\mathbf{M}\mathbf{v_2} - 1) - \mu_2\mathbf{v_2}^T\mathbf{M}\mathbf{v_1} +\] + +\begin{align*} + & \frac{\partial \mathcal{L}}{\partial \mathbf{v_2}} = 0 \\ += \{& \text{Par définition de } \mathcal{L} \text{ et car A et M sont symétriques.}\} \\ + & 2\mathbf{A}\mathbf{v_2} - 2\lambda_2\mathbf{M}\mathbf{v_2} - \mu_2\mathbf{M}\mathbf{v_1}= 0 \\ +\Rightarrow \{& \quad \text{Multiplier à gauche chaque membre par } \mathbf{v_1}^T \\ +\phantom{=}\phantom{ }\phantom{\{} & \quad \mathbf{v_1}^T\mathbf{A}\mathbf{v_2}=0 \; ; \; \mathbf{v_1}^T\mathbf{M}\mathbf{v_2}=0 \; ; \; \mathbf{v_1}^T\mathbf{M}\mathbf{v_1}=1 \quad \} \\ + & \mu_2 = 0 +\end{align*} + +Nous avons utilisé ci-dessus la propriété $\mathbf{v_1}^T\mathbf{A}\mathbf{v_2}=0$ que nous prouvons ci-dessous. + +\begin{align*} + & \mathbf{v_1}^T\mathbf{A}\mathbf{v_2} \\ += \{& \mathbf{A} \text{ est symétrique.} \} \\ + & \mathbf{v_2}^T\mathbf{A}\mathbf{v_1} \\ += \{& \mathbf{A}\mathbf{v_1}=\lambda_1\mathbf{M}\mathbf{v_1} \text{ car } \mathbf{v_1} \text{ est vecteur propre de } \mathbf{M}^{-1}\mathbf{A} \} \\ + & \lambda_1 \mathbf{v_2}^T\mathbf{M}\mathbf{v_1} \\ += \{& \mathbf{v_1} \text{ et } \mathbf{v_2} \text{ sont M-orthogonaux.} \} \\ + 0 +\end{align*} + +Ainsi, nous trouvons : + +\begin{align*} + & \mathbf{A}\mathbf{v_2} = \lambda_2 \mathbf{M}\mathbf{v_2} \\ += \{& \mathbf{M} \text{ est définie non négative donc inversible.} \} \\ + & \mathbf{M}^{-1}\mathbf{A}\mathbf{v_2} = \lambda_2 \mathbf{v_2} +\end{align*} + +Donc, $\mathbf{v_2}$ est le vecteur propre de $\mathbf{M}^{-1}\mathbf{A}$ qui correspond à la seconde plus grande valeur propre $\lambda_2$. Le raisonnement est similaire pour $\mathbf{v_3}$, etc. \ No newline at end of file diff --git a/09_puissance_iteree_valeur_propre.Rmd b/09_puissance_iteree_valeur_propre.Rmd new file mode 100644 index 0000000..208c9da --- /dev/null +++ b/09_puissance_iteree_valeur_propre.Rmd @@ -0,0 +1,100 @@ +--- +title: "ML 09 Plus grande valeur propre et puissance itérée" +output: + bookdown::pdf_document2: + number_section: yes +toc: false +classoption: fleqn +--- + +Dans l'intention de comprendre plus tard comment calculer en pratique une décomposition en valeurs singulières, nous commençons par présenter la méthode dite de la puissance itérée pour calculer la plus grande valeur propre et son vecteur propre associé pour une matrice carrée. + +# Décomposition en valeurs propres d'une matrice carrée + +Dire que $\mathbf{u}$ est un vecteur propre de la matrice $\mathbf{A}$ associé à la valeur propre $\lambda$ signifie que $\mathbf{A}\mathbf{u}=\lambda\mathbf{u}$. Remarquons que tout multiple de $\mathbf{u}$ est également un vecteur propre associé à la même valeur propre $\lambda$. Considérons l'exemple de la matrice ci-dessous. +\[ +\mathbf{A}= +\left( \begin{array}{cccc} +1 & 2 & 0 \\ +2 & 1 & 0 \\ +0 & 0 & -1 +\end{array} \right) +\] +Notons $\mathbf{e_1}^T=(1 \; 0 \; 0)$, $\mathbf{e_2}^T=(0 \; 1 \; 0)$ et $\mathbf{e_3}^T=(0 \; 0 \; 1)$. Nous avons : + +\begin{align*} +\mathbf{A}\mathbf{e_3} &= -\mathbf{e_3} \\ +\mathbf{A}(\mathbf{e_1}+\mathbf{e_2}) &= 3(\mathbf{e_1}+\mathbf{e_2}) \\ +\mathbf{A}(\mathbf{e_1}-\mathbf{e_2}) &= -(\mathbf{e_1}-\mathbf{e_2}) \\ +\end{align*} + +Nous observons en particulier que la valeur propre $-1$ est associée à deux vecteurs propres linéairement indépendants : $\mathbf{e_3}$ et $(\mathbf{e_1}-\mathbf{e_2})$. + +# Plus grande valeur propre et méthode de la puissance itérée + +Posons $\mathbf{A}\in\mathbb{R}^{N\times N}$ et $\mathbf{u^{(0)}}\in\mathbb{R}^N$ avec $\|\mathbf{u^{(0)}}\|=1$. Supposons que $\mathbf{u^{(0)}}$ puisse se décomposer en une somme de vecteurs propres de $\mathbf{A}$. + +\begin{align*} +\mathbf{u^{(0)}} &= \mathbf{u_1}+\mathbf{u_2}+\dots +\mathbf{u_r} \\ +\text{Avec } \mathbf{A}\mathbf{u_i} &= \lambda_i \mathbf{u_i}\\ +\end{align*} + +Considérons la séquence $\mathbf{u^{(0)}}, \mathbf{A}\mathbf{u^{(0)}}, \mathbf{A}^2\mathbf{u^{(0)}}, \dots$. Nous avons : + +\begin{equation} +\mathbf{A}^m\mathbf{u^{(0)}} = \lambda_1^m \mathbf{u_1} + \lambda_2^m \mathbf{u_2} + \dots + \lambda_r^m \mathbf{u_r} +\label{eq:1} +\end{equation} + +Supposons par ailleurs que les valeurs propres soient ordonnées et que la première valeur propre soit strictement la plus grande en valeur absolue : $|\lambda_1|>|\lambda_2|\geq\dots\geq|\lambda_r|$. En divisant les deux membres de l'égalité (\ref{eq:1}) par $\lambda_1^m$, nous avons : +\[ +(\lambda_1^{-1}\mathbf{A})^m\mathbf{u^{(0)}}=\mathbf{u_1} + \sum_{i=2}^{r}\left(\frac{\lambda_i}{\lambda_1}\right)^m\mathbf{u_i} +\] +Or, $\left|\frac{\lambda_i}{\lambda_1}\right|<1$ pour $i=2,\dots,r$, donc nous pouvons conclure : +\[ +\lim\limits_{m \to \infty}(\lambda_1^{-1}\mathbf{A})^m\mathbf{u^{(0)}}=\mathbf{u_1} +\] + +En résumé, si $\mathbf{u^{(0)}}$ peut s'écrire comme une somme de vecteurs propres de $\mathbf{A}$ et que la valeur propre $\lambda_1$ associée au vecteur propre $\mathbf{u_1}$ est plus grande en valeur absolue que les autres valeurs propres, alors, à un facteur multiplicatif près, $\mathbf{A}^m\mathbf{u^{(0)}}$ converge vers le vecteur propre $\mathbf{u_1}$. + +Pour gérer ce facteur multiplicatif, nous considérons la séquence suivante : +\[ +\mathbf{w^{(k)}}=\mathbf{A}\mathbf{u^{(k-1)}} \quad ; \quad \mathbf{u^{(k)}}=\frac{\mathbf{w^{(k)}}}{\|\mathbf{w^{(k)}}\|} \quad ; \quad \lambda^{(k)}=\mathbf{u^{(k)}}^T(\mathbf{A}\mathbf{u^{(k)}}) +\] + +Par induction, nous obtenons : + +\begin{align*} + & \mathbf{u^{(k)}} \\ += \{& \gamma_k = \| \lambda_1^k\mathbf{u_1} + \sum_{j=2}^{r} \lambda_j^k\mathbf{u_j} \| \} \\ + & \frac{1}{\gamma_k} \left(\lambda_1^k\mathbf{u_1} + \sum_{j=2}^{r} \lambda_j^k\mathbf{u_j}\right) \\ += \phantom{\{}& \\ + & \frac{\lambda_1^k}{\gamma_k} \left(\mathbf{u_1} + \sum_{j=2}^{r} \left(\frac{\lambda_j}{\lambda_1}\right)^k \mathbf{u_j}\right) \\ +\end{align*} + +Par ailleurs, $\mathbf{u^{(k)}}$ étant à chaque itération normalisé, nous avons : + +\begin{align*} + & \|\mathbf{u^{(k)}}\| = 1 \\ += \phantom{\{}& \\ + & \frac{|\lambda_1^k|}{\gamma_k} = \frac{1}{\|\mathbf{u_1} + \sum_{j=2}^{r} \left(\frac{\lambda_j}{\lambda_1}\right)^k \mathbf{u_j}\|} +\end{align*} + +D'où : +\[ +\lim\limits_{k \to \infty} \frac{|\lambda_1^k|}{\gamma_k} = \frac{1}{\|\mathbf{u_1}\|} +\] + +Finalement : + +\[ +\lim\limits_{k \to \infty} \mathbf{u^{(k)}} = \frac{\mathbf{u_1}}{\|\mathbf{u_1}\|} +\] + +Et : + +\begin{align*} + & \lim\limits_{k \to \infty} \mathbf{u^{(k)}}^T(\mathbf{A}\mathbf{u^{(k)}}) \\ += \{& \mathbf{A}\mathbf{u_1}=\lambda_1\mathbf{u_1} \} \\ + & \lambda_1 \\ +\end{align*} \ No newline at end of file diff --git a/10_puissance_iteree_svd.Rmd b/10_puissance_iteree_svd.Rmd new file mode 100644 index 0000000..db514ff --- /dev/null +++ b/10_puissance_iteree_svd.Rmd @@ -0,0 +1,88 @@ +--- +title: "ML 10 Méthode de la puissance itérée et SVD" +output: + bookdown::pdf_document2: + number_section: yes +toc: false +classoption: fleqn +--- + +Nous montrons comment adapter la méthode de la puissance itérée au calcul de la plus grande valeur singulière et de ses deux vecteurs singuliers associés. Nous considérons une matrice de données rectangulaire $\mathbf{A}\in\mathbb{R}^{N\times M}$. Nous notons $\mathbf{U}\in\mathbb{R}^{N\times N}$ et $\mathbf{V}\in\mathbb{R}^{M\times M}$ les deux matrices orthogonales qui contiennent les vecteurs singuliers en colonnes. Nous notons enfin $\mathbf{\Sigma} = diag(\sigma_1,\sigma_2,\dots,\sigma_r,0,\dots)\in\mathbb{R}^{N\times M}$ la matrice qui contient sur sa diagonale principale les valeurs singulières, avec $r$ le range de $\mathbf{A}$. Nous avons donc : + +\begin{align*} +\mathbf{A} &= \mathbf{U}\mathbf{\Sigma}\mathbf{V}^T & \text{$r$ est le rang de $\mathbf{A}$.} \\ +\mathbf{U} &= [\mathbf{u_1}, \mathbf{u_2},\dots,\mathbf{u_N}] & \text{avec $\mathbf{u_i}\in\mathbb{R}^N$, base orthonormale de $\mathbb{R}^N$} \\ +\mathbf{V} &= [\mathbf{v_1}, \mathbf{v_2},\dots,\mathbf{v_M}] & \text{avec $\mathbf{v_j}\in\mathbb{R}^M$, base orthonormale de $\mathbb{R}^M$} \\ +\mathbf{A} &= \sum_{k=1}^{r} \sigma_k \mathbf{u_k}\mathbf{v_k}^T & +\end{align*} + +Nous allons montrer que la séquence ci-dessous converge vers $\sigma_1$, $\mathbf{u}$ et $\mathbf{v}$ tels que $\mathbf{A}\mathbf{v}=\sigma_1\mathbf{u}$. + +\begin{align*} +\mathbf{w^{(k)}} &= \mathbf{A}\mathbf{v^{(k-1)}} & \alpha_k &= \|\mathbf{w^{(k)}}\|^{-1} & \mathbf{u^{(k)}} &= \alpha_k \mathbf{w^{(k)}} \\ +\mathbf{z^{(k)}} &= \mathbf{A}^T\mathbf{u^{(k)}} & \beta_k &= \|\mathbf{z^{(k)}}\|^{-1} & \mathbf{v^{(k)}} &= \beta_k \mathbf{z^{(k)}} \\ +err &= \|\mathbf{A}\mathbf{v^{(k)}} - \beta_k \mathbf{u^{(k)}}\| & \sigma_1 &= \beta_k & & +\end{align*} + +Nous décomposons $\mathbf{v^{(0)}}\in\mathbb{R}^M$ dans la base orthonormée formée par les vecteurs $\mathbf{v_j}$ : +\[ +\mathbf{v^{(0)}} = \sum_{j=1}^{M} y_j\mathbf{v_j} \quad \text{avec $y_j = \mathbf{v_j}^T\mathbf{v^{(0)}}$} +\] + +Partant de ce $\mathbf{v^{(0)}}$, nous calculons les deux premières itérations : + +\begin{align*} +\mathbf{w^{(1)}} &= \mathbf{A}\mathbf{v^{(0)}} = \sum_{j=1}^{r}\sigma_j y_j \mathbf{u_j} & \mathbf{u^{(1)}} &= \alpha_1 \sum_{j=1}^{r}\sigma_j y_j \mathbf{u_j} \\ +\mathbf{z^{(1)}} &= \mathbf{A}^T\mathbf{u^{(1)}} = \alpha_1 \sum_{j=1}^{r}\sigma_j^2 y_j \mathbf{v_j} & \mathbf{v^{(1)}} &= \alpha_1 \beta_1 \sum_{j=1}^{r}\sigma_j^2 y_j \mathbf{v_j} \\ +\mathbf{w^{(2)}} &= \mathbf{A}\mathbf{v^{(1)}} = \alpha_1 \beta_1 \sum_{j=1}^{r}\sigma_j^3 y_j \mathbf{u_j} & \mathbf{u^{(2)}} &= \alpha_2 \alpha_1 \beta_1 \sum_{j=1}^{r}\sigma_j^3 y_j \mathbf{u_j} \\ +\mathbf{z^{(2)}} &= \mathbf{A}^T\mathbf{u^{(2)}} = \alpha_2 \alpha_1 \beta_1 \sum_{j=1}^{r}\sigma_j^4 y_j \mathbf{v_j} & \mathbf{v^{(2)}} &= \alpha_2 \beta_2 \alpha_1 \beta_1 \sum_{j=1}^{r}\sigma_j^4 y_j \mathbf{v_j} \\ +\end{align*} + +Par induction, nous avons : + +\begin{align*} +\mathbf{v^{(k)}} &= \delta_{2k} \sum_{j=1}^{r} \sigma_j^{2k} y_j \mathbf{v_j} & \text{Avec } \delta_{2k} &= \alpha_k \beta_k \alpha_{k-1} \beta_{k-1} \dots \alpha_1 \beta_1 \\ +\mathbf{u^{(k)}} &= \delta_{2k-1} \sum_{j=1}^{r} \sigma_j^{2k-1} y_j \mathbf{u_j} & \text{Avec } \delta_{2k-1} &= \alpha_{k+1} \alpha_k \beta_k \alpha_{k-1} \beta_{k-1} \dots \alpha_1 \beta_1 \\ +\end{align*} + +Par ailleurs, nous avons : + +\begin{align*} + & 1 \\ += \{& \text{A chaque itération, la norme de $\mathbf{u^{(k)}}$ est maintenue égale à 1.} \} \\ + & \|\mathbf{u^{(k)}}\|^2 \\ += \{& \text{Par définition de $\mathbf{u^{(k)}}$ et car } \|\mathbf{u_j}\| = 1 \} \\ + & \delta_{2k-1}^2 \sum_{j=1}^{r} \sigma_j^{4k-2} y_j^2 \\ +\end{align*} + +De même : + +\begin{align*} + & 1 \\ += \phantom{\{}& \\ + & \|\mathbf{v^{(k)}}\|^2 \\ += \phantom{\{}& \\ + & \delta_{2k}^2 \sum_{j=1}^{r} \sigma_j^{4k} y_j^2 \\ +\end{align*} + +Donc : +\[ +\frac{\|\mathbf{u^{(k+1)}}\|^2}{\|\mathbf{v^{(k)}}\|^2} = 1 = +\sigma_1^2 \left(\frac{\delta_{2k+1}^2}{\delta_{2k}^2}\right) +\left(\frac{\sum_{j=1}^{n_1} y_j^2 \; + \; \sum_{j=n_1}^{r}\left(\frac{\sigma_j}{\sigma_1}\right)^{4k+2}y_j^2} +{\sum_{j=1}^{n_1} y_j^2 \; + \; \sum_{j=n_1}^{r}\left(\frac{\sigma_j}{\sigma_1}\right)^{4k}y_j^2}\right) +\] + +Ci-dessus, $n_1$ est la multiplicité de $\sigma_1$, c'est-à-dire le nombre de vecteurs singuliers indépendants associés à cette valeur singulière. En pratique, sur de "vraies" données, cette multiplicité sera presque toujours égale à $1$. + +Nous avons ainsi montré que : +\[ +\lim\limits_{k \to \infty} \frac{\delta_{2k+1}}{\delta_{2k}} = \sigma_1 +\] + +Et : +\begin{align*} + & \lim\limits_{k \to \infty} \|\mathbf{A}\mathbf{v^{(k)}} - \sigma_1 \mathbf{u^{(k)}}\| \\ += \{& \mathbf{A}\mathbf{v^{(k)}} = \frac{\delta_{2k+1}}{\delta_{2k}} \mathbf{u^{(k+1)}}\} \\ + & 0 +\end{align*} \ No newline at end of file diff --git a/11_projecteurs.Rmd b/11_projecteurs.Rmd new file mode 100644 index 0000000..b7fa80c --- /dev/null +++ b/11_projecteurs.Rmd @@ -0,0 +1,138 @@ +--- +title: "ML 11 Projecteurs orthogonaux et SVD" +output: + bookdown::pdf_document2: + number_section: yes +toc: false +classoption: fleqn +--- + +Nous introduisons la notion de projecteur et nous montrons son lien avec la décomposition en valeurs signulières (SVD). Ainsi, nous comprendrons mieux la géométrie du SVD et nous serons plus tard en mesure de construire un algorithme qui est une extension de la méthode de la puissance itérée pour le calcul des $k$ premières valeurs singulières et vecteurs singuliers. + +# Projecteur + +Un projecteur est une matrice carrée idempotente : $\mathbf{P}^2=\mathbf{P}$. Rappelons que l'espace des colonnes de $\mathbf{P}$ formé de l'ensemble des vecteurs qui s'expriment comme une combinaison linéaire des colonnes de $\mathbf{P}$, c'est-à-dire tous les vecteurs $\mathbf{v}$ qui s'écrivent sous la forme $\mathbf{P}\mathbf{x}$. L'espace des colonnes de $\mathbf{P}$ est également appelé l'image de $\mathbf{P}$ et se note $Im(\mathbf{P})$. Un vecteur $\mathbf{v}$ qui appartient à l'image d'un projecteur $\mathbf{P}$ est tel que $\mathbf{P}\mathbf{v}=\mathbf{P}$ : + +\begin{align*} + & \mathbf{v} \text{appartient à l'image de $\mathbf{P}$} \\ += \{& \text{Par définition de l'image d'une transformation linéaire.} \} \\ + & \mathbf{v} = \mathbf{P}\mathbf{x} \\ += \{& \text{Multiplication des deux membres de l'égalité par $\mathbf{P}$} \} \\ + & \mathbf{P}\mathbf{v} = \mathbf{P^2}\mathbf{x} \\ += \{& \text{Idempotence d'un projecteur.} \} \\ + & \mathbf{P}\mathbf{v} = \mathbf{P}\mathbf{x} \\ += \{& \mathbf{v} = \mathbf{P}\mathbf{x} \} \\ + & \mathbf{P}\mathbf{v} = \mathbf{v} \\ +\end{align*} + +On appelle noyau d'une transformation linéaire, l'ensemble des vecteurs que cette dernière transforme en $\mathbf{0}$, c'est-à-dire les vecteurs qui montrent une dépendance linéaire entre des colonnes de la représentation matricielle de la transformation. Dans le cas d'un projecteur $\mathbf{P}$, la projection sur l'image de $\mathbf{P}$ d'un vecteur $\mathbf{v}$ se fait selon la direction $\mathbf{P}\mathbf{v}-\mathbf{v}$ qui appartient au noyau de $\mathbf{P}$ noté $Ker(\mathbf{P})$ (voir la figure \@ref(fig:projecteur)) : +\[ +\mathbf{P}(\mathbf{P}\mathbf{v}-\mathbf{v}) = \mathbf{P}^2\mathbf{v}-\mathbf{P}\mathbf{v} = \mathbf{P}\mathbf{v}-\mathbf{P}\mathbf{v} = \mathbf{0} +\] + +```{r projecteur, fig.cap= "Schéma d'un projecteur", echo=FALSE} +knitr::include_graphics("images/projecteur1.jpg") +``` + +Le projecteur complémentaire à $\mathbf{P}$ est $\mathbf{I-P}$, il projette sur le noyau de $\mathbf{P}$. C'est-à-dire que $Im(\mathbf{I-P})=Ker(\mathbf{P})$. En effet : + +\begin{align*} + & \mathbf{v} \in Ker(\mathbf{P}) \\ += \{& \text{Par définition du noyau} \} \\ + & \mathbf{P}\mathbf{v} = \mathbf{0} \\ +\Rightarrow \phantom{\{}& \\ + & (\mathbf{I-P})\mathbf{v} = \mathbf{v} \\ += \{& \text{Par définition de l'image} \} \\ + & \mathbf{v} \in Im(\mathbf{I-P}) \\ +\end{align*} + +Donc, $Ker(\mathbf{P}) \subseteq Im(\mathbf{I-P})$. Par ailleurs, un élément de l'image de $\mathbf{I-P}$ peut se noter $(\mathbf{I-P})\mathbf{v}$, ce qui est égal à $\mathbf{v}-\mathbf{Pv}$ et appartient donc au noyau de $\mathbf{P}$. Donc, nous avons aussi $Im(\mathbf{I-P}) \subseteq Ker(\mathbf{P})$. Ainsi : $Ker(\mathbf{P}) = Im(\mathbf{I-P})$ (de même : $Ker(\mathbf{I-P}) = Im(\mathbf{P})$). + +Donc, un projecteur $\mathbf{P}$ sépare l'espace initial en deux sous-espaces $S1 = Im(\mathbf{P})$ et $S2 = Ker(\mathbf{P})$. $\mathbf{P}$ projette **sur** $S1$ **parallèlement à** (ou **le long de**) $S2$. + +# Projecteur orthogonal + +Un projecteur est dit orthogonal quand $S1$ et $S2$ sont orthogonaux. Montrons qu'une projection $\mathbf{P}$ est orthogonale si et seulement si $\mathbf{P}=\mathbf{P}^T$. Montrons d'abord que $\mathbf{P}=\mathbf{P}^T \Rightarrow \text{orthogonalité}$. + +\begin{align*} + & (\mathbf{Px})^T(\mathbf{I-P})\mathbf{y} \quad \text{pour tout $\mathbf{x}$ et $\mathbf{y}$} \\ += \{& \text{Propriété de la transposition.} \} \\ + & \mathbf{x}^T\mathbf{P}^T(\mathbf{I-P})\mathbf{y} \\ += \{& \mathbf{P}=\mathbf{P}^T \} \\ + & \mathbf{x}(\mathbf{P}-\mathbf{P}^2)\mathbf{y} \\ += \{& \mathbf{P} \text{ est un projecteur donc } \mathbf{P}=\mathbf{P}^2\} \\ + & 0 +\end{align*} + +Montrons ensuite que $\text{orthogonalité} \Rightarrow \mathbf{P}=\mathbf{P}^T$. $\mathbf{P}$ projette sur $S1$ le long de $S2$ avec $S1 \bot S2$. Soit $\{\mathbf{q_1}, \mathbf{q_2}, \dots, \mathbf{q_n}\}$ une base orthonormale de $S1$ et $\{\mathbf{q_{n+1}}, \mathbf{q_{n+2}}, \dots, \mathbf{q_m}\}$ une base de $S2$. Soit $Q$ la matrice orthogonale dont la jème colonne est $\mathbf{q_j}$. Nous avons : $\mathbf{PQ} = (\mathbf{q_1}, \mathbf{q_2}, \dots, \mathbf{q_n}, 0, \dots, 0)$. Donc, $\mathbf{Q^TPQ} = diag(1, 1, \dots, 1, 0, \dots, 0)=\mathbf{\Sigma}$. Ainsi, $\mathbf{P} = \mathbf{Q\Sigma Q^T} = \mathbf{P}^T$. + +Etant donné un vecteur unitaire $\mathbf{q}$. Le projecteur orthogonal sur l'espace à une dimension défini par $\mathbf{q}$ est la matrice de rang 1 $\mathbf{P_q}=\mathbf{q}\mathbf{q}^T$. Le projecteur complémentaire est $\mathbf{P_{\bot q}}=\mathbf{I}-\mathbf{P_q}$ de rang $m-1$ (si $\mathbf{q}$ est de dimension $m$). + +Pour toute matrice $\mathbf{Q} \in \mathbb{R}^{M \times N}$ dont les colonnes $\mathbf{q_j}$ sont orthonormales, $\mathbf{P}=\mathbf{QQ^T}=\sum \mathbf{q_j}\mathbf{q_j}^T$ est un projecteur orthogonal sur $Im(\mathbf{Q})$. + +Rappelons la décomposition en valeurs singulière d'une matrice $\mathbf{A} \in \mathbb{R}^{M \times N}$ de rang $r$ : $\mathbf{A} = \mathbf{U} \mathbf{\Sigma}\mathbf{V}^T$ avec $\mathbf{U} \in \mathbb{R}^{M \times M}$, $\mathbf{V} \in \mathbb{R}^{N \times N}$ des matrices aux colonnes orthonormales et $\mathbf{\Sigma}=diag(\sigma_1,\dots,\sigma_r) \in \mathbb{R}^{M \times N}$. Nous pouvons partitionner le SVD en blocs pour obtenir la forme réduite présentée sur la figure \@ref(fig:svdreduit). + +```{r svdreduit, fig.cap= "SVD réduit", echo=FALSE} +knitr::include_graphics("images/SVD_reduit.jpg") +``` + +Nous avons les propriétés fondamentales suivantes : + +\begin{align*} +Im(\mathbf{A}) &= Im(\mathbf{U_r}) & \mathbf{U_r}\mathbf{U_r}^T &: \text{ est un projecteur orthogonal sur } Im(\mathbf{A}) \\ +Ker(\mathbf{A}) &= Im(\tilde{\mathbf{V_r}}) & \tilde{\mathbf{V_r}}\tilde{\mathbf{V_r}}^T &: \text{ est un projecteur orthogonal sur } Ker(\mathbf{A}) \\ +Im(\mathbf{A}^T) &= Im(\mathbf{V_r}) & \mathbf{V_r}\mathbf{V_r}^T &: \text{ est un projecteur orthogonal sur } Im(\mathbf{A}^T) \\ +Ker(\mathbf{A}^T) &= Im(\tilde{\mathbf{U_r}}) & \tilde{\mathbf{U_r}}\tilde{\mathbf{U_r}}^T &: \text{ est un projecteur orthogonal sur } Ker(\mathbf{A}^T) \\ +\end{align*} + +Montrons par exemple que $Ker(\mathbf{A})=Im(\tilde{\mathbf{V_r}})$, c'est-à-dire que les $N-r$ dernières colonnes de $\tilde{\mathbf{V_r}}$ forment une base orthonormale pour le noyau de $\mathbf{A}$, c'est-à-dire pour l'ensemble des vecteurs $\mathbf{x}$ tels que $\mathbf{A}\mathbf{x}=0$. + +Considérons $\mathbf{x} \in Im(\tilde{\mathbf{V_r}})$, ce vecteur s'écrit sous la forme $\mathbf{x} = \sum_{i=r+1}^{N}w_i\mathbf{v_i}$ avec $w_i \in \mathbb{R}$ et $\mathbf{v_i}$ la i-ème colonne de $\tilde{\mathbf{V_r}}$. Nous devons montrer que $\mathbf{A}\mathbf{x}=0$. + +\begin{align*} + & \mathbf{A}\mathbf{x} \\ += \{& \text{Par définition de $\mathbf{x}$.} \} \\ + & \mathbf{A}\sum_{i=r+1}^{N}w_i\mathbf{v_i} \\ += \{& \text{Par linéarité.} \} \\ + & \sum_{i=r+1}^{N}w_i(\mathbf{A}\mathbf{v_i}) \\ +\end{align*} + +Il suffit donc de montrer que chaque $\mathbf{A}\mathbf{v_i}$ est nul. + +\begin{align*} + & \mathbf{A}\mathbf{v_i} \\ += \{& \text{SVD de $\mathbf{A}$} \} \\ + & (\mathbf{U}\mathbf{\Sigma}\mathbf{V}^T) \mathbf{v_i} \\ += \{& \text{Les colonnes de $\mathbf{V}$ sont orthonormales et $\mathbf{v_i}$ est une colonne de $\mathbf{V}$.} \\ +\phantom{=}\phantom{\{}& \text{Donc $\mathbf{V}^T \mathbf{v_i} = \mathbf{e_i}$ le vecteur avec des 0 partout sauf le 1 en position $i$.} \} \\ + & \mathbf{U}\mathbf{\Sigma}\mathbf{e_i} \\ += \{& \mathbf{\Sigma}\mathbf{e_i} = 0 \text{ car à partir de la (r+1)ème ligne, toutes les lignes de $\mathbf{\Sigma}$ sont nulles et $i \geq r+1$.} \} \\ + & 0 \\ +\end{align*} + +# Projection sur une base quelconque + +Nous cherchons à projeter un vecteur $\mathbf{v}$ sur l'espace des colonnes d'une matrice $\mathbf{A}$. Notons le résultat de cette projection $\mathbf{y} \in Im(\mathbf{A})$. $\mathbf{y}-\mathbf{v}$ est orthogonal à $Im(\mathbf{A})$ (c'est-à-dire qu'il appartient à $Ker(\mathbf{A})$). Autrement dit, en posant $\mathbf{y}=\mathbf{A}\mathbf{x}$, nous avons : + +\begin{align*} + & \mathbf{a_j}^T (\mathbf{y}-\mathbf{v}) = 0, \forall j \\ += \{& \text{Posons $\mathbf{y}=\mathbf{A}\mathbf{x}$.} \} \\ + & \mathbf{a_j}^T (\mathbf{A}\mathbf{x}-\mathbf{v}) = 0, \forall j \\ += \{& \text{Linéarité de la multiplication matricielle.} \} \\ + & \mathbf{A}^T (\mathbf{A}\mathbf{x}-\mathbf{v}) = \mathbf{0} \\ += \phantom{\{}& \\ + & \mathbf{A}^T\mathbf{A}\mathbf{x} = \mathbf{A}^T\mathbf{v} \\ += \{& \text{Si $\mathbf{A}^T\mathbf{A}$ est inversible.} \} \\ + & \mathbf{x} = (\mathbf{A}^T\mathbf{A})^{-1}\mathbf{A}^T\mathbf{v} \\ +\end{align*} + +Ainsi, la projection $\mathbf{y} = \mathbf{A}\mathbf{x} = \mathbf{A}(\mathbf{A}^T\mathbf{A})^{-1}\mathbf{A}^T\mathbf{v}$. Nous trouvons donc que $\mathbf{A}(\mathbf{A}^T\mathbf{A})^{-1}\mathbf{A}^T$ est un projecteur orthogonal sur l'espace des colonnes de $\mathbf{A}$. En fait, il est simple de vérifier que ce projecteur est $\mathbf{U_r}\mathbf{U_r}^T$ : + +\begin{align*} +\mathbf{A} &= \mathbf{U_r}\mathbf{\Sigma_r}\mathbf{V_r}^T \\ +\mathbf{A}^T &= \mathbf{V_r}\mathbf{\Sigma_r}\mathbf{U_r}^T \\ +\mathbf{A}^T\mathbf{A} &= \mathbf{V_r}\mathbf{\Sigma_r}^2\mathbf{V_r}^T \\ +(\mathbf{A}^T\mathbf{A})^{-1} &= \mathbf{V_r}\mathbf{\Sigma_r}^{-2}\mathbf{V_r}^T \\ +(\mathbf{A}^T\mathbf{A})^{-1}\mathbf{A}^T &= \mathbf{V_r}\mathbf{\Sigma_r}^{-1}\mathbf{U_r}^T \\ +\mathbf{A}(\mathbf{A}^T\mathbf{A})^{-1}\mathbf{A}^T &= \mathbf{U_r}\mathbf{U_r}^T \\ +\end{align*} \ No newline at end of file diff --git a/12_factorisation_qr.R b/12_factorisation_qr.R new file mode 100644 index 0000000..20945d4 --- /dev/null +++ b/12_factorisation_qr.R @@ -0,0 +1,28 @@ +# Norm 2 +Norm <- function(x) { + stopifnot(is.numeric(x)) + sqrt(sum(x^2)) +} + +# Modified Gram-Schmidt +mgs <- function(A, tol = .Machine$double.eps^0.5) { + stopifnot(is.numeric(A), is.matrix(A)) + M <- nrow(A); N <- ncol(A) + if(M &: \text{l'espace (la ligne) couvert par la première colonne de $\mathbf{A}$} \\ +<\mathbf{a_1}, \mathbf{a_2}> &: \text{l'espace (le plan) couvert par les deux premières colonne de $\mathbf{A}$} \\ +\text{Etc.} &: \text{Etc.} \\ +\end{align*} + +Nous avons $<\mathbf{a_1}> \subseteq <\mathbf{a_1}, \mathbf{a_2}> \subseteq \dots \subseteq <\mathbf{a_1}, \mathbf{a_2}, \dots, \mathbf{a_N}>$. + +Nous cherchons des vecteurs orthogonaux $\mathbf{q_i}$ tels que : +\[ +<\mathbf{q_1}, \mathbf{q_2}, \dots, \mathbf{q_j}> = <\mathbf{a_1}, \mathbf{a_2}, \dots, \mathbf{a_j}> \quad \text{pour } j=1,2,\dots,N +\] + +Ce qui peut s'écrire sous une forme matricielle : +\[ +(\mathbf{a_1} \; \mathbf{a_2} \; \dots \; \mathbf{a_N}) = +(\mathbf{q_1} \; \mathbf{q_2} \; \dots \; \mathbf{q_N}) +\left( \begin{array}{cccc} +r_{11} & r_{12} & \dots & r_{1N} \\ +0 & r_{22} & \dots & r_{2N} \\ +\dots & \dots & \dots & \dots \\ +0 & 0 & \dots & r_{NN} +\end{array} \right) +\] + +Sans perte de généralité, considérons la situation où $M>N$. Nous appelons alors cette expression, $\mathbf{A}=\mathbf{\hat{Q}}\mathbf{\hat{R}}$ (avec $\mathbf{\hat{Q}} \in \mathbb{R}^{M \times N}$ et $\mathbf{\hat{R}} \in \mathbb{R}^{N \times N}$), la forme réduite de la factorisation QR. + +Nous pouvons aussi considérer la forme dite complète de la factorisation QR, avec $\mathbf{Q} \in \mathbb{R}^{M \times M}$ et $\mathbf{R} \in \mathbb{R}^{M \times N}$. Dans ce cas, les $M-N$ dernières lignes de $\mathbf{R}$ sont nulles et il suffit de compléter $\mathbf{\hat{Q}}$ en lui ajoutant $M-N$ colonnes orthogonales. + +# Orthogonalisation de Gram-Schmidt + +Nous présentons l'algorithme de Gram-Schmidt pour la décomposition QR. + +Le principe est de construire $\mathbf{q_j}$ orthogonal à $\mathbf{q_1}, \mathbf{q_2}, \dots, \mathbf{q_{j-1}}$ en retirant à la colonne $\mathbf{a_j}$ ses composantes selon les directions de $\mathbf{q_1}, \mathbf{q_2}, \dots, \mathbf{q_{j-1}}$ : + +\begin{align*} +\mathbf{v_j} &= \mathbf{a_j} - (\mathbf{q_1}^T\mathbf{a_j})\mathbf{q_1} - (\mathbf{q_2}^T\mathbf{a_j})\mathbf{q_2} - \dots - (\mathbf{q_{j-1}}^T\mathbf{a_j})\mathbf{q_{j-1}} \\ +\mathbf{q_j} &= \mathbf{v_j} / \|\mathbf{v_j}\| \\ +\end{align*} + +Nous obtenons ainsi la matrice $\mathbf{\hat{Q}}$. Pour $\mathbf{\hat{R}}$, nous avons, par construction des $\mathbf{q_j}$ : + +\begin{align*} +r_{ij} &= \mathbf{q_i}^T \mathbf{a_j} &\phantom{,} i \neq j \\ +|r_{jj}| &= \| \mathbf{a_j} - \sum_{i=1}^{j-1} r_{ij} \mathbf{q_i} \| &\phantom{,} \text{par convention, nous choisissons $r_{jj}>0$} \\ +\end{align*} + +Nous avons ainsi dérivé l'algorithme classique de Gram-Schmidt : + +\begin{algorithm}[H] +\DontPrintSemicolon +\NoCaptionOfAlgo +\SetAlgoLined +\SetKwInput{Res}{Résultat} +\SetKwInOut{Input}{Entrée}\SetKwInOut{Output}{Sortie} +\SetKwFor{Pour}{pour}{faire}{finpour} +\SetKw{KwTo}{à} +\Res{Décomposition QR de la matrice A} +\Input{$\mathbf{A} \in \mathbb{R}^{M \times N}$, avec $M \geq N$} +\Output{$\mathbf{Q} \in \mathbb{R}^{M \times N}$ une matrice orthonormale \\ +$\mathbf{R} \in \mathbb{R}^{N \times N}$ une matrice triangulaire supérieure \\ +telles que $\mathbf{A} = \mathbf{Q} \mathbf{R}$} +\BlankLine +\Pour{$j \gets 1$ \KwTo $N$}{ + $\mathbf{v_j} \gets \mathbf{a_j}$ \; + \Pour{$i \gets 1$ \KwTo $j-1$}{ + $r_{ij} \gets \mathbf{q_i}^T\mathbf{a_j}$ \; + $\mathbf{v_j} \gets \mathbf{v_j} - r_{ij}\mathbf{q_i}$ \; + } + $r_{jj} \gets \| \mathbf{v_j} \|$ \; + $\mathbf{q_j} \gets \mathbf{v_j} / r_{jj}$ \; +} +\caption{Gram-Schmidt classique (GS)} +\end{algorithm} + +# Algorithme de Gram-Schmidt modifié + +## Dérivation de l'algorithme de Gram-Schmidt modifié + +Nous présentons une autre manière de voir ce résultat qui nous mènera à un algorithme équivalent mais numériquement plus stable. + +Nous pouvons considérer que $\mathbf{q_j}$ est obtenu en projetant $\mathbf{a_j}$ sur l'espace orthogonal à $<\mathbf{a_1}, \mathbf{a_2}, \dots, \mathbf{a_{j-1}}>$, cet espace étant par construction égal à $<\mathbf{q_1}, \mathbf{q_2}, \dots, \mathbf{q_{j-1}}>$ : +\[ +\mathbf{q_j} = \frac{\mathbf{P_j}\mathbf{a_j}}{\|\mathbf{P_j}\mathbf{a_j}\|} \quad ; \quad +\mathbf{P_j} = \mathbf{I} - \mathbf{Q_{j-1}}\mathbf{Q_{j-1}}^T \; \text{avec } +\mathbf{Q_{j-1}} = (\mathbf{q_1} \; \mathbf{q_2} \; \dots \mathbf{q_{j-1}}) +\] + +De façon équivalente, $\mathbf{P_j}$ peut s'exprimer comme le produit de projecteurs sur les espaces orthogonaux à $\mathbf{q_1}$, puis $\mathbf{q_2}$, etc. jusqu'à $\mathbf{q_{j-1}}$ : +\[ +\mathbf{P_j} = \mathbf{P_{\bot \mathbf{q_{j-1}}}} \mathbf{P_{\bot \mathbf{q_{j-2}}}} \dots + \mathbf{P_{\bot \mathbf{q_{2}}}} \mathbf{P_{\bot \mathbf{q_{1}}}} \quad + \text{avec } \mathbf{P_{\bot \mathbf{q}}} = \mathbf{I} - \mathbf{q}\mathbf{q}^T +\] + +Dans la version classique de Gram-Schmidt, $\mathbf{v_j}=\mathbf{P_j}\mathbf{a_j}$. + +Dans la version modifiée de Gram-Schmidt : +\[ +\mathbf{v_j}=\mathbf{P_{\bot \mathbf{q_{j-1}}}} \mathbf{P_{\bot \mathbf{q_{j-2}}}} \dots + \mathbf{P_{\bot \mathbf{q_{2}}}} \mathbf{P_{\bot \mathbf{q_{1}}}} + \mathbf{a_j} +\] + +Ainsi, nous observons que dès que $\mathbf{q_1}$ est connu ($\mathbf{q_1} = \mathbf{v_1} / \|\mathbf{v_1}\|$, $\mathbf{v_1} = \mathbf{a_1}$), nous pouvons appliquer $\mathbf{P_{\bot \mathbf{q_{1}}}}$ à $\mathbf{a_2}, \mathbf{a_3}, \dots, \mathbf{a_N}$. Puis, quand $\mathbf{q_2}$ est connu, nous pouvons appliquer $\mathbf{P_{\bot \mathbf{q_{2}}}}$ à $\mathbf{P_{\bot \mathbf{q_{1}}}}\mathbf{a_3}, \mathbf{P_{\bot \mathbf{q_{1}}}}\mathbf{a_4}, \dots, \mathbf{P_{\bot \mathbf{q_{1}}}}\mathbf{a_N}$. Etc. Nous en déduisons l'algorithme de Gram-Schmidt modifié. + +\begin{algorithm}[H] +\DontPrintSemicolon +\NoCaptionOfAlgo +\SetAlgoLined +\SetKwInput{Res}{Résultat} +\SetKwInOut{Input}{Entrée}\SetKwInOut{Output}{Sortie} +\SetKwFor{Pour}{pour}{faire}{finpour} +\SetKw{KwTo}{à} +\Res{Décomposition QR de la matrice A} +\Input{$\mathbf{A} \in \mathbb{R}^{M \times N}$, avec $M \geq N$} +\Output{$\mathbf{Q} \in \mathbb{R}^{M \times N}$ une matrice orthonormale \\ +$\mathbf{R} \in \mathbb{R}^{N \times N}$ une matrice triangulaire supérieure \\ +telles que $\mathbf{A} = \mathbf{Q} \mathbf{R}$} +\BlankLine +\Pour{$i \gets 1$ \KwTo $N$}{ + $\mathbf{v_i} \gets \mathbf{a_i}$ \; +} +\Pour{$i \gets 1$ \KwTo $N$}{ + $r_{ii} \gets \|\mathbf{v_i}\|$ \; + $\mathbf{q_i} \gets \mathbf{v_i} / r_{ii}$ \; + \Pour{$j \gets i+1$ \KwTo $N$}{ + $r_{ij} \gets \mathbf{q_i}^T\mathbf{v_j}$ \; + $\mathbf{v_j} \gets \mathbf{v_j} - r_{ij}\mathbf{q_i}$ \; + } +} +\caption{Gram-Schmidt modifié (MGS)} +\end{algorithm} + +## Test de l'algorithme + +Nous comparons notre implémentation (voir la fonction `mgs` dans le code source associé à ce module) avec le résultat d'un petit exemple issu de [Wikipedia](https://en.wikipedia.org/wiki/Gram–Schmidt_process#Example). + +```{r} +mat <- matrix(c(3,2,1,2),nrow=2,ncol=2,byrow = TRUE) +qr <- mgs(mat) +all.equal(qr$Q, (1/(sqrt(10))) * matrix(c(3,1,-1,3),nrow=2,ncol=2)) +``` + +## Comparaison de la stabilité numérique des deux versions de l'algorithme de Gram-Schmidt + +Observons comment de petites erreurs de calculs se propagent dans le cas de l'algorithme de Gram-Schmidt (GS). Mettons qu'une petite erreur ait été commise sur le calcul de $\mathbf{q_2}$ qui s'en trouve de ce fait n'être pas parfaitement perpendiculaire à $\mathbf{q_1}$ : +\[ +\mathbf{q_1}^T\mathbf{q_2} = \epsilon > 0 +\] + +Quelle est l'effet de cette erreur sur la valeur de $\mathbf{v_3}$ ? + +\begin{align*} +\mathbf{v_3} &= \mathbf{a_3} - (\mathbf{q_1}^T\mathbf{a_3})\mathbf{q_1} - + (\mathbf{q_2}^T\mathbf{a_3})\mathbf{q_2} \\ +\mathbf{q_2}^T\mathbf{v_3} &= \mathbf{q_2}^T\mathbf{a_3} - (\mathbf{q_1}^T\mathbf{a_3})\epsilon - + (\mathbf{q_2}^T\mathbf{a_3}) = - (\mathbf{q_1}^T\mathbf{a_3})\epsilon \\ +\mathbf{q_1}^T\mathbf{v_3} &= \mathbf{q_1}^T\mathbf{a_3} - (\mathbf{q_1}^T\mathbf{a_3}) - + (\mathbf{q_2}^T\mathbf{a_3})\epsilon = - (\mathbf{q_2}^T\mathbf{a_3})\epsilon \\ +\end{align*} + +L'erreur commise sur $\mathbf{q_2}$ s'est transmise à $\mathbf{q_3}$ qui n'est plus parfaitement orthogonal à $\mathbf{q_1}$ et $\mathbf{q_2}$. + +Pour l'algorithme de Gram-Schmidt modifié (MGS), nous avons : + +\begin{align*} +\mathbf{v_3}^{(0)} &= \mathbf{a_3} \\ +\mathbf{v_3}^{(1)} &= \mathbf{v_3}^{(0)} - \mathbf{q_1}^T\mathbf{v_3}^{(0)}\mathbf{q_1} \\ +\mathbf{v_3} &= \mathbf{v_3}^{(1)} - \mathbf{q_2}^T\mathbf{v_3}^{(1)}\mathbf{q_2} \\ +\mathbf{q_2}^T\mathbf{v_3} &= \mathbf{q_2}^T\mathbf{v_3}^{(1)} - \mathbf{q_2}^T\mathbf{v_3}^{(1)} = 0 \\ +\end{align*} + +Il n'y a pas d'erreur pour l'orthogonalité entre $\mathbf{q_3}$ et $\mathbf{q_2}$. Etudions maintenant l'erreur commise sur l'orthogonalité entre $\mathbf{q_3}$ et $\mathbf{q_1}$. + +\begin{align*} +\mathbf{q_1}^T\mathbf{v_3} &= \mathbf{q_1}^T\mathbf{v_3}^{(1)} - \mathbf{q_2}^T\mathbf{v_3}^{(1)}\epsilon \\ +\mathbf{q_2}^T\mathbf{v_3}^{(1)} &= \mathbf{q_2}^T\mathbf{v_3}^{(0)} - \mathbf{q_1}^T\mathbf{v_3}^{(0)}\epsilon \\ +\mathbf{q_1}^T\mathbf{v_3}^{(1)} &= \mathbf{q_1}^T\mathbf{v_3}^{(0)} - \mathbf{q_1}^T\mathbf{v_3}^{(0)} = 0 \\ +\mathbf{q_1}^T\mathbf{v_3} &= -(\mathbf{q_2}^T\mathbf{v_3}^{(0)} - \mathbf{q_1}^T\mathbf{v_3}^{(0)}\epsilon) = + -\mathbf{q_2}^T\mathbf{v_3}^{(0)}\epsilon + \mathbf{q_1}^T\mathbf{v_3}^{(0)}\epsilon^2 +\end{align*} + +Ainsi, comme pour la version classique de Gram-Schmidt, une erreur de l'ordre de grandeur de $\epsilon$ est commise sur l'orthogonalité entre $\mathbf{q_1}$ et $\mathbf{q_3}$. Par contre, il n'y a pas de propagation de l'erreur pour l'orthogonalité entre $\mathbf{q_2}$ et $\mathbf{q_3}$. C'est pourquoi l'algorithme de Gram-Schmidt modifié est numériquement plus stable que l'algorithme de Gram-Schmidt classique. + +# Autres algorithmes + +Notons qu'une autre approche très souvent employée pour le calcul de la décomposition QR est la méthode dite de Householder. Elle assure le maintien d'une quasi parfaite orthogonalité des colonnes de $\mathbf{Q}$. Il existe également l'approche dite des rotations de Givens, plus facilement parallélisable et permettant la mise à jour de la décomposition suite à l'ajout d'une ligne à la matrice $\mathbf{X}$. + +Par contre, Gram-Schmidt permet de calculer la forme réduite de la décomposition QR, alors que les algorithmes de Givens et Householder calculent la forme complète. Ainsi, Gram-Schmidt n'est pas tout à fait aussi stable mais peut parfois être plus efficace. \ No newline at end of file diff --git a/13_puissance_iteree_par_blocs.Rmd b/13_puissance_iteree_par_blocs.Rmd new file mode 100644 index 0000000..c838040 --- /dev/null +++ b/13_puissance_iteree_par_blocs.Rmd @@ -0,0 +1,72 @@ +--- +title: "ML 13 Puissance itérée par bloc et SVD" +output: + bookdown::pdf_document2: + number_section: yes + extra_dependencies: + algorithm2e: [ruled,vlined,linesnumbered] +toc: false +classoption: fleqn +--- + +# Principe et application à la décomposition en valeurs propres + +Nous avons vu précédemment que la méthode de la puissance itérée permet de calculer la plus grande valeur propre et un vecteur propre associé pour une matrice carrée $\mathbf{A}$. La décomposition QR est au cœur d'une approche élégante pour étendre la méthode de la puissance itérée au calcul des $S$ plus grandes valeurs propres et leurs vecteurs propres associés. + +Une itération consiste à multiplier $\mathbf{A}$ par un bloc $\mathbf{V} \in \mathbb{R}^{N \times S}$ de $S$ vecteurs colonnes. Si ce processus est répété, nous nous attendons à trouver le même résultat que pour la méthode de la puissance itérée : chaque colonne de $\mathbf{V}$ convergera vers le (même) plus grand vecteur propre. + +Mais si, par une décomposition QR, nous forçons les colonnes de $\mathbf{V}$ à rester orthogonales, pour une matrice $\mathbf{A}$ symétrique, le processus itératif fera tendre ces colonnes vers différents vecteurs propres et la diagonale de $\mathbf{R}$ contiendra les valeurs propres correspondantes. + +\begin{algorithm}[H] +\DontPrintSemicolon +\NoCaptionOfAlgo +\SetAlgoLined +\SetKwInput{Res}{Résultat} +\SetKwInOut{Input}{Entrée}\SetKwInOut{Output}{Sortie} +\SetKwFor{Tq}{tant que}{faire}{fintq} +\SetKw{KwTo}{à} +\Res{Vecteurs propres d'une matrice symétrique par la méthode de la puissance itérée par bloc} +\Input{$\mathbf{A} \in \mathbb{R}^{N \times N}$ symétrique \\ + $\mathbf{V} \in \mathbb{R}^{N \times S}$} +\Output{Les colonnes de $\mathbf{V}$ sont les $S$ premiers vecteurs propres. \\ +Les valeurs propres correspondantes sont sur la diagonale de $\mathbf{\Lambda}$.} +\BlankLine +\Tq{$err > \epsilon$}{ + $\mathbf{B} \gets \mathbf{A}\mathbf{V}$ \; + $\mathbf{B} = \mathbf{Q}\mathbf{R}$ \; + $\mathbf{V} \gets S \text{ premières colonnes de } \mathbf{Q}$ \; + $\mathbf{\Lambda} \gets S \text{ premières valeurs diagonales de } R$ \; + $err \gets \|\mathbf{A}\mathbf{V} - \mathbf{V}\mathbf{\Lambda}\|$ \; +} +\caption{Puissance itérée par bloc pour les vecteurs propres} +\end{algorithm} + +# Application à la décomposition en valeurs singulières + +Nous pouvons de même adapter l'algorithme de la puissance itérée au cas du SVD pour calculer les $S$ plus grandes valeurs singulières et les vecteurs singuliers associés pour une matrice $\mathbf{A} \in \mathbb{R}^{M \times N}$. + +\begin{algorithm}[H] +\DontPrintSemicolon +\NoCaptionOfAlgo +\SetAlgoLined +\SetKwInput{Res}{Résultat} +\SetKwInOut{Input}{Entrée}\SetKwInOut{Output}{Sortie} +\SetKwFor{Tq}{tant que}{faire}{fintq} +\SetKw{KwTo}{à} +\Res{Valeurs singulières d'une matrice symétrique par la méthode de la puissance itérée par bloc} +\Input{$\mathbf{A} \in \mathbb{R}^{M \times N}$ \\ + $\mathbf{V} \in \mathbb{R}^{N \times S}$ peut être initialisée avec des $1$ sur la diagonale principale et des $0$ ailleurs.} +\Output{$\mathbf{U} \in \mathbb{R}^{M \times S}, \mathbf{V} \in \mathbb{R}^{N \times S}$ \\ + $\mathbf{\Sigma} = diag(\sigma_1, \sigma_2, \dots, \sigma_S)$ \\ + telles que $\mathbf{A}\mathbf{V} = \mathbf{U}\mathbf{\Sigma}$} +\BlankLine +\Tq{$err > \epsilon$}{ + $\mathbf{A}\mathbf{V} = \mathbf{Q}\mathbf{R}$ \; + $\mathbf{U} \gets S \text{ premières colonnes de } \mathbf{Q}$ \; + $\mathbf{A}^T\mathbf{U} = \mathbf{Q}\mathbf{R}$ \; + $\mathbf{V} \gets S \text{ premières colonnes de } \mathbf{Q}$ \; + $\mathbf{\Sigma} \gets S \text{ premières lignes et colonnes de } R$ \; + $err \gets \|\mathbf{A}\mathbf{V} - \mathbf{U}\mathbf{\Sigma}\|$ \; +} +\caption{Puissance itérée par bloc pour le SVD} +\end{algorithm} \ No newline at end of file diff --git a/14_geometrie_ridge_svd.Rmd b/14_geometrie_ridge_svd.Rmd new file mode 100644 index 0000000..399535c --- /dev/null +++ b/14_geometrie_ridge_svd.Rmd @@ -0,0 +1,54 @@ +--- +title: "ML 14 Géométrie de la régression ridge et SVD" +output: + bookdown::pdf_document2: + number_section: yes + extra_dependencies: + algorithm2e: [ruled,vlined,linesnumbered] +toc: false +classoption: fleqn +--- + +# Coefficients de la régression ridge en fonction du SVD + +Exprimons le calcul des coefficients d'une régression ridge en fonction de la décomposition en valeurs singulières de la matrice des données observées $\mathbf{X} \in \mathbb{R}^{M \times N}$. + +\begin{align*} + & \mathbf{\hat{\beta}_\lambda} \\ + = \{& \text{Voir la dérivation de la régression ridge dans un précédent module.} \} \\ + & (\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^T\mathbf{y} \\ += \{& \text{SVD de $\mathbf{X}$ : } \mathbf{U}\mathbf{D}\mathbf{V}^T \quad + \mathbf{U} \in \mathbb{R}^{M \times M} \; , \; + \mathbf{D} \in \mathbb{R}^{M \times N} \; , \; + \mathbf{V} \in \mathbb{R}^{N \times N} \} \\ + & (\mathbf{V}\mathbf{D}^T\mathbf{U}^T\mathbf{U}\mathbf{D}\mathbf{V}^T + \lambda\mathbf{I})^{-1} + \mathbf{V}\mathbf{D}^T\mathbf{U}^T \mathbf{y} \\ += \{& \text{$\mathbf{U}$ et $\mathbf{V}$ sont orthogonales : } \mathbf{I} = \mathbf{V}\mathbf{V}^T \} \\ + & (\mathbf{V}\mathbf{D}^T\mathbf{D}\mathbf{V}^T + \lambda\mathbf{V}\mathbf{V}^T)^{-1} \mathbf{V}\mathbf{D}^T\mathbf{U}^T \mathbf{y} \\ += \phantom{\{}& \\ + & (\mathbf{V}(\mathbf{D}^T\mathbf{D} + \lambda\mathbf{I})\mathbf{V}^T)^{-1} \mathbf{V}\mathbf{D}^T\mathbf{U}^T \mathbf{y} \\ += \{& (\mathbf{X}\mathbf{Y})^{-1} = \mathbf{Y}^{-1} \mathbf{X}^{-1} \quad + \text{$\mathbf{V}$ est orthogonale : } \mathbf{V}^{-1} = \mathbf{V}^T\} \\ + & \mathbf{V}(\mathbf{D}^T\mathbf{D} + \lambda\mathbf{I})^{-1}\mathbf{V}^T\mathbf{V}\mathbf{D}^T\mathbf{U}^T \mathbf{y} \\ += \phantom{\{}& \\ + & \mathbf{V}(\mathbf{D}^T\mathbf{D} + \lambda\mathbf{I})^{-1}\mathbf{D}^T\mathbf{U}^T \mathbf{y} \\ += \{& \text{Soit $d_j$ le jème élément sur la diagonale de $\mathbf{D}$, + $\mathbf{u_j}$ et $\mathbf{v_j}$ les jèmes colonnes de respectivement $\mathbf{U}$ et $\mathbf{V}$} \} \\ + & \sum_{d_j>0} \mathbf{v_j} \frac{d_j}{d_j^2 + \lambda} \mathbf{u_j}^T\mathbf{y} \\ +\end{align*} + +Remarquons ainsi que la décomposition en valeurs singulières de $\mathbf{X}$ donne les coefficients de la régression ridge pour toutes les valeurs possibles du coefficient de régularisation $\lambda$. + +# Régression ridge et géométrie + +Observons la relation entre les étiquettes prédites $\mathbf{\hat{y}_\lambda}$ et les étiquettes observées $\mathbf{y}$. + +\begin{align*} +\mathbf{\hat{y}_\lambda} &= \mathbf{X} \mathbf{\hat{\beta}_\lambda} \\ + &= \mathbf{U}\mathbf{D}\mathbf{V}^T \mathbf{V}(\mathbf{D}^T\mathbf{D} + \lambda\mathbf{I})^{-1}\mathbf{D}^T\mathbf{U}^T \mathbf{y} \\ + &= \sum_{d_j>0} \mathbf{u_j} \frac{d_j^2}{d_j^2 + \lambda} \mathbf{u_j}^T\mathbf{y} \\ +\end{align*} + +Nous remarquons qu'en labsence de régularisation, $\lambda = 0$, les valeurs estimées $\mathbf{\hat{y}}$ sont les projections sur les axes principaux $\mathbf{u_j}$ -- qui couvrent l'espace des colonnes de $\mathbf{X}$, i.e. $Im(\mathbf{X})$ -- des valeurs observées $\mathbf{y}$. + +En présence de régularisation, $\lambda > 0$, les coordonnées, sur les axes principaux, de l'estimation $\mathbf{\hat{y}_\lambda}$ sont de plus en plus contractées lorsqu'on progresse vers les axes qui expliquent de moins en moins la variabilités des données. \ No newline at end of file diff --git a/images/SVD_reduit.jpg b/images/SVD_reduit.jpg new file mode 100644 index 0000000..d542120 Binary files /dev/null and b/images/SVD_reduit.jpg differ diff --git a/images/hortensia.jpg b/images/hortensia.jpg new file mode 100644 index 0000000..377b11f Binary files /dev/null and b/images/hortensia.jpg differ diff --git a/images/hortensia.pgm b/images/hortensia.pgm new file mode 100644 index 0000000..7605460 --- /dev/null +++ b/images/hortensia.pgm @@ -0,0 +1,118 @@ +P5 +256 256 +255 +/Ù-8\bYZbe{py|xwgce~e[WW_xlnzW:YMMJMe{v{}dRNTKFN?>7[{6?~ChpiZ1:+(~ɪ=+ 15[J2A:;Wۨl^rkТ779szoommlkf`^g4!! 5=3(!!(/367=9- "!,.)%%"+&!/".-  i8!^Ƽvx~yxzs{xeaMFhsq|gOHidicc`uuiOCBBA7KvEITE-7@H:Bun14'R¢1+*Aa10¢NL`ѫgbջ]8Dmxm}}xqonjcA& &!!"7>9'""&*,274.( ''*11'()! #+/)%% ˰ɩ% +!I~zttylmfoYJVuibfk`a\s`}RG?Ig{E];y#)=8&4S{M0$2AõN1'X:&^ye5CQxƠe@M;=?Ե¼;4nr[fgljdm{|t"!)6==*#!!34.-.$%((+0/%&(" !"! ĝЯμnùuv~y|oh}i]bxbAPzofkXNK>qa|~ol`Gq6++-6A81u%).5Ÿ;QZ|B.'iC"3@6űŬwCLghl~~^,Li|}sb<%!07:@2&"$76/+.'#+$#%),*&$(9)( + NİѢ̿ڽA\ϱ»tjeqokquZjycUM\y~lgLWKQM3h_isvtQsGy( 01/6--;EGΒ̛īȟ,;+%'00>AijݹɼE7MJ_YW\''04('($)%20199+*,&-31/.$%'%#"&*-,!64+>N=0  +  ðê¾ͳǸ}o~ms~}Tnp``ylesx|{igTZPUO4c_TnjsНYjks6%%-#/0Yi*&C<ϧִܿ?F=,Au̍ʵ<*=ܼȻ/3 %PI07'%&# *%!'!!(:1+.(&,23%)&%%!@`dYYVG0 ŸṦƓþzpdkovXsxcglg}hyhvunlhYJFNG<:/+@61265.-2+88%00'(!!&,+arhSND>H7-  DzʡYqyjkYsnzzjlxsj|bFcb_aQBJVURWITnQlԜs$#%-[sgcaViikso?-:iϱFV+)H_1,0?85&()/3 #()&"(<92.04-#$$!.8* )30$$&BbhcWKFGC:EE'@ϑ~mvnsfdhmxxxmcytkgIw~y~dlaULNQSRNB5Czɇ{WPpŒR01W}q_oq}ht\T73%0=l`KjKZ$-/0KHE69878  #&(('& !!-8==:61((62#"!%--#!"*I\__ULIC=65$ !"##%'+(#!$#!)6AHKE=60*,4,%)&$:ZcWXIDHE9:H7'(Tÿ¾vplun|ne{k~aeTk|^p`\GUaONSCJ49bz~XFFRwh`tf6pzBe<#1*c~ܗ>d9".D74;29zεʩ&hm%/2)Hje[M\|Σ°!!#(.+<48;,*)#$##*))**-7ANKB984.#$!]~TMWamlsyo~t}mynJdtxkwgsXRV`T@GN>?6;P~mRB4:19RUnm`S||íoW93!#45*%''"$+9A>BMMF>BA=<=3%)!,3"$SϮ¼ZXKRXhfxpsxz{jttg_xWpiKTTSOGC?9?GPR|~rxHF5:05DKbŬ^O\wƱja-n.,7{O/1./+;wRCɡgƘ%NB00$4;<><5-*)%)7/"KʱĿR8Y^M\]aw}wjuufde\owsz|sdXXPGKOA252GWNvxzwzh8I;6-RVeőUPaam±Yf}0rT;>B4"%5*#"ww׸üȸzPKMRTNJS`vnje_Y_ZW`q~}y|_y|WiWH?CHA60QF>CWl{}vqlC+B73)c[pƠRYg\ƪoa`*l>[yRD"* q|3]aptakypBpx+ >64/"!,/)%)( )661/.+*$'-/07=#4+3SH=:>AA>=>;/()($ ī»ͧzEKWNMPRPOZflg`[YWUSNMOZq|~{j{wXaUMD>738AIPJ9?Znqtkke8):3,' QPqRh^^xΨ``l,a8EC4H&]PRIOLKMPT]fXWTPNMPRCIMUkwwxsy}r[UXN@5/2AOJOH5,Ad}saJ.*10%'2;U}eh\jtùĸn_g|7]b;AB>H3v~q{}y9#58. )+/.#&))$"'+.,'.!0=4(28:94DBCA?<;?B62*#"$$ ~íɼ͡WL88;PPMJKOSX[UWUPIHMRGSWT_rwntuuirbLUG938CKPRC;7++X|_M@+0*,#+E/H•abywlYtzPfJE>**_jvheW]Yn&%>* )*0- ""#!'*(#(#$-5<5=HJD@I?@@>;9:;*(&')%kǪαYMN^YTONLMORQORTSLDCFLLY^UUcoqzus^}sgLB?<=BIKLQ?9:,!@nQ84*4%(#.L*$dͳw`w{dgihlsoo6A0FW~albdzfg{Ķ4!+~) //#',+36-+..'%$&*1/'%'$#(03(+7<>@@8:BRRDDKFDA>=;74'$$(*#q*½ɸ{FP=N?OPONORPLMNNIECHNKTZUU_oz}xp\wpdL2ANKB?BFBEB5%&0db|xF+')0 %(/A)!:g·Ugy{wsza^W\xl(97mztzt_f5;3!$%#"0..)#*4?HG@82/7FLLHCA?<52/)%%%Ge;K|ƳuFYEWFHLMIJPRPLMMKHJOTMMRWY_m{utj[iaTE=MYPB:=0$!!?47;/'*,(+ '310,2*1`ĚVXfjy~quabZ^|b1Ow}omΩ}N/.:!*),+.BLC:=>;2)!$2>?<<92-%#*267:G;1+&'6KXVOD5&184HG6DIIIFA;5323-++&ds tvQXk¼aOJKNL@FFBDNTULLJIILQTTLMZ^_frjndYZQB=WZXMDA;67.,/*!(9721+%/1$))#+=4$/.0%b`U^jy|Wcae}ZDprmll{İshO;'+0*'&-CJ<77776666&6BA=>>;C::CMRW]UKIRXYam^XOGBADF:<1BB3CDBEF>4-.10.24-!%"}s!(utbxg]QN@LQZPCA?@=9B<@LROPWYNLZeaYVPRQMLOPNQI>67=@A:3----3<40../* *z@]orxOSJ@F7.KOPQQKFFIBISXWROMPJbiujtl^QHbcoXCl_H@IF732B"6AEF1761382$%!,4,3(ƸλeXkmhcbhupdkuqfk|{q{ЮteZc@;)*,5%(%"2=D?77?F=/%-211589;752)#-=OSYWJ9239988GID7/1%%*(*.2-)0?JC7/9<3'&*/)+,("7PVR?* #"!"*.'(-04;C93/-+*,071)" #(+& %++**e3u_bVd;RHh{p9N/( #$'1;?MJQ]`UNPJUq{qb\RIXmvklh^ERI:?B5'=@G=A+fY*G; ?Ȼ־¾fr}dgpz}nvms|ҤڻùƳokslpuga[BKSD7+%%)/0)&)+('/BQJHOQPRY\fsqTHEBQkr]^cSESJ45:/!*+=N12pJ9,5ũοf~vnhzx~{⳴ɿqdirqvym_f\YM**2$1-##-4EA=94-%" 5GRT:1!"8B?CEEBCHHE-&"*10-,7640,'$" %%%$$$!Cb\GFhjks`aEJEK_ykBF3&+ " <7>H&,_|:J;жӽºȽakwxrmnpowίӶZp|wrprsfg_\U)'+--($+;B>E@<4'"$:MN=,%,05=AA=@@=>BED.09:;<==72.)#!$$¼°ɺoecpywuxfpͯy[bsrqksn]MNZ0-*O:'*9B?92,&%# .CLG9($+7@B@BGKKJHJLJD.':<53>BB<83,$ ""&'"tvp{dVtXBOSG8-'  !CK>4HSNNZ^VaXCF6OLF`JRLM5CG8BBCPRHEN-@}ȘMȽɿоmddkfZ]kumcg縪̱~v^cfvjx|mpeWKEN%2=D<<=-26*!  "% :GL9 '8?76@FHNVTRMKMMD8!>TP@7=BD>84,$%+!%&$(-) aKs[~qJOX83)  +CKW^^ZVUdHIHISMN[YJH=+6;RVQFGQOC@6BxxǽǵѺgǾߣarca_as[AGLS`v}hv}eVLEDCA?BB,# &%#=D4*+)'.:6-&%"$.:LWWMMXUGA;>LSKB=?DA4)$ ),'*+(# ޖjVDHGPeCB0 c|koHGNLJHYTT><>7GTWTG:;GIAD5\ĥpvb`j^kx\tfipNTMFBEHGE6% $  #2;;2*$$$##& "#(5DOQIMD6:JJ=7<9)*.20.+'$?zPPAPAR_Ftagh¸7 aZYQLHVEYPXbNFF[UI@CMMD/6A@~ŚVWKDDJ9(Eu~ƤiohdlU7SQG@EONE0=E2%" !5A=3((#   #+14HPC:>A;6:7(*062-,+(!m~~co6URJflwgkpvg}"{jXLNEEF\LJP@>CKHDEJLD9E8@&Ҹcc5\lOKE>=Pcp̖~op}m]`7=TTJ?CRQC(7I? %9F4&+%+%"    %4DG>9<8:7-""))$ %-,+,-+&"l|lI;VLG\gAmv|%QbG<540=TOQI@DRQC2$" &+%2JE'$*)! + + "'(*>FB>750-*((*)  #(,.+***xn{gfTR`-ZVisnpvo# A-BGJDAI:TKFHH@<=D:4=2ļTPO5FD.BSBRX^gĺxyX{fiWDGC>@IHA0/.  0.;C;/' *-#"! +&)$.#,,& !C.8.(+("%.,3T:%+÷~ӆED5C@FQORRNSX_¾n@lynG::9620/-&3$ &#C\J/).2>50//.07#*/245XC@UT^uu~U+~|LMLD;DjǴx +991 )I")"!(&$042'1fůȸºfί2EGJICWODZLV`fiʹ|BD9cA;400-%%"+FxU/"(0,#/"$+! -?G6-:9.-5830/,-6C +   #0897_MGKA^K^wopwkedse?Rhµy~l ,%*85B4$$#" !)1=4*![DǼĺĽ̜9HM]IQY^\[]_hkifuǫzg'9A6.3+"*+')Jee7 +76+&'%%,0* &%##3GS= &1545<@P_bT@1) )27ulF[XIG9( './,)llR>EO2.[aSrxn|<%#5;-($"A2$$"&*&!,?=0!O*K<1:1vrŶ¹ϸF+*/07fhkjlnu{zuv}3,7JQD-  $"!:VS:0)*(!%-+3$#-,)()-+ +() (4>;403/594( "'"!-8<9/'qlz}L-QN{*Uqbb`xqnul& 76"A3&&$'-&*A8;!b(D?>).<#3Dŷض,.) .1'>Cagouy|~{sqw~ĽQ3OQIF/!$#3D:'0*(')58(0" !%*. &5<;400123-% !'.331.}wxAqzwks\WZAMRrj{yeo|q 91 +B2&*'),$)?34$4DFC/2*$AƵȻڴ). 4*>gbfmsy|}xx}soóĻÀ+GCF>1,.(!#*23-+)-:9&(")/'"-08?1 ,21/-+-,)% "$&yE(=t|{[VYXArzuy{wsav}A)<. C2%--,("".=;214)B<@A7"$̳ɶԴбM( +*Hdkjiimsy|xlh~ĵӺ*GAG>=?=4' 1! *36;5+18/*"%/($0;>58>2!! &++,//-(%   qM72$rhuyO(72ʏost}~! 3=) B1'/1-$!*7?B99$h-EFKB.6)oɡȱ&/,Nclkghkpuxusu|ŬJ+IDBDIJF>:0&/&(#!#55:8*+21+' "+=EC,*,& ./32,&',.+#   CI"9XQHN1%/#i{t{;@ 7<%@2*04.!!/>A?7+%x OlDIH].Vmȿƽٷ!-'4fffjnruutrmйѶ!6=?3# 1;;7,+.'')%%%!$(&(7CD?#"" *54;4&$%! !"   +74-30wC73}R,0'~vyg  8/#-: 765/!4X9<6")xTL=2Ŀ-ѿ%$ %%-Sjhtaw|{y~Ʒĉ5+DHIE'81;%3! '=NC*+@1'( .3#*A;@:,01%*49%# &!  +"[G1#*VerVuvfp#HajBS^s*   A63%##",*3=2!;I@Ndi?0k\H<}u$qBοð!!!%/AVabfw}ʂB,+\59Zx: #/46- 1<7.*))5;'&#'*3LA?55062580(!'' $*    pX:/Xeuxit17IlO=$${~$uG (2&%*%)6F%-:,Kg2,-(RH9///& 3R]`eMPOJGIJJ@pXZIcs01('xʼǹm)L//)a3*.R8&3.V4>k`T.xulPyAvu{>W@fygЦ+$+N$"578*$(*,IJLQappg]TJB==BKU\zv /X$.MtŏdzbpKGYf]d*/ZT%0I,1)R2%t}ZOCO7eT?w&OKxaaq6<1,,$24>A4/'->.',M %CRRHK]jcVWdqlkeYMJPWwq^PbpDvɺŹɯԻ @eYjV3.:,G7G?/(.6v^qgKpoxq\TDJzq`JEvUEDK/.Lpbk,Wq[nOfOioknlpwhyzyuİΦm}oB%!AXTKWSIA/+AqL30'06E=/2+);5) :78D[kgYNaoj[SUYVXVOIJQYGiguP̵ɶ°κo^22(60%1. $M\V`pk\O]gbUMHEEHIGFGOVqwf >ŽŵƷìɜܿ#YJtaQYYIlC1-<*$%xqOA9*J>p3KmtpXtcVz6N1/;*NE*Phs'y5~}|}Zrikya|«oB$>4.)$'!2:0.54.)%4/$) !Qg^fsqiOTUNGFED@CFECDINkV)GʶʪħȸD QW!f;@1VhZ_e^3/OGK.[Qw=.\ygzw/do|_:x4.=X0-p?6{QQUzYx|h|gwl{{{¿űɰür[B"OS0 #"#+8,334+ # 2<3041*-'+((0  5NZeifbMMIDCEGFBDEECBDGfhxźźǏJxSS.sF&`lci;]L&t&5$';67:?=XIvqGX\x.SBaTcGZq` ?q?{N%l6=X*Ĉxauaerfoyvó}zwmWX7$! #+0&"00&'!%-+.7645(%#,;&  #BdeYRHIKIFFHFB@BCCBBDFǺWͱȾųλŹɼ)t>D%\7Ntg]dw`c,I$2o~=io'!@:EQ͊bObiU}JZY{JT|Z'Mjn#wrW;l`bgf[}whT5M÷üùyxqfwQ7-)' +"*(,+&)19580#'&)8) '3'PziWUFEHHEEGC>CDEEEFHI1VûͿſjGLQHdXSP ]Lb;k'LX`e[#118?F9&a79wMv#UL%WIR4cnVE3D"F7S>BCAIIIIIIJKƱN¿TE8BB9^aWV0[ZLOAuR]EvjH'H+ADe6=pQ{S-DT,3HmdF}~*}vfoXvvWwySec|p[^9:DDm7PTA,#,Ľ}wonkcC(%'!6$& +7:4"$'#%+-'6H7!%B]OBJDDFGHHFDC@BFEBGOôsǿͩJ&M>LQJnY9C9@].K6YUU*oj{Y7F.6Q9Q[NRUM1b5c]oEfy~.39EpdZ"A)g–zb!oq|`_VT_WA4/'YUW5I/)ȽĿx|yyvrpZoS "' *#043* )!%#5>1 44BKK@FAABCEECB>=@CA?ENúqQ4$='NWRH'Ia[V3;'-6'qpsU1/#\j{,qzjg9Hg{Cewmp5AVfR2PTJhlqaOW]TM]+RAE6_S_{ǧľҰvkuU# $,)$!" ;*0663,&$ !0.6<178=7KDM@@ABCEFFCDFFBBIT~UI+2$/)`"8E}~FZJU)zI_wo{_{hC=>P^{U%B;S3=fNcCú˻xlmC!7$!# $(5&''07;<:1'!$#:C@-+29+D=LDEFFGIJLHJLHCFQ\¾K82C@f*@dqn[V-=7YT;/B9DCIJFDAACIKKP[fkvijn}K8db=TcF?bnAY2I6:(7FFDCCIOUX_giehu˰mwpsp>a>iMI+=O@]AV`Wd4282?n(\udC$|j}dQllpouzCTtPı\XvX.+[όz\Uo7/1+eW6J/1,AOXvqj\n(ʡ{/QNPdU)&u`EHh|yvz~{`ABQUKBFSL`ŽB>RKH7"&%" $39' ",yC2" 3BIKKMKS[^df_U|ЄcQz}}wge?INSU<:@N/\XI=\cxt+>-309lJ3OfbnA]lT|Q2$QrgQE?J Tz^L{+53=#ASeXjkxr[sJa\PSPb^dXyCMEPHW}zwI7A̵ǩoGIU>  %% .A155!$@>%AQ43"!5;ATISMQ_if^[ҪmudIE9KCf}mXfz+&*A@5;_:EQX0]SK>&7BVG2YHJ|YiQ>DN-IsQZoF42D5:N\z~~zhLBFRTFH,&)ɾžϪJG;!#+$(-*98 !'$ AP32)11FKNOUahe]VõalRUggO4PEFGpALG|woooob[{E+Jd4/XD=B@ZdZ0$ 'I9HqzfCiHMKF$?;+8.LTc[@bXhkgXZ{bbYj^Uc[depYvU\"O^T^D6pztw~tiAIBDEBF?*""X˻ȳǣ˵h.A#/- ?90"=L  # !%CN01 #))939=DKS[]ZeelgG>TB<<7P1;#ukjirz{me=$}I;63XGLjOI~RuM_b}#R}Vy^hwRXgTOjµda^QBNP>x{VW}[KZYLS7BTE;zhym.N|_]oj[djwgA}fXI_N>vqouzpS71G;:87F63%/2ʽǽɲǫ*B;uY<#'%05&#$',)FM-1%3F/.1119IVZù!^ruli|[rYeZ|?{Rsqd@xKW%/sV;)@'3;UYoxRqn1/>TGo&bUmsRK¶{~U=WIvYPgb{|gOGK8|@73R`~Hu>ah=jAMpxQ}blY]Al|xd1EPcso_K9,2D47::F/1#+3ļƹã´M@=2 4'63,% &.IM+0!(#2OA77304?HMĹſ!5[gXAt}IOEpL=rI(TUoy^G\DETSdl{PcaiYaoic36FNI@<7166,08:;->ŵ¹δ|Ưö("+30!,JN+/" .(!34867>BCBD2 +;=rvqiSh/*A@D@v\vxluinpgQco6/HBQH07#W_dmY/0J>5Hq#9'Y}0Jwx¾Ȧ+H3E?iYuSGGI-@Sf+XPsn\e@ijWBB32,'/?D>A*%%*/1B")KūΊàE/$%77%%FO*- 1/$24=59KRIDGĽֻ)8Ixwxj7:SdhcsIaMf_M6B4xP\V4"k7dĢwWNaoENDcy_icGV>efTP96`gf{FOdbrO_@FY|lBXui;\F\m\F>M/76+);OYQ*&'2_I@9B̸Ŝ^,9O8 CO+, %:]]UDATZKDJǵжĜ #(&PU}pGJ>7T4meUPcVn[dgHcGNPlMBUmWLqv{vctevL}OAWa]wwgCF6cJu5üǟeM!!,I4/+-','!tJCO]ǿБw 296!EG,194..*2QsgUKFKPOIFſȺͻw1"uDJv}|MZ;KuD#%brnBS|Uz<>URK4z>C]6naPiWiwTYPn]`CvfV,f,T}n(,VY=p0jZ?`^\>C^$XqY0[4A-UGMhlq|ig!ue|D>h^MSrWOF_O3.())251D)8/n?"/GŻ¸c;2-II8'/,+71-DSlj^NCBIKKJа̼+ #' f%0KUou+kBRg_SWcfyfG])`SS16`5F8BJ;AD#X>*&'6?=Dg9QKƾr(/8JPR!*).E:.CSUdibPA?FKLLȼҳ&&#6rF^ebO>}``g4L=f?k^FnV{tbsL{IB_l^QldžN_nuUt:?BU@/FaSWMLM P?QntBP946.zWw`{s`C50O&]3p_qX_EQrZV3tQTVMQO06$:xW#JQIHY:_GalV2C[r%).8PD;EQW`d]NCDLPNLʺͽEg + #7~-xW\I.@*]iSHMXMHP@=wBJxTye/Tm^YZgs[?gs5F dZ^Z]lmdi-9GJ_rSfwC=IUD^wbsYDd7yIRCZJiUF:-JjS ;?@NO=`K?RIM2.$6_-(.:VND=I]cc`QFGMQTVȳѵ#)M<1$T=DAvD\MHj:C7:BH^W\}vpiDVvlbi5G#@4iQwYd\>;(i}{vZjx61BQGr=R=U:/vo>ZKz:d-NZgsQd@DI7&=HQ6!);}9m7:PYZSG61»""/ (Uz4&+3SPOAIbegiZLGJP[eö˺$ +  \LXFI>3h_GT3Fn=Jk:lvWv_39)5B^I9C%LAMHNJN<t`!#&,'B^8*/.EFYRR`_nnbWPNQ]kȽ† + + ygpw=evfkU9QiDkzq|C`CT:asa^HWvkZ^BQecLGK0+RS%)ƶo^b]ZnREw{vqgf6_wi?JgaOhO>P`LVRpI3JSg]Y*N?BQf>/&GXBT`AbP.Z;nIR#'(1G:1608;^b\[VqmgaZUT\fǘ  DX8,FJy}ktwmeFnr~~xWyuvdEDei%Zo\bd\29Ug~`YM|˻wi|Xy{i`JGSgSkR[o}WIOe/0~L*`>=2"XauvfM[O]EH$NȸvTPYdQXaPokmE̱YN2 !%1778856G]jmgcbm^fyoZU[ü#L + #ī{q}jq;gbhkr`=M_^">{9][TH^iw*XzxfPgC]F<lŽ}m6qzSDGKSlŧrqg9r\}:iiŻSENs:bSxyPRYvOUb/IDK|;nGYIuPzil22E|?oGa~slW?FkY?EPO6(0yǭáNUVI:4-*44993//10,4=GWdf`dqvyek4  +NvkhzjfdeAhMyVa8~kmjLduo:8>Jqikj^hyv{|˲yhlj~CsG|)+wg9euPsMlyo=n`6If9s7S`l51ijhOoSe`:>J?]c.75j­IGKNHA;/(0,21,'&%#),-1BWa`e{}pxY +  % k3Sf\h}m=TPyoSlra7]mNC<[V|\|DfNXsn'K|YQ{&j_^tw}~|o^oMh?/rbRFUWJAybutxbi=f<%BMLAl\QHLE`\,8MzGo;r|ryf6tPLe{z3$.:IGIEED2!#$&$ &**/139CPXZbxyvU & Jhci}[BgJ|f?P^Q?LuI\q\}q[Fysa4[=cbesctylt~]{jh[T}a]-|\Bv[bUX7cgos{DGG6GTXpXJaW0K#.pNIKXdR70p,NBg8@i).<Φ@SDMJKGJL5$1886?OWROOYdk}x2 @3o}ahhwO]n_`ZppGrOA]]9I{i]r9u`UtlxyPTwnX`tYxwxx]BNO\Ruy3{\ghTs([k|czBNIP3h5e00K5iZFX4BO1E@3V3zb5W%d7`zj;S?K^=KJFJI0 #-03)2O\OEFOM\|}~0   ("zD2Rm*C@IEuB[Kt?O}_^@4aD lĸ¯OAPRGLCNFOKMMHGCEA' $$"'C99L^y}6. +5  #XulV<)t~osrTqXWXlg@2sJfp\gn{rWtP\~ ;-)gyr}t~=b_NYx__jpc?QSkۄYJ^awTKMt+;TBpmy}`\Yfvrb|ZsPM]mtU&'9<=;1'&,#*#/37.91/-NFE8F@ANTxvu\<_5b 1?"Xhot?pXdXAih^0:mT50|wy}vh}6Q268!,+-*d#&s^JGobdՏ~W|mKLiQd|yPD=[DT|or_RAXBRztoiqr}lR!!&(-+),"%)9?;5<("69O=B;A@FDmxcTdYGVbDso5No]Hwst~q{8@Qml. K)WbkPuhe\GoJ%): 87&6pSgV@1Nyuuhj|kA|t^f]dnB9X@Fx*ezB5O@>!jbSXqzxikHe{qi6!!%"!'10/0.*((*$0EPH>56-1&BBE2?;::gx~~TWIFYJhXIhbI)?G!/}zkquchxc[xYG1 sl5cQ~j}cy~uzw5tu0,**))N|^vI6YVwaoTr^LKiz{qF4TKM恄HM=|~sLh]Y[sjjh|YP:loe{zuxuAdVgup~"&-#1.,,+((,#$ #$6nmyi\7?syd,+,%+.!*3#\dDktPbșMpa]lysViJdo`aoELc]nsprq9Mco"(\V{@]o^e`xqzm,+.wnkik~*7'3($!& & -_D>]UORV}~VECGLKIVVvd{HjptmqKO]x}U`uZoaISN`bt@'Yf}Ru}c/C#&(-%#4ZdYkUnoyldd~j^_ALgn^e[^eIj=8oDw.HBM`O~z~|p?  " +,9B456;?,GQZU/Eppt|o]`[TVYRNWF@fo|YICV)rvqfuev5Sc+vU?d{&a_a{Q/%dkˆLjaT2(SR ) /!M(,=jZc*)hg08RlumsFVFUUklXjllMB\?Wtc~qsuS?X'fdaCM?FNSrc~SMygjfs:k# 0 YMȿP9:'AanmmsugSSWLPTHEhuhYHXspĊY>f-tu]TmXK&3Dc;=@Ad;]?Y[eRYwWPXlktT n^+L ȽR-(Oljmpqwvdc`OLL>MDE?D@iwkt;EKL*WX\_dbx|b|Y!9*1I.;9SnXBE%339I;~q0/(_n,$X|mfb'WB!$}}`/q;w +eMQ{U>tutQlimFKW_e~]JZA2\blV?r×X~snBqYFyoPj\Kc4rz`_L u}T,$/aqhosovybWevV?QEPj`gӗDB<:G rwy^I~\Q1uZA!$?zHM ^eDW:.kUgHQgJ|g^l)'0+MgjG!='-ywDhmW%,5eUYBpz]T?NyNR|bMGDq{w#TYDDIE\cs\]SCJH``1TNUd)zLKe}?C9Cqz}O{wwe_xǕv&+M]6dGykzRo3_hokYg7ra1Mc[kjÚ©½NB?dnnmpkit@^GNUZefURTIpYn)^3,*C?M;R\xƉ=>RF_GUD#71{=M;eEK.+.VMpwGfpT\eaaolm{sOpSS2+/7`*"TR~'`s|Mcfdg4gd=jW}vSTr:-gVMTJC;BDJODYuv}Oh_UYIBHY.d3GU&H}ƻþJKMjfipogkyNEsjKNL]P\RKbkmzwgY+G?JB8:xR^dyr-"(+cSKWA>Z85I:paY`cg|i`pi_k@HXV.I3v7-jd`V^(,zteXeKm`RG]OYlLdU5Vfd=GOgHoGj[E|oSm @~%F/4JĹýGHZdgijjklmIaXjUbITMOPSJ_yZ.@)N9S7IK/q{bXk6/scpgACvK$)X{]čxiysZ5C̺[d_ZgZC,cbNĎ~&[4J8=|ns>gXZqJ)inM0~X Pz0?.4^3.hu4o}Yn=OsHr@yc,):~jsEii3%&+wv1ddl;05n]edeLkVhZ`fU>~FOP9CL>]ich\\lD5\^DS#3)&108EdgijjjjlnlOYZb|VVISF;SIBdw=ZISGEywNb|~}UJBt>>/D"\+rqZHw;*]-ML;nt_bMRs~aZCoVhQ-5j~Zt67!]wpud;b[afruZghY^M4}KE_qnmSR?~2[xiv8W#4@j8^KNH)Unf^xbg~|#fxd=Qpatk!{D#FdfHaZTL[kjeZSHR}Bsw``C6!Pb*6Vd,bJ(}huo\^Zlwhdqu͝*4gW|ikljjjloaa[Li`qmiKWWGLNMW3XZL1_E~[z&[Nbg466IzuXYJ>;hYLbpLd]Cl/1_Z|77N4xj\A-`LTRV,*wtn^13q}lvg^P9EB[9_XGrt]}n(a92H}54@PnP6KPo`fjU;`xm2süúӔKdkSjklkijmpW`Y_GMxnPXTCJORNEbHN:\AFbvuuldU^l+@K9?{;ObGXlQ`|љJ W@1@k>;9W}/h*qraatfZYj`^0$.mv~y~m?[CTuKMNRiBRIiRagk|g"!,T}tEF==E(EMfoB@u0hd2 o)5Ĵë̃WqjlljijmpKjSW]V`ci[VBARKFCVI4GSPgA=_>\prXazwd]s1_w``yTcRgXhy>IgWYrMB1RF82Q:Fqxa[xb[XcbiC\z}D(RZ,,znj>HE]~s3j?FOCQmfVQS1c4yOFkMKRD91m%PiMdq/3i22 @ĺ˻x]illihimpMkTOXSccRYd[D?FEERKbHf-BI>EjfWrQq%?^|mN)$6"sJla~|AM_RZm{|HNdOLK4gq`O\Rpmvx]wL:u_XJ'`z/zwItmcD)4JWRD|376$3qTZȡ}}~lB\Wn/>MR[F{JNtMcZ^&BE-!¢Եjmfdu\dgnAZQMWX]`TMQXPGFGBV"GzIJA.|O`n9&G)~†s;1kp~O=qR>PetOXXt`muL`9,A3l2&RTmQzNGis݇r=e@:F^{Z@XMlJ^ɀfT|sGmOjng^pyFZ^]qPkP.ڂna^O^_EVcnx`GkzV3[]쉀=j|}|wA8?ѳɥfPndm@]A-JlxbvsizoY.#"7S/Ξûeo]]naniKLaTSZDJaQDL^ZNJH]<<Ϥ~-5swdp_Gp25hJffcfk:\XywKmt]3FW:TOcHgzD^j[E``J\M*Rsetkr^rq{gts#akW=h=@HJQ{nMl,_etQ{inuS6ag[S0:MOE,V4W(nUU#uEEbL86-ΗjcXhdolGAXQOZJNOSDkgfcm:whTQHiHjPLb|laBcZ}u_UHxJ-^iN:9:]dSv7c}jj>/tnM{k{5q9rS2}~W+`>#~V4*51^@% TF4O[yhZ_`N6_967*˲^ldfonlH@OOINIE>UL8CMMTOOD@/n/v*+*$=.eVF a5Z[qu~i9PdU1;8CMYdU|iL#yEr_1uSti]PzMSFs?B3fcsppr5Uuڢ'9RtcZrO\VWʥp}[55AQ@]HqsE^2Fm,r,8=Ad=;dz(Ragr<^r&1&ReaaggbpgeIDGQFAJ97UR@FDE\CR@>C0O~|@$13*@/,blH^ hxuTf[yh6C^yZkK)\uK)aN]VxiwduOcv[svhkudȠPEoZJz%-@~gz_y_386pm|rlvxWiS>YOªDJMJMIQ'A\q[R;dLITfl8,RTO mSBSvyId_2:Z_0_4ML,H%SI@T<6QPFF76X@FO3*:>&$1b*4(#6&3\T-|>A%4Voo"2uCRaNe3=A334$3"kH.&//k2BV^Y@v1*vP>kzxt83nH]vl=D\cQHo(ObX^1R^MbNV^lQ_Z)STY`fWabEp"7L+5?^ZFCOfSGlQy`cc06BHo"<8<];C7Rv>8E>GrҪˋ]cgnopK>@GDCA743FDAF6'H1'(+6:. 0&4$/'4rYp^\gcS>1'جvusE3+]NSajXSFNT\N/p[]1y`h^fLyv|^[_Y^pwT;GSNOԴvw{b]n@gyWyyl?TIeUz='cpQoQC9]0sw\iz"JzP+r/D>r9SKA8,AgQ2XH6FE6LūѨ`baoqmMC@BDEA:31AGA<0$F6(&(02.3%%0ARMEIQjO_YYhuWs03L>|vB]Dcg=a_8f=XILSn(wciwuSYb\UON[uaujHNhws^shBfsPztueOXeMoynTI7]`6AVka+ZWVa:>][YDH[\`P[ggu{9ZZR$OX0.P5V#&V*"IC9d_YqqlKF=R^suDLtկ#OI<(3*cKg$-HF<7UlicE}`oO: + +&-=%'V1.@)XZPqkkEK:2DE642/$67/:'/1(,2HpRhOZWZ_YOXghibqef_bl\%/!0.-56;J&,Jh`XsF^DA~nUzqhbbjZVkhUnaaNYdHjI[NowsQyezumeyRt~/2!03[h(pls'"*4?*6T(S__MMgRncbaJP#s$B"$+#5#1?B%#5F1!BTɻS\TsiiB?94006<.-.364,%!%,0**?YbXRhV]Za[_``fpkaiiYeS*11-4(9'9l(dyM]s/]G_`PP~yhcpxljcZfvqXkUuxDZ%Q^vln}ymbpmiaqcjTxNkC2y|AfdMO**GT):]daZ\uLdx|;C.w7)HwB(!%()()+.K>3*!#:V:'7&&7#2q8WɷgfedggfeM9**241.90(&-8@D9ARee\[crUYV\bWRV\fonfbc_^\i]:2156+0%C;%A`t_hN\pHNQ\vyXjtrbOnJ@\vvn:u|srv|r{tGRWWXKadS}kamj]apjWg[gb`^AA5.8-+/ 0')/('(V0;^Q1`@/ONMJf]{kkj\[YPoikTWd~t~t}c[M_8{`}{ȵ{.N1FG4 ^$0e6f@YsITT !B#$%&(()*.'! (KU@B;).,Lv{rvh^`dc``ea^dlnh_B0'.1,)-=9657=ADBBHSYXTPPOQnjbjqd^jh^j`_cC<5-:+++-0"|!'F3HF=^I7<8JYqc4IswP:47fqqqwxgxs`usrz{vav9sgl~nkQ|ĮfONYWRlTTR4Vrtp!5PhW6|2~j;#"'))(,'"4FIC9%)HnXW\hpngab^]dlnh_VI=4.,08>:5239@E8?GNPNKGB[SΦڻqTnv\]aiddcG3.'2))'(6ـ&Ty$J[2o99DORbzmI^DCPKNRER{ɝ^q~}vhljY{F|qd}X{kxowQU+kNrTgeW'e^Zn2A[#I8"$$""#$" /BANbV9&Kp|qmhefgfdb`]`fnoia^eaM96=D345679;ec|aDT~pptw}}ƻ\GQ4W@O/;jf_z_uW3d1i4k&"+'%  &+%9 +9 .b^fmljigdefhlprokb]SH=6446=8@@:?9ZIKURMPT}c\soitH0&+'&!QT4A(#$+SiK65')/&IkllQQ=(CRm^GJ^ksjae{Gewńgpþ;TNXjq]lKw|O}M[0tsdVkXh+$--4PtnY$   + # T}[djkjieajjlprqnl_WJ?97778;:ANH\XPWVf»ǽpo^pykF+,5# $MQ89'$/YU..) +#mc^W+31'"*>-JRynF@DxvLr̢opo|¾NY=W*<:{ZzcWRV7\Vi<&* ~r|0   !#2GV]ejospknlmqqljmID=:<><:@8BA?>:HXDTXYQVΣŰʯa`l]RA0.83($! 6:06( ';VE((&'Fm^S7%$!.5-'gb0mwZq`_%CayxXv×NFQHWBUGb_5 %vv|wk'%-tx  +  %%4<[_ehmqpimikqpjjq@A@>::=@G9E@=@;PWKVT\QSƼ١K:"'% """+..3$%&GN2$!(.cWa6*& ,1ʄCtlV'9kh_Ir~Yior±^T\:PMT3CU{9& +>@_beglpoikgjqpjlu;BGB96>GM?=:888CJGBA>CRGVVN[V_ѽĻŻF&#(4/',#4=5*%"XeX8($ ''*&.'+0?vBWPWuf0~yN{Vo^k:x`Mcs®POMD4097zw" + 29   +  5?_hmhgnx}qnsyo`ds=>??>?BDFLGA@>DUSXWRXKZʿžϪT4$""" )73*)(9?2"!"N[O2!&(#%#+%, z\:WOlfzÖwNʱdX`ŕoEXcP\˺MYO|T$=Gj|o; !5:<1'  4=^hmjglrtnjo{ymjt==<;;>DHLPJEEAGZRPQSXI]¾¾O_eY:'%!"44*)19+,8SQXZoFGfah7ɹ-g[`?8L;||~>79A28570  + +39W`hgegjklfk{tllC@<97:@EKNJGG@DXQNRQVNZȻ`XLA8.%%.#$!)63)*+2$ 6ITB$#!$*.!  @QƾhÐm|?6NJJ@]YfZQGLv»ZjOHvWxyz\G.B22,   +  07QW^`_`dgjeivynecJF?98PTVYORQMlſúMPC:?FJQXHJEBB5),,/ #7F>4%%-('# *#/ IUrwvk}YvzCEAJY`SjRVk|shÿŬ].+WWK_`}uS'1M%%! +-5TVY[[\`dgfipndacFB<874>SVJRVPJE5*.$%"1?D,5>F< & !&C)NL;m{Zįj|z~tscu_ufKKhJMO9]Vsjƽ׾erN:,nUFaH"&D6!'$. +*5XWY\\[]daeimicciC?:668CY_WKMHDD7+..AF@&1CSG! ( 8",*H=$hWrrQlr~g:zoYy[reCL;wjZ2Hyoxf}'1 '@- 4B5& *5PPS\^\^c^ahlidehHE@;:@JPLKIML>@W\NVWSRRi_nſƱi>F?/-%  5FA3!!.GP5$2Y7 '&GR47E4X]rM:}^W\b|}OwfQJd|z}skcgwYUYPqre$]w]|=~|{Z'%-/r.)44/"#*)),,+/5AHHJTWYb^`lncbhfHGHA;EPIQPNMMJGClMS^εʽije54$7[SOXGAIOPQPB2% )?H;*#2&+bKrGe z*&<ʹwLEEV,e*{euwYnVhpsivh0N;d}zu4$Mj^im ji5K=54X:1?.11,$#'.349>=869HVSLOW]d`ajmdcgc?BFA=G?GNONLD;*(4E?%#/?.W|Aod@ 54pحĨx<|0,,&UplyúxU`}wik]WDyZns*9-*fks9P~t}\< >8=@DBBEGD=;=J`_PPX\bbailddf`9G`eWVYY\a_gkcbe_<=A>+#;?DIKJFA4?JF3 &%%!*hyAn2fz$--D4><487):SVOONMJIFEFB?@>;?FCYb]\YUX]]gj``d_@>@=><>OYPTRPOPONLͼh)+#,9lqo]^nG[!%1:AD9DB1$#" $%2"dL^d?2w{δ»׼¥4ŠðƧĻqaG)=D$nYo{`"N-E=uc9~c`Vi\j[dalgshc\TNMMNMKJMPLEBD>JT]`VNVUYeg\ZaaJ=M@:OPZPUNLOOPD³Жr)FHW`mtcEZfb+3 &5?4AC?%_+idsn|D3=)sSìۿӨ°ϴnzwV4011(&>]K}mRzR xSprhSxZ`kdlllnrvwxz{|wod^ZSNIFHHB?B8H^VU]YRZkb_f[[K@MC@TX^[WQUSJOXͶ̴RE"Aw6Bxbqu&gLRD6N PC;@#6,4,AyNd^NlR)6'0Š~ȶƸ̘jf&)25,%.>/L&| Y=hhP:~kghggjmpsvvxyyxtmhc`[TQOLE?5-AVSRYVRVd\\f[ZMAICCW]_WQPYYKWy׼ĵԼҼdNA7"`mG6P_\o+/-[pH503>:) +CFABYvc@1)ɸ݈{yԮϷѸũ~s'&/41.025&2S}4/9vE#uSZ_Ycjkorux{|{zy{zvpjokfca\TN>6@MLLPQSU`VVaWVLBCBCT]ZSPTVYPVĿűҾɺͺöe,ET6Nsc/)S;Xf`[][.OQB&1|#!%h{YGoyhM_VKO`irwyxy{}|{}}ytpppqomkiejediaSQHN\MENLUPTGLGH`WXX`NTNHȹƾŻq_B?$$\!]<(4`bZ[h7&@Bu?/]tttd64"JBcM-*7)cŶĶķµ~о-./0#2$ )@az}~sfXf[z}SɢcUV`mvwy{}yz|~~~|{{{yvqmkj^hccpj[\NOXJEOHKNTGLFEaXd\eVWEB¿δ°ɭƿT74"?;F?+!TCqE^J ?Q8cK-SjnkIzmT$6SR#%)03kơuֳ¾§~؅[38%#,-"%,Fnxyoe}wMθŠǷxTHVdnrtw{y|~~~}zuokiifma[kj_eYTWIGPB=WSUPNAjX[`TNZFmt5a_;1<)7"1<7I7pKVL*9@d.7dq^)3zt>y??H56"',QFʫ̾ҶՋ}ѱɸy0,'''%(/45=gn~sucϥҵmWNURTgsuyz|yuxy{|ukcbdZadbcfc\ZeVDFCD9&"%¥Ƹȉt{׳|(%%)*)*-'B7Idp`vıA3OMGKVV[kjt|}}~~ywwupkgdcZ_a_`eea[cRACAD@>/+8R_mvzz{z|yuqmjhgabeegjkjikaL<6?MWVUUSTWh^Y[BlstbZQQUX[^U`f]PPYdWDHQaԿiGnC%2@JB.S?Blr]b@Z[_Kj AaXRv?iDE/R1_*ٛľĸĿʭzhlr}źž˻Ѱq&"$*0012. )7P~+4=650*+,9Qjy{vr|}{tlfeghhlpqpomglkW<4I^ZZWXTYOj_]UIlqhdadihb\]ZVWWTLEH[Y[RD̽9p[EfuByo}wReIkhl@4ibpLQ8(# +"U_y*CjJ(EVԟʾ¸߹ϳ½ͿtooxӼü&261,+%4"+ƣɼ!#0'>=..1-)$'2G]lxw|}wmffhoorwwsoosuzd>4JaZ[WYU]Joa\Ni|iia_[Z[_aclbXQQTY\\ASSW:BM?u>js{{~TA<1icPhc@vacKz~S@e4cZil-RrYaa>q~JyOeI4qwzZ^[L*,.0/4B-_'uȷʯӷ.154.'(/.((Aovy2-4-/'4404=5))+*+,*)2?tuy|{wzustsmkopjvsY\{eQ\cPd^zcVgpmsifaZWWXZ[TQSQOQW`QZXMA#+f+Y_lY$39/#pmpz L^lW;n`{WD|+4g>pvcg6.47Ŭ̿ǿ̽ſ̫2.-5/'&+L~vlsoA?;L)1*(?44%%iدt%%6mĵ˼8)"#/ ÷cĽJ4=7-0A/!.97:>6(%+,*)(*/48018I_osp~wrommmsjqmX~{jcdWk]WOQzmnm_`^\ZZ\]`WVX]^\Z\\SUO? Za !G' O&sYaTE{ii?$5Y5T ++ʁ4DdRgXZ;&+41 C{̽ŸO*"#!+ľE)A9244-'/21*!#+*((*/36&3><38Ogipy~}wplc\b]PzSlVhb_ZKuejl^b_ZWWY\_TW[aa^ZXYUUI0 -   +\}h.(}ck_GBeo/luW~W EbXsjl@L=wCO)OʽûC%(  +! JƭǼK1=/*-2&*23'!&()+,.--,*66552026@QdpqkaXKK[VTtIrZ_\b\Kypeega`\WSSW]aY]adda]YRVC!!,   + ++igM_mQ#c>-E j 3nQCإLhU?J{LU|k\R·<*ATü'+'*-۾Ѩ#-0 '&/2' $('++.02200338><4//7>B?>BEDDLd\`vSxea[dZ]rpgaa`YWTTVZ_cW\__`_\W^D# #  ("'$ W !onqGjYcgsNB|0yF;"154&febuWU^ָ2=,CĹԽΞE&0#1˟M60*&,(02*&))$+**,/5<@1:GH9-3B8<8//7??HTqcjvaqihcdYroe^^[TVX[]^^^TY[Z[]\XY*+, ,k0* CbyiJeRKbs1ylsQCA;.%9&0L]da][׺R* 2/"&ɵ˽ɻϿ?I!&%)7,' 9Ĝ´ѼȬN6&4<:>AAA@@>;10//28=AP=3lBzZ5'8<){Rbdw"U&9D9'!:>a]ZĮO(+,$,*,­Χ,#7A]*!?08^ŲǼȿ4?@?<:8E>=A=34=8=;47?<3;Uq{quZNohgRyf^\YWWbaaa`^\\UfThfW[0+ "  !(BT]Qzy9dn[TWVgVlkBdfwpm4',!?);HRkZ6*)&)'$"""$$ŪՋ$GɭY',4%' qL+# ) +!Buo$!?ýɸF;KLGJLLLORPLD=9@QXF3BFDVB/A0=:8888658Nxtys~STyjhU~pa]]]]^`__`_^]]f_[aYX ) ' \o{toxr .ls_yx}fO?uUG,cP~Ww_l?s"/&2@JTXL:92,'$$',&$$$%&()µğM'$$;_5,2%&-!(")#k6. _%#   %zȤVJXUQOTUTSRMH=<42?PK>=N|xyv{^^`lnui^]___`^_^]]]]\[SgC% (  !yt|K?cbjbI8/19FTcwjyLHC'2L #5('*6GQQM@0*36-%&''')+++*Oİ)'*'()'+-=!&/!+- -M-n# a^ !!ĸLAJF@DNVXVUQMEE>37EG@IYvYDC1E@9118BHK?H}|zzvyahIxsyuja``][\_a`]]__\gi{l. + , +nQW,dlrm0 +=6?e{I +>OUZ) +YbInGBIGI*$+=;*!.=O_XI>==5-*,--/-,*(1ĺ\1'"!'.3$$,&%&%+,;J#" + PR% !&#!ȻɻSKPKF;IW[XVSPHGEA<=BH[jhC:@9MB?;877:HFBh}uaYaysmea`_\[]W^b_bih`]5ºy nu]sO"W?:{n AMS5!NX>$M*$/=7/('+0.*/9OXR9(')&"$(*++*)3PǺSa&0/($))"(+6\'.3/(#'%Cǵ¥IJq:<4762-2?MEBS|kXZkd^[\\\]`Zcfaae`U0) +(:@;bDͼ̤lwg5i[@I;[D"91-.644-*&# !%+->MK=1./&&(*+--,.Dþ˙EWHIG&+ 5%%"!''%$!%10U<83&+!)#06;>CILJUOJHIO_neQ?AEKN?JE;0(&',588~|_f[mx`[WXYY[acllc__UG]ï6 % {@[@Pv:0̿njqh8 8\~L45;\!"%/!+-!.)$.&*;$.#Uak?G81UNƽƵiG.F3!)/%'("&-"  /4+ #Z̽P38 ,)!,%+4-C??Qbcbgf_J>IWTG>DB8488?Rpi=UZ_tYW[]\^`ZUie^UPNML HƶP#27qr{`clĹ0=Gjqn6DgTR)By V!1%%3& 0@$StO\!1-.7%7%5%:*~M12L\=BB:,%'- #$$! +   "45 &2ʿI*8**A7+-B=')' -8,Ĭ!*/03+',326rmqnd``]eicIPX[TKJNILB:>=;APESB@[\ap͒ud[`mog[U[fkime\TMJHE&"Ƕ³Û¸VUW_.U3973!'pQ2$(87:?C@5,! $"&!$2R1% ##") +K\ SJ30# &"%NwdurpgfghppbLTZZSKIKFH=4:=:=FELDBJO`¿~nbsqNB1)1CORea\TMHGG+qýž³Ѿٱ۾BUGU4)%1&p343oS|xvL *@<77@G@+ !  !("00&"@N,"& !.#&,&6LŧNRC(=G(wtms{ijdeffkeQTX\ZTNJHFI;.047=8JA?HFO\Μrm|w.&,.+(,02RVXUPIEC9) »{ìŵDA>/B/=(32l2Fy pq_i+*3?IF7(""!#$!!"!&('.Ze[EH %.!!"("71URAXN?0:Hq_uyhpeb_djmrjSZ\]ZUPLGJM=*&+2<6X<9IBJMֳo3>#)-(!"'/6=CHIC:9 ð~ƫ33.yW21%&2<&3=~X'ur#%_91A?J-' "# <<#--!ZUD  "'%"*,,))-,(On_+<`UH91)xrompuumefkkjij`U\WTWUNLPEB=743584=CCABDDpz,"+.-,+$"25534=CDP"3¾Ƹּ@$>4*!J%.?`<<>?>935>FEDFGFk#3#&+7CC6*%)2851368@ ( PɴƷy\/<2;E7 5KTVZdY|j_e9el^<[39CI67/7?EEDFGFo7HLO:/'  !&,(5>:1--- (ιžʣu3-(A59+=@[cc.DtVgUQg4=?35>@3)% !"# (% $,")dV5" +  +))# &:/:90xr{qsnqrmhgmtqmhjomcYNNOPOKIH?:312443;?CAAEGGf)H6<94029>?BB>>CGHS|2W "!"%'-300353 +ȿσC:C%C**/ZqX`__i1kzFD&r75+404%!!!#$# " #(*,"&NiD   + #"It{eqhsolloqpnlklrodYRQQNHHIJI;@DB=Gn76C244WC)  !"!  -0"'('(0#( :{=!#$!! %Btszssjrnqtutrpnpopri[UVQUTMKLICA@?@AAA@28;;;:2)mc=0" +8>D!< ˻ƿƾаŰڢ~5=^6 3Pp_@EcH6S-;sqy*=<}k<:;85105:$% #%&$&"##(/.'$'&,,kK/&&$')#)6b|~srkiwlrusmjot{qsvdWUUPYTKLGCE=BGGC<9:2:=7.' ^eE:ze'46FT*7  ̬˿½ûƧK-3L%UqKCmrkE)fDS6$ +kd>*r $ )Q2<  ѺöɺĻĬ¸ΜX9˚r?,E02+=@0GYZWD2(^a)?.e."HCA,0147763080 %#%(!#.50*) &,-++-%0fN416,%k{ln{xmot}zvtvutsoqxq^XWPITRIOO@5047742259-$!hV<2 + J2:0 +ϺƼȷĹɴPѾ7.I*8(2&3[O_?:; T~uw||6&B2238;<72.*!&.!)89.$%'))))*,&/+ %xA"1).+)+-+'#,'(5E[fxnzjgu{~yutsropu|p]ZYMNWODKM?3(-230-++&, p`6)y Q2.5+UƵòúƶžE٣P/M!7+20:EU`_CF-,8xV]twg[E:88648;>=:;>CF9)%!#)%"#&'$!+7CIKG@9- .+1vL# (4AJKJ=F=RJQllprlpytvxwttuw}uqgXVWMEOPQM;+*1-% #&'!\X?-xz  %`K7;6<,ǼDzŗl˼*U%.A*(0/-8ƞQ4.$6DD*2&|rI%439U*D668667:<@@0!"$"&1&!*- #09CLPNH/ 7hd<$$#&4DQUPJMMKJRUZqluts~zsyxx{{zwu~ka\SRSNHLNOC(.%#"(/$ jHE%03  JyV'&ӶƿCî݇|XqjK@2*73$(%-5txr:8yS܂lt.$%y¶wôI>6=<353303774358E@+$!"$'$#&0CRWPUY0AhM'GDABEDB?xm0C9/:2()"$:(SWnE,[RyXcbJ*-tzŔȹ1A?(7B70/2:;879:N> &' #%###"%4FVSSYZK-TLQ)8]MIPVRLFGJAF,+;:6599669@DEF;$%$(#!!/D[f_VSSG2x>4K<\PRTTQMLJUNWK\Yfbw}z}|xx|SJGKMLPXNM8A<;P)"$%# #64$%6``F#-srejvv60UªĴ1YRVkfjx@ia9F! ,IK* )$-dxnC=X=E8' #'&$(-% ! + #>`zoeGZU@,aWY\\]]^`[aZhYd`lfp~ug`acgp{ziWFMQNLMOMV\LFF90&%#$+;B:00/,LߒcF';ؿemhicw}s[QJ$.27 +'^l;'*6 Cͮ\qt~p}}tQnj@C7-;K;.<0~lleo8dCHuwr˷ȾDzE198;@?:67L=@G4$'*$',<>6% + +9Tn}n.rWkY[^^[[\^[W]N\T_Yhwxj_YVZ[^dmri]HIJJIHHJOIDFKG7%&#"(6;2-/-&bB**]uxrevnqJJ +0 +K..?" +"-ιРKDqqta]9m"SV<(JqW3+(O8w6udRy]WewuɿžĹάȾ7:9PT6"%,**Fbm>d`9Ekheccb_]a]^I]_gWvZTY[bfhhjjdZWMILJDDGNCPTIJE9+%%%*/.*(*ARU1$XopqXw|uuMR:' + +'" R$7PյDOB}[Q|X?DyhkLV=uCM?SyYZC(f1b~^kUvst̼øū:88;AC@@B:B?H_S.$#)0,   ' lh1yVabk~pu^\WugWl_]\_gnnlf\S]OIPRJJSU=KPEMP@0)  !%14&$2xaD;)?jowodg~vlve;2?!   .#   +Tj>Mwy'EgUJ;8biOp]sQ`ZU]1"mf +. )IhJ6ouMGeŹſ÷@@;8===>?CGIOJPZ3&%*(23)#)9L]C4\d~x^b`ĝ{Viih`dikg^QGXdjbTKJJJLQRNFCED)"3"$  "#',pa8'/Lzr{sjj`ssmaH3-&    &PA>+ (ciN%[xn}v]}v=9uBx":a[V~Kz_'>IK2->f0<>DJGWHZMT+#).351'$+/&%HJ$SAeYX@wlG_}}?"i;u)Ph?Tr»k?>GIBDDCABFHD^QT$6#)CA>94/,)%!#$%%L4SD3PBl|Xpizogǔ|u~ogjqporwrg__ca\^isuttspl\E#    K87BZ4}̣x^^nyimqleswhWC (IP]tcr+$mJ,D'g4pQ/-#CyKw~|pkLS}+^QdZq|IqcO%~pUS$f;K=8I>D@DHHFEIMLUVQ "(.!#%%5@;]=0oft\^UUäh\phygrijilpxwtoptrmmou{wrrvvgK  + +6D;KY1zgabzrhjkjZkvfH20  +BKCDjtpO1+BtZC/('/*rqwm_z[qB F7XhaL~E_xc37A.CpELP\W;990 =KEGFHX]aX7 -CAOjj:j[)4!p͜rdnT(EHc[no`z4/b<",MF>6F>M~«?A?;53699<@BBDHLDLSfV$;8JVOC><6 bOHG7mkJst|yqttnntusrv||xssvxwuuy}}ud6$ @5BQ:#Kgiulng_sxrg;<>- +-2DF?ABA@;60#!1?K:W4,BD<($-Ia|(Way;TGUCMpyw#E&VFUGJROwjz?<<:9::9969<>>AEIO]WK9!G,=NWK:- !(e:+Xie|Ժ|wqortsrtwtrtxyxvxxwustwz}zrH,(  eYA:%#zolerxjhkp[RLC=./@8:AC?CE:2)!!(.G:L^L4@5H=8!(.Jh#+cocjGZ1c|wml`oro"':dLr\Oc?lqi]kudO=847<@<7579<=@DGNYTG9&PF9>9*O[\Z6=EFaǨb͝hybfjlqwwtvusrtvxz|yusstuuxywR62#"'"! }m5("$?rufdiTZ[m`bG:>HB9503:45>C<567:7-2;VJ<9@C?8779<@DGHOMRQQ()) +  d4\@oM^±~ưWRWSY`ǒbXW_a_covsrtwxxwxytssuwvplnomH02.4153+)-(sh55##$?O6;@8=4?D\RFABC@<668;@EGGITL$*%(ĩ:yI=6ONdSbZLRK͢^_]__]_hpsmorvxyxwqrtvwtolnsvQ86.1.43,)+)!vn@8)#*K6LQKILE4>=;5A6.-+,0-&3FL8J8@DFK\9 # +Ǵjc]x96[`F2Kfe}Q,6M3;BHJV=<999=CHJICRE *( +  .3NK*\{~jMSOWSPiX`ta_^``bjrkjjntwwusttspnmmntxU;8/3120,)))*ejH?4.-[YEKMDDHXG=8%8,2534840>HA3K:EHGJd. + ɳHkOni>MKTZQP@Y-4[l;L>-lFVA]|T\L~7$9.)+s1EB*31:,?9?0?ACA<878><<@GLNMLSH.   \]ë6V;&fSEwuWFW`Qr`Yc^^dc_eqlgehovwuxxvqkiknmqsM562820,($#',SG^F507DJD=AML=IH"/¼s([0sgdRRF~yKD;*X&fsij"->?DINOOYLV,<.  +! +g !\DI%nM3y{~j\RKMQRQQR[ef\Wcvtgemnho|zwtqonnmqY>7964A9/)())(,!DMTX?POOONMKIIAC>73786CCCIPPKGKCFGL?/Mj(Dƴűfz%oY@}^^cFK>>J6&4EmT/tyNbZJb+ oB,;}=)VG@KAI@49<73989>?<;>A@@BFKOQ*;.2(! ('  %EE;K'=PqYQKGGEFJRSZdhc]ajuhbfhgkrpooonlihddQ<:;88;9510.,*2QdcFRRQONLJHHIJF=779;MT[^[WTSTNI9Ig{̶ŵy6)&K;3ZTFPmT^X[G.&!x9!-. 9Ep8A?<+FHMA,=4327953535:;9:=BDGIJKNO:#",1-253%'  +&-OQWh^j]}ΚHHHFAAN[UYbkle^]xpfehllihjmpomif]WD9<;7;148:74109O^`OYTROLKJKKQQLA75;AX^egb]ZYFN`eʹ֐)&(`D""SjH oJ[}rb"&% 6" JG_],,5/;6C?5/6+2I.863<<763127:89=BGMQQPPP\) $9>!"% gbT/0JIY{VTZCDGGDEQ_U[dmmg]Xrsohjomdbeimnmji^Q=4841805983+'&,JSW^WYTRNIIJLOQOKC<;CK]_`bb`]\_Ro¤μԮ*SQ%D)+ K ]#WgQ% Dj;D-# 7+ 1Zs>0@.<%C>006859,544>?89646;=<>BEGKORUY[(+eln&B *% + Gcθ5;<+MTr^J=Y>CABDCCGNS\hlic^[anslhkib`abdeeeddU<01-(/-00,"4S\^`\UURLHFGJLIGGHGHNT\ZXZ]`_]sg»Ѿѹ5$8r_+ +% F3N*@bdղ|h0.YH>8IG%D"FD'* 8?#Y7#@BA03<;*3.+7;417716745:=>?:;ERUNHFPTX[^`cehMh´,Zug !AEFGNVXL3 #Woy, "').1 D*:==DC?LN6B@@ACELUMZdgjke_a`bfhebb_befc]VRHQD003.*  \ No newline at end of file diff --git a/images/lagrange.jpg b/images/lagrange.jpg new file mode 100644 index 0000000..d21e849 Binary files /dev/null and b/images/lagrange.jpg differ diff --git a/images/projecteur1.jpg b/images/projecteur1.jpg new file mode 100644 index 0000000..4fe114d Binary files /dev/null and b/images/projecteur1.jpg differ diff --git a/images/svd1.jpeg b/images/svd1.jpeg new file mode 100644 index 0000000..8c7c90d Binary files /dev/null and b/images/svd1.jpeg differ diff --git a/intro_to_ml.Rproj b/intro_to_ml.Rproj new file mode 100644 index 0000000..a680331 --- /dev/null +++ b/intro_to_ml.Rproj @@ -0,0 +1,15 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX + +SpellingDictionary: fr_FR