init: initial commit

1 files changed, 365 insertions, 0 deletions

diff --git a/PVCM/cama/fr/ma60 Méthode du gradient conjugué.ipynb b/PVCM/cama/fr/ma60 Méthode du gradient conjugué.ipynb
new file mode 100644
index 0000000..7b01941
--- /dev/null
+++ b/PVCM/cama/fr/ma60 Méthode du gradient conjugué.ipynb
@@ -0,0 +1,365 @@
+{
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4,
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "Cette leçon a 3 effets positifs :\n",
+    "\n",
+    "* manipuler des bases de sous-espaces\n",
+    "* mélanger les dérivées partielles et les vecteurs\n",
+    "* trouver un algorithme plus performant que la méthode du gradiant pour résoudre \n",
+    "  $A \\, {\\bf x} = {\\bf b}$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "# Méthode du gradient conjugué\n",
+    "\n",
+    "On a vu que dans la méthode du gradiant on peut calculer le µ optimal pour\n",
+    "aller directement au minimum de la fonctionnelle $J$ dans le long de la ligne\n",
+    "${\\bf x}^k \\; + µ \\, \\nabla J({\\bf x}^k)$ avec µ qui varie de $- \\infty$\n",
+    "à $+ \\infty$ (cf ma33).\n",
+    "\n",
+    "On a aussi vu que si on a le µ optimal, alors la plus forte pente suivante, \n",
+    "$\\nabla J ({\\bf x}^{k+1})$, sera orthogonale à la pente qui définie la droite sur laquelle on cherche µ. On a donc\n",
+    "\n",
+    "$$\\nabla J ({\\bf x}^{k+1})^T \\, . \\, \\nabla J ({\\bf x}^k) = 0$$\n",
+    "\n",
+    "C'est bien car cela veut dire que le minimum suivant, ${\\bf x}^{k+2}$, sera le minimum de l'espace généré par $\\nabla J ({\\bf x}^{k+1})$ et $\\nabla J ({\\bf x}^k)$\n",
+    "(ce qui est mieux que d'être seulement le minimum le long d'une ligne). \n",
+    "\n",
+    "Malheureusement rien me dit que ${\\bf x}^{k+3}$ ne sera pas calulé le long de notre\n",
+    "première direction, $\\nabla J({\\bf x}^k)$. On a vu dans le cas\n",
+    "du gradiant sans recherche du µ optimal qu'on peut osciller entre 2 valeurs. Cela\n",
+    "n'est pas possible avec le µ optimal mais on peut tourner en rond sur quelques valeurs.\n",
+    "\n",
+    "Même si on ne tourne pas en rond et qu'on converge est il probable qu'on cherche\n",
+    "à nouveau un minimum dans une partie de l'espace $ℝ^n$ qu'on a déjà explorée.\n",
+    "\n",
+    "Une recherche optimal du minimum d'une fonction convexe dans un espace $ℝ^n$ ne\n",
+    "devrait pas prendre plus de $n$ itération si on est toujours capable de caluler\n",
+    "le µ optimal dans la direction choisie. Pour s'en convaincre il suffit de choisir\n",
+    "comme directions les vecteurs de la base de notre espace $ℝ^n$. Une fois qu'on\n",
+    "a trouvé le minimum dans toutes les directions de la base, on a le minimum global.\n",
+    "Il n'y a plus aucune direction possible de recherche qui ne serait pas la combinaison linéaire des vecteurs\n",
+    "déjà utilisés."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "#### Cas 2D\n",
+    "<table><tr>\n",
+    "    <th>\n",
+    "        \n",
+    "Si on cherche le point ${\\bf x} = (x_1, x_2)$ qui minimise $J({\\bf x})$ (convexe)\n",
+    "et qu'on sait calculer le µ optimal de notre méthode de gradient,\n",
+    "alors on trouve le minimum globale en 2 itérations au maximum.\n",
+    "\n",
+    "Mais on \n",
+    "prend un µ trop grand qui nous fasse passer au dessus du minimum\n",
+    "global, alors à la prochaine itération la plus grande pente sera exactement\n",
+    "l'opposée de celle qu'on vient d'utiliser et on oscillera.\n",
+    "\n",
+    "Sur le dessin on voit bien qu'on trouve le minimum global en 2 itérations\n",
+    "mais on se dit que c'est un coup de chance et d'ailleurs est-ce vraiment\n",
+    "en 2 itérations puisqu'il est écrit ${\\bf x}^k$ ? C'est toute la limite\n",
+    "d'un dessin qui veut expliquer ce qui se passe dans $ℝ^n$ à des humains\n",
+    "habitués à $ℝ^3$ (cf cette vidéo sur <a href=\"https://www.youtube.com/watch?v=LQFkUjYzOn8\">la 4e dimension</a> en dehors\n",
+    "du cours).\n",
+    "\n",
+    "</th><th>\n",
+    "<img src=\"images/gradient.png\" width = \"400px\"/>\n",
+    "</th></tr></table>\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "## Générer une base de $ℝ^n$\n",
+    "\n",
+    "Si on veut trouver notre minimum global en $n$ itération au maximum, il faut que\n",
+    "nos directions ne soit pas redondante et donc que les $n$ premières directions\n",
+    "génèrent $ℝ^n$ ou forment une base de $ℝ^n$.\n",
+    "\n",
+    "Il faut donc que la nouvelle direction ${\\bf d}^k$ soit orthogonales à toutes les directions\n",
+    "précédentes ou qu'elle permette de trouver une base qui génère un\n",
+    "espace de dimension $k+1$ (on commence avec ${\\bf d}^0$).\n",
+    "\n",
+    "Est-ce possible ?\n",
+    "\n",
+    "### Le cas $A {\\bf x} = {\\bf b}$\n",
+    "\n",
+    "Dans ce cas la fonctionnelle à minimiser est\n",
+    "\n",
+    "$$\n",
+    "J({\\bf x}) = \\frac{1}{2} \\, {\\bf x}^T \\, A \\, {\\bf x} - {\\bf b}\\, . {\\bf x}\n",
+    "$$\n",
+    "\n",
+    "et son gradient est $\\nabla J({\\bf x}) = A \\; {\\bf x} - {\\bf b}$ si A est \n",
+    "symétrique.\n",
+    "\n",
+    "\n",
+    "On sait que si on calcule ${\\bf x}^k$ comme avant on a seulement d'orthogonalité\n",
+    "de deux directions successives. Cela arrive lorsque\n",
+    "µ est choisi pour minimiser $J$ le long de la droite de direction \n",
+    " $\\nabla J({\\bf x}^k)$.\n",
+    " \n",
+    "Que se passe-t-il si ${\\bf x}^{k+1}$ **minimise $J$ dans l'espace \n",
+    "${\\bf G_k}$ généré par\n",
+    "toutes les directions précédentes** ? \n",
+    "\n",
+    "$$\n",
+    "J({\\bf x}^{k+1}) = \\min_{\\bf v \\in G_k} J({\\bf x}^k + {\\bf v})\n",
+    "$$\n",
+    "avec ${\\bf G_k} = span\\{{\\bf d}^0, {\\bf d}^1,\\dots, {\\bf d}^k\\} \n",
+    "                =  \\left\\{ {\\bf v} = \\sum_{i=0}^{k} \\alpha_i \\, {\\bf d}^i \\quad \\forall \\alpha_i \\in ℝ \\right\\}$.\n",
+    "\n",
+    "Alors toutes les dérivées partielles par rapport aux vecteurs de ${\\bf G_k}$ \n",
+    "sont nulles (tout comme $\\partial J / \\partial x$ et $\\partial J / \\partial y$ \n",
+    "sont nulles au minumum de J dans le dessin ci-dessus). Vectoriellement\n",
+    "cela s'écrit ainsi :\n",
+    "\n",
+    "$$\n",
+    "\\nabla J({\\bf x}^{k+1}) \\, . \\, {\\bf w} = 0 \\quad \\forall {\\bf w} \\in {\\bf G_k} \\qquad (0)\n",
+    "$$\n",
+    "\n",
+    "ce qui se vérifie facilement si ${\\bf w}$ est un des vecteurs de la base :\n",
+    "\n",
+    "$$\n",
+    "\\nabla J({\\bf x}^{k+1}) \\, . \\, {\\bf e}_i = \n",
+    "\\begin{bmatrix}\n",
+    "\\partial J / \\partial x_{1} \\\\\n",
+    "\\partial J / \\partial x_{2} \\\\\n",
+    "\\vdots \\\\\n",
+    "\\partial J / \\partial x_{i} \\\\\n",
+    "\\vdots \\\\\n",
+    "\\partial J / \\partial x_{n} \\\\\n",
+    "\\end{bmatrix}\n",
+    "\\, . \\,\n",
+    "\\begin{bmatrix}\n",
+    "0 \\\\\n",
+    "0 \\\\\n",
+    "\\vdots \\\\\n",
+    "1 \\\\\n",
+    "\\vdots \\\\\n",
+    "0 \\\\\n",
+    "\\end{bmatrix}\n",
+    "=\n",
+    "\\frac{\\partial J}{\\partial x_i}({\\bf x}^{k+1}) \n",
+    "$$\n",
+    "\n",
+    "On note que dire la dérivée partielle de $J$ dans une direction ${\\bf w}$\n",
+    "de ${\\bf G_k}$ est nulle,\n",
+    "revient à dire que $\\nabla J({\\bf x}^{k+1})$ est orthogonal à ${\\bf w}$ \n",
+    "(ou nul si $\\nabla J({\\bf x}^{k+1})$ la même dimension que  ${\\bf G_k}$).\n",
+    "\n",
+    "#### Générons les directions ${\\bf d}^i$\n",
+    "\n",
+    "Ainsi notre formule itérative devient :\n",
+    "\n",
+    "$$\n",
+    "{\\bf x}^{k+1} =  {\\bf x}^k + µ^k\\, {\\bf d}^k\n",
+    "$$\n",
+    "\n",
+    "Pour calculer les ${\\bf d}^k$ on utilise la formule qui nous dit que les\n",
+    "dérivées partielles de $J$ par rapport à un vecteur ${\\bf w \\in G_k}$ sont\n",
+    "nulles. Et comme les directions ${\\bf d}^i$ génèrent l'espace ${\\bf G_k}$\n",
+    "il suffit que les dérivées partielles de $J$ par rapport à tous les ${\\bf d}^i$\n",
+    "soient nulles :\n",
+    "\n",
+    "$$\n",
+    "\\nabla J({\\bf x}^{k+1}) \\, . \\, {\\bf d}^i = 0 \\quad \\forall i \\le k \\qquad (1)\n",
+    "$$\n",
+    "\n",
+    "On déroule les calculs\n",
+    "\n",
+    "$$\n",
+    "\\begin{align}\n",
+    "(A\\, {\\bf x}^{k+1} - b) \\, . \\, {\\bf d}^i &= 0 \\quad \\forall i \\le k \\\\\n",
+    "(A\\, ({\\bf x}^{k} + µ^k \\, {\\bf d}^k) - b) \\, . \\, {\\bf d}^i &= 0 \\quad \\forall i \\le k \\\\\n",
+    "(A\\, {\\bf x}^{k} - b) \\, . \\, {\\bf d}^i + µ^k \\, A \\, {\\bf d}^k \\, . \\, {\\bf d}^i &= 0 \\quad \\forall i \\le k \\\\\n",
+    "\\nabla J({\\bf x}^{k}) \\, . \\, {\\bf d}^i + µ^k \\, A \\, {\\bf d}^k \\, . \\, {\\bf d}^i &= 0 \\quad \\forall i \\le k \\\\\n",
+    "\\end{align}\n",
+    "$$\n",
+    "\n",
+    "Si $i < k$ alors le premier terme est nul ce qui donne\n",
+    "\n",
+    "$$\n",
+    "A \\, {\\bf d}^k \\, . \\, {\\bf d}^i = 0 \\quad \\forall i < k \\qquad (2)\n",
+    "$$\n",
+    "\n",
+    "Comme A est symétrique, on a aussi $A \\, {\\bf d}^i \\, . \\, {\\bf d}^k = 0 \\quad \\forall i < k\\quad$ ($\\sum_{m=1}^N \\sum_{n=1}^N a_{m,n}\\, d^i_n\\, d^k_m = \\sum_{n=1}^N \\sum_{m=1}^N a_{n,m}\\, d^k_m\\, d^i_n$). \n",
+    "\n",
+    "Ainsi on a des conditions pour construire la nouvelle direction ${\\bf d}^k$.\n",
+    "\n",
+    "On suppose que la nouvelle direction est le gradiant + une correction qui est une combinaison linéaire des directions précédantes :\n",
+    "$$\n",
+    "{\\bf d}^k = - \\nabla J({\\bf x}^k) + \\sum_{j=0}^{k-1} β_j {\\bf d}^j\n",
+    "$$\n",
+    "alors l'équation (2) devient\n",
+    "$$\n",
+    "A \\, \\left( - \\nabla J({\\bf x}^k) + \\sum_{j=0}^{k-1} β_j {\\bf d}^j \\right) . {\\bf d}^i = 0 \\quad \\forall i < k\n",
+    "$$\n",
+    "La somme va faire apparaître des $A\\, {\\bf d}^j {\\bf d}^i$ qui sont tous nul d'après (2) et sa version symétrique sauf $A\\, {\\bf d}^i {\\bf d}^i$. Donc\n",
+    "\n",
+    "$$\n",
+    "β_i = \\frac{A\\, \\nabla J({\\bf x}^k) \\, . \\, {\\bf d}^i}{A\\, {\\bf d}^i \\, . \\, {\\bf d}^i} \n",
+    "$$\n",
+    "\n",
+    "et ainsi\n",
+    "$$\n",
+    "{\\bf d}^k = - \\nabla J({\\bf x}^k) \n",
+    "            + \\sum_{i=0}^{k-1} \\frac{A\\, \\nabla J({\\bf x}^k) \\, . \\, {\\bf d}^i}\n",
+    "                                    {A\\, {\\bf d}^i \\, . \\, {\\bf d}^i} \\; {\\bf d}^i\n",
+    "$$                                   \n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "## 2e tentative\n",
+    "\n",
+    "La méthode trouvée est bien car elle montre qu'on peut optimiser la descente de gradiant mais le calcul\n",
+    "de ${\\bf d}^k$ est bien trop lourd. Aussi on va tout recommencer pour trouver une autre facon plus efficace.\n",
+    "\n",
+    "### Travaillons dans la base des $\\nabla J({\\bf x}^i)$\n",
+    "\n",
+    "On repart de l'équation (0) qui permet de \n",
+    "montrer que les $\\nabla J({\\bf x}^i), i \\le k$ forment une bases du sous espace $\\bf G_k$. Aussi tout vecteur s'exprime dans cette base et donc\n",
+    "\n",
+    "$$\n",
+    "{\\bf x}^{k+1} = {\\bf x}^k + μ^k \\, {\\bf d}^k \\quad \\textrm{avec} \\quad {\\bf d}^k = - \\nabla J({\\bf x}^k) + \\sum_{i=0}^{k-1} γ^i \\, \\nabla J({\\bf x}^i)\n",
+    "$$\n",
+    "\n",
+    "On note que ces ${\\bf d}^i, i \\le k$ forment aussi une base de $\\bf G_k$.\n",
+    "\n",
+    "L'équation (2) nous dit que $A \\, {\\bf d}^k \\, . \\, {\\bf d}^i = 0 \\quad \\forall i < k$, aussi\n",
+    "\n",
+    "$$\n",
+    "\\begin{align}\n",
+    "< A \\, {\\bf d}^k, {\\bf x}^{j+1} - {\\bf x}^j > &= 0 \\quad \\forall j < k \\\\ \n",
+    "< {\\bf d}^k, A \\, ({\\bf x}^{j+1} - {\\bf x}^j) > &= 0 \\quad \\textrm{car A est symétrique} \\\\\n",
+    "< {\\bf d}^k, A \\, {\\bf x}^{j+1} - {\\bf b} - A {\\bf x}^j + {\\bf b} > &= 0 \\\\\n",
+    "< {\\bf d}^k,  \\nabla J({\\bf x}^{j+1}) - \\nabla J({\\bf x}^j) > &= 0 \\\\\n",
+    "< - \\nabla J({\\bf x}^k) + \\sum_{i=0}^{k-1} γ^i \\, \\nabla J({\\bf x}^i) ,  \\nabla J({\\bf x}^{j+1}) - \\nabla J({\\bf x}^j) > &= 0 \\quad (3) \\\\\n",
+    "\\end{align}\n",
+    "$$\n",
+    "\n",
+    "Si j = k-1 alors\n",
+    "$$\n",
+    "\\begin{align}\n",
+    "< - \\nabla J({\\bf x}^k), \\nabla J({\\bf x}^k) > + \\, γ^{k-1} < \\nabla J({\\bf x}^{k-1}), -\\nabla J({\\bf x}^{k-1}) > &= 0 \\\\\n",
+    "γ^{k-1} = -\\frac{< \\nabla J({\\bf x}^{k}), \\nabla J({\\bf x}^{k}) >} \n",
+    "             {< \\nabla J({\\bf x}^{k-1}), \\nabla J({\\bf x}^{k-1}) >}\n",
+    "\\end{align}\n",
+    "$$\n",
+    "\n",
+    "Si j = k-2 alors\n",
+    "$$\n",
+    "\\begin{align}\n",
+    "γ^{k-1} \\, < \\nabla J({\\bf x}^{k-1}), \\nabla J({\\bf x}^{k-1}) > + \\, γ^{k-2} < \\nabla J({\\bf x}^{k-2}), -\\nabla J({\\bf x}^{k-2}) > &= 0 \\\\\n",
+    "γ^{k-2} = γ^{k-1} \\frac{< \\nabla J({\\bf x}^{k-1}), \\nabla J({\\bf x}^{k-1}) >} \n",
+    "             {< \\nabla J({\\bf x}^{k-2}), \\nabla J({\\bf x}^{k-2}) >}\n",
+    "        = \\frac{< \\nabla J({\\bf x}^{k}), \\nabla J({\\bf x}^{k}) >} \n",
+    "             {< \\nabla J({\\bf x}^{k-2}), \\nabla J({\\bf x}^{k-2}) >}\n",
+    "\\end{align}\n",
+    "$$\n",
+    "et en continuant les calculs, on a $\\forall j < k-1$\n",
+    "$$\n",
+    "    γ^j = \\frac{< \\nabla J({\\bf x}^{k}), \\nabla J({\\bf x}^{k}) >} \n",
+    "             {< \\nabla J({\\bf x}^{j}), \\nabla J({\\bf x}^{j}) >}\n",
+    "$$\n",
+    "\n",
+    "La formule générale de ${\\bf d}^k$ est donc\n",
+    "$$\n",
+    "\\begin{align}\n",
+    "{\\bf d}^k &= -\\nabla J({\\bf x}^k) \n",
+    "             - \\frac{< \\nabla J({\\bf x}^{k}), \\nabla J({\\bf x}^{k}) >} \n",
+    "                    {< \\nabla J({\\bf x}^{k-1}), \\nabla J({\\bf x}^{k-1}) >} \\, \\nabla J({\\bf x}^{k-1})\n",
+    "             + \\sum_{i=0}^{k-2}  \\frac{< \\nabla J({\\bf x}^{k}), \\nabla J({\\bf x}^{k}) >} \n",
+    "             {< \\nabla J({\\bf x}^{i}), \\nabla J({\\bf x}^{i}) >} \\, \\nabla J({\\bf x}^i) \\\\\n",
+    "         &= -\\nabla J({\\bf x}^k) + β \\left( -\\nabla J({\\bf x}^{k-1}) \n",
+    "                                            + \\sum_{i=0}^{k-2}  \\frac{< \\nabla J({\\bf x}^{k-1}), \\nabla J({\\bf x}^{k-1}) >} \n",
+    "             {< \\nabla J({\\bf x}^{i}), \\nabla J({\\bf x}^{i}) >} \\, \\nabla J({\\bf x}^i) \\right) \\\\\n",
+    "         &= -\\nabla J({\\bf x}^k) + β \\, {\\bf d}^{k-1} \\quad \\textrm{avec} \\quad β = \\frac{< \\nabla J({\\bf x}^{k}), \\nabla J({\\bf x}^{k}) >} \n",
+    "                    {< \\nabla J({\\bf x}^{k-1}), \\nabla J({\\bf x}^{k-1}) >} \n",
+    "\\end{align}\n",
+    "$$\n",
+    "\n",
+    "On a donc réussi à exprimer la nouvelle direction seulement en fonction du gradient courant et de la direction précédente."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "lang": "fr"
+   },
+   "source": [
+    "### Nouveau calcul de μ\n",
+    "\n",
+    "Le calcul précédent de μ reste valable puisque ce calcul ne dépend pas de la facon d'exprimer ${\\bf d}^k$.\n",
+    "\n",
+    "$$\n",
+    "µ^k = -\\frac{\\nabla J({\\bf x}^{k}) \\, . \\, {\\bf d}^k}\n",
+    "            {A \\, {\\bf d}^k \\, . \\, {\\bf d}^k}\n",
+    "\\,         \n",
+    "$$\n",
+    "\n",
+    "Avec notre choix de ${\\bf d}^k$ cela donne\n",
+    "\n",
+    "$$\n",
+    "µ^k = -\\frac{<\\nabla J({\\bf x}^{k}), - \\nabla J({\\bf x}^{k}) + β \\, {\\bf d}^{k-1}>}\n",
+    "            {<A \\, {\\bf d}^k,  {\\bf d}^k>}\n",
+    "    = \\frac{<\\nabla J({\\bf x}^{k}), \\nabla J({\\bf x}^{k})>}\n",
+    "            {<A \\, {\\bf d}^k,  {\\bf d}^k>}\n",
+    "$$\n",
+    "puisque ${\\bf d}^{k-1} \\perp \\nabla J({\\bf x}^{k})$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ]
+}
\ No newline at end of file