path: root/PVCM/cama/fr/ma60 Méthode du gradient conjugué.ipynb

{
 "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": []
  }
 ]
}