diff options
Diffstat (limited to 'PVCM/cama/en/ma62 Gradient conjugué -- Exercice.ipynb')
| -rw-r--r-- | PVCM/cama/en/ma62 Gradient conjugué -- Exercice.ipynb | 409 |
1 files changed, 409 insertions, 0 deletions
diff --git a/PVCM/cama/en/ma62 Gradient conjugué -- Exercice.ipynb b/PVCM/cama/en/ma62 Gradient conjugué -- Exercice.ipynb new file mode 100644 index 0000000..0bca368 --- /dev/null +++ b/PVCM/cama/en/ma62 Gradient conjugué -- Exercice.ipynb @@ -0,0 +1,409 @@ +{ + "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": {}, + "source": [ + "## Programmer le gradient conjugué\n", + "\n", + "A partir de ce [cours sur le gradient conjugué](http://perso.unifr.ch/ales.janka/numeroptim/07_conjgrad.pdf) programmez en Python + Numpy le gradient conjugué en exploitant les astuces mathématiques indiquées pour optimiser\n", + "votre code.\n", + "\n", + "* Effectuez des tests pour valider votre code. \n", + "* Comparez la vitesse de convergence à celle du gradient avec μ optimal. Tracez des courbes de convergence (cf la feuille qui en parle)\n", + "* Comparez les temps de calcul.\n", + "\n", + "\n", + "Note : Veuillez écrire des fonctions les plus propres possibles, en particulier qui n'utilisent pas des variables globales comme c'est le cas dans ma correction du gradient (ma33)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import numpy.linalg as lin\n", + "import matplotlib.pylab as plt\n", + "\n", + "%matplotlib inline\n", + "%config InlineBackend.figure_format = 'retina'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def make_system(N):\n", + " A = 1.0 * np.random.randint(-10, 10, size=(N,N))\n", + " A[np.diag_indices(N)] = 0.1 + np.abs(A).sum(axis=0) # diag dominante\n", + " A = A + A.T # symétrique\n", + " A = A / np.abs(A).sum(axis=0).mean()\n", + " b = 1.0 * np.random.randint(-10,10,size=(N))\n", + " return A, b\n", + "\n", + "A, b = make_system(2)\n", + "print(A, \"\\n\\n\", b)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def gradient_conjugué(A, x0, b, error=1E-9, convergence=False):\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# vérfions que ca marche\n", + "\n", + "# combien faut-il d'itérations\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def compute_error(N, method=gradient_conjugué):\n", + " A, b = make_system(N)\n", + " x = method(A, np.zeros(N), b)\n", + " err = A @ x - b\n", + " return np.sqrt(err @ err)\n", + "\n", + "compute_error(10)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Comparons avec le gradient simple" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def gradient(A, x0, b, e = 1E-9, convergence=False, max_iterations=1000):\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# vérifions que ca marche\n", + "\n", + "compute_error(10, method=gradient)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Perfs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# comparons les performances\n", + "\n", + "seed = 123\n", + "np.random.seed(seed)\n", + "\n", + "%timeit compute_error(1000, method=gradient)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "seed = 123\n", + "np.random.seed(seed)\n", + "\n", + "%timeit compute_error(1000, method=gradient_conjugué)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Nombre d'iteration dans les 2 cas\n", + "\n", + "On teste avec N=1000 puis N=10000." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "N = 1000\n", + "A,b = make_system(N)\n", + "x0 = np.zeros(N)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Pour le gradient simple" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "err = gradient(A, x0, b, convergence=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#plot" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Pour le gradient conjugué" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "err = gradient_conjugué(A, x0, b, convergence=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#plot" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Un cas réel\n", + "\n", + "Logiquement vous devriez être décu aussi on va tester avec un problème réel qui correspond à cet exemple\n", + "de l'[équation de Poisson](https://doc.freefem.org/tutorials/poisson.html). Le système matriciel de ce problème est téléchargeable \n", + "[ici](https://www.lrde.epita.fr/~ricou/cama/data/Ab.npz). Une fois le fichier sauvé, pour récupérer\n", + "A et b faites :\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "npz = np.load('/tmp/Ab.npz',allow_pickle=True)\n", + "A = npz['A']\n", + "b = npz['b']" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "* Faites une étude rapide de A, indiquez quel pourcentage des valeurs de A est différent de 0. Afficher l'image de la matrice avec `plt.imshow(A)` (faire une grande image pour voir quelque chose).\n", + "* Refaites la comparaison entre les deux méthodes avec ce système matriciel.\n", + "* Regardez la documentation de `lin.solve` (en particulier les options) et comparer `lin.solve` à vos deux algorithmes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Comparaison gradient simple et conjugué" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%time err = gradient_conjugué(A, np.zeros(len(A)), b, convergence=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#plot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# le gradient simple\n", + "\n", + "%time err = gradient(A, np.zeros(len(A)), b, convergence=True, max_iterations=10000)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "On voit la supériorité du gradient conjugué tant en nombre d'itérations (175 contre 7800) qu'en temps de calcul (environ 40 fois plus rapide)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Comparaison avec `lin.solve` de Scipy\n", + "\n", + "Le solveur de Scipy a des option que celui de Numpy n'a pas." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# regarder la doc pour avoir les options optimales\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# calculer le résidu\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "On note aussi `lin.solve` est plus rapide et sa solution est nettement meilleure... `lin.solve` utilise une méthode directe ici. Cela est dû au fait que Scipy utilise la bibliothèque Lapack (qui est imbatable)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Le gradient conjugué de Scipy (avec Lapack)\n", + "\n", + "Le gradient conjugué à tout son sens pour les matrices creuses aussi il est dans la partie \"sparse\" de Scipy.\n", + "On a vu que notre matrice à plus de 99 % de valeur nulles ce qui en fait bien une matrice creuse. Aussi je\n", + "la charge dans le format COO (téléchargez __[Acoo.npz](https://www.lrde.epita.fr/~ricou/cama/data/Acoo.npz)__) qui ne stocke que les valeurs non nulles et" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import scipy.sparse as sparse\n", + "from scipy.sparse.linalg import cg\n", + "\n", + "Ac = sparse.load_npz('/tmp/Acoo.npz')\n", + "%time x = cg(Ac, b)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "On gagne un ordre de grandeur." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ] +}
\ No newline at end of file |