{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Le pivot de Gauss" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Marc Lorenzi\n", "\n", "21 février 2020" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import random\n", "import math\n", "import time\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.rcParams['figure.figsize'] = (8, 4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "L'algorithme du pivot de Gauss est un vaste sujet. Nous allons dans ce notebook nous intéresser à cet algorithme dans un cas particulier, celui des matrices inversibles. \n", "\n", "Soit $A\\in GL_n(\\mathbb R)$ une matrice inversible. Soit $B\\in\\mathcal M_{n,q}(\\mathbb R)$. Considérons l'équation d'inconnue $X\\in\\mathcal M_{nq}(\\mathbb R)$\n", "\n", "$$AX=B$$\n", "\n", "Cette équation admet une unique solution, à savoir $X=A^{-1}B$.\n", "\n", "Nous allons voir comment l'algorithme du pivot permet d'obtenir $X$ efficacement. Nous verrons également deux effets secondaires de cet algorithme :\n", "- Le calcul de $A^{-1}$.\n", "- Le calcul de $\\det A$.\n", "\n", "Dans ce qui va suivre nous travaillerons sur des matrices de flottants, représentées par des listes de listes en Python. __Pour nous conformer aux conventions du langage Python, les indices de lignes et de colonnes des matrices commenceront systématiquement à 0.__\n", "\n", "Avant toute chose, mettons en place quelques fonctions qui nous seront fort utiles." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Quelques fonctions utiles" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.1 Nombre de lignes, de colonnes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Les fonctions `nb_lig` et `nb_col` renvoient le nombre de lignes et colonnes de la matrice $A$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def nb_lig(A): return len(A)\n", "def nb_col(A): return len(A[0])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = [[1, 2, 3], [4, 5, 6]]\n", "print(nb_lig(A), nb_col(A))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.2 Afficher une matrice" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici une fonction qui affiche de façon lisible une matrice de taille raisonnable. Par défaut, les coefficients de $A$ sont affichés avec 5 chiffres après la virgule. Changez ce nombre si cela vous fait plaisir." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def prnt(A):\n", " p = nb_lig(A)\n", " q = nb_col(A)\n", " s = q * ('+' + 12 * '-') + '+'\n", " for i in range(p):\n", " print(s)\n", " for j in range(q):\n", " print('|%+10.5f ' % A[i][j], end='')\n", " print('|')\n", " print(s)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "prnt(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.3 Matrices aléatoires" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Il nous sera fort utile de disposer de matrices \"aléatoires\" pour tester nos fonctions. La fonction `randmat` renvoie une matrice de taille $p\\times q$ dont les coefficients sont des réels suivant une loi uniforme sur l'intervalle $[-r,r]$ (où, par défaut, $r=1$)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def randmat(p, q, r=1):\n", " A = [q * [0] for i in range(p)]\n", " for i in range(p):\n", " for j in range(q):\n", " A[i][j] = random.uniform(-r, r)\n", " return A" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "prnt(randmat(3, 5))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.4 Opérations matricielles" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction `eye` renvoie la matrice identité d'ordre $n$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def eye(n):\n", " A = [n * [0] for i in range(n)]\n", " for i in range(n): A[i][i] = 1\n", " return A" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "prnt(eye(4))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction `prodmat` renvoie le produit des matrices $A$ et $B$. Si $A\\in\\mathcal M_{pq}(\\mathbb R)$ et $b\\in\\mathcal M_{qr}(\\mathbb R)$, la fonction `prodmat` effectue $2pqr$ opérations sur des flottants. En particulier, lorsque $p=q=r=n$, la complexité de cette fonction est $2n^3$. Cette complexité sera notre valeur de référence : nous verrons que l'algorithme du pivot possède une complexité comparable. Résoudre un système, calculer un déterminant, inverser une matrice, sont des opérations dont le coût est analogue à celui de la multiplication matricielle." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def prodmat(A, B):\n", " if nb_col(A) != nb_lig(B):\n", " raise Exception('Matrices incompatibles')\n", " p = nb_lig(A)\n", " q = nb_col(A)\n", " r = nb_col(B)\n", " C = [r * [0] for i in range(p)]\n", " for i in range(p):\n", " for j in range(r):\n", " s = 0\n", " for k in range(q):\n", " s += A[i][k] * B[k][j]\n", " C[i][j] = s\n", " return C" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "prnt(prodmat([[1, 3], [2, 4]], [[5, 7], [6, 8]]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction `submat` renvoie la différence des matrices $A$ et $B$. Elle effectue $pq$ soustractions pour soustraire deux matrices de taille $p\\times q$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def submat(A, B):\n", " p = nb_lig(A)\n", " q = nb_col(A)\n", " if (p != nb_lig(B)) or (q != nb_col(B)):\n", " raise Exception('Matrices incompatibles')\n", " C = [q * [0] for i in range(p)]\n", " for i in range(p):\n", " for j in range(q):\n", " C[i][j] = A[i][j] - B[i][j]\n", " return C" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "prnt(submat([[1, 2, 3], [4, 5, 6]], [[2, 1, 3], [5, 4, 6]]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Opérations sur les lignes d'une matrice" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.1 Introduction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Le principe de l'algorithme du pivot est d'effectuer des opérations sur les lignes d'une matrice. Ces opérations sont au nombre de trois :\n", "\n", "1. Échanger deux lignes.\n", "2. Multiplier une ligne par un réel non nul.\n", "3. Ajouter à une ligne un multiple d'une autre ligne.\n", "\n", "__Pourquoi ce opérations ?__\n", "\n", "__Proposition__ : Soient $A\\in\\mathcal M_{n}(\\mathbb R)$ $X\\in\\mathcal M_{n,q}(\\mathbb R)$ et $B\\in\\mathcal M_{n,q}(\\mathbb R)$. Soit $A'$ la matrice obtenue à partir de $A$ en effectuant l'une des opérations ci-dessus. Soit $B'$ la matrice obtenue en effectuant sur $B$ la même opération que celle faite sur $A$. On a alors\n", "\n", "- ${\\rm rg} A'={\\rm rg}A$\n", "- $AX=B\\Leftrightarrow A'X=B'$\n", "\n", "__Démonstration__ : Faite dans tous les bon cours d'algèbre linéaire.\n", "\n", "On peut de plus remarquer que\n", "\n", "1. Pour l'opération d'échange de lignes, $\\det A'=\\det A$.\n", "2. Pour l'opération multipliant une ligne de $A$ par le réel $\\lambda$, $\\det A'=\\lambda \\det A$\n", "3. Pour l'opération d'ajout à une ligne de $A$ d'un multiple d'une autre ligne de $A$, on a $\\det A'=\\det A$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2 Echange des lignes $i$ et $j$ : $L_i \\leftrightarrow L_j$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $A\\in\\mathcal M_{pq}(\\mathbb R)$. Échanger les lignes $i$ et $j$ de la matrice $A$ nécessite $q$ échanges de coefficients, où $q$ est le nombre de colonnes de $A$. Elle ne fait en revanche aucune opération sur des flottants." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def echanger_lignes(A, i, j):\n", " q = nb_col(A)\n", " for k in range(q):\n", " A[i][k], A[j][k] = A[j][k], A[i][k]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = [[1, 2, 3], [4, 5, 6]]\n", "prnt(A)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "echanger_lignes(A, 0, 1)\n", "prnt(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Remarque__ : Si l'on connaît un peu la façon dont Python implémente les listes, il s'avère que cette opération s'exécute en complexité $O(1)$. Voici la fonction magique `echanger_lignes` :" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def echanger_lignes(A, i, j):\n", " A[i], A[j] = A[j], A[i]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.3 Multiplier la ligne $i$ par le réel $t$ : $L_i \\leftarrow t L_i$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Multiplier la ligne $i$ de $A\\in\\mathcal M_{pq}(\\mathbb R)$ par le réel $t$ nécessite $q$ multiplications de flottants." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def multiplier_ligne(A, i, t):\n", " q = nb_col(A)\n", " for j in range(q): A[i][j] *= t" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "multiplier_ligne(A, 0, 4)\n", "prnt(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.4 Ajouter à la ligne $k$ $t$ fois la ligne $i$ : $L_k \\leftarrow L_k + t L_i$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ajouter à la ligne $k$ de $A\\in\\mathcal M_{pq}(\\mathbb R)$ $t$ fois la ligne $i$ nécessite $2q$ opérations sur des flottants." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def combiner_lignes(A, k, i, t):\n", " q = nb_col(A)\n", " for j in range(q): A[k][j] += t * A[i][j]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sur notre exemple, si nous voulons annuler $A_{00}$, il nous suffit d'effectuer l'opération $L_0\\leftarrow L_0-16L_1$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "combiner_lignes(A, 0, 1, -16)\n", "prnt(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. L'algorithme" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On s'intéresse ici au cas où la matrice $A$ est inversible (et donc carrée). La version que nous allons étudier est appelée __pivot total__. Elle est environ trois fois moins efficace que celle du __pivot partiel__ mais elle est un peu plus simple à décrire. Nous examinerons l'algorithme du pivot partiel dans la dernière partie du notebook." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.1 Introduction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici un petit exemple qui va permettre de comprendre l'algorithme du pivot de Gauss. Résolvons le système\n", "\n", "$$\\left\\lbrace\\begin{array}{ccc}\n", "x+2y+3z&=&8\\\\\n", "4x+5y+6z&=&17\\\\\n", "7x+8y+8z&=&24\n", "\\end{array}\\right.$$\n", "\n", "L'unique solution de ce système est $(0, 1, 2)$. Ce système s'écrit de façon plus compacte $AX=B$, où $A=\\begin{pmatrix}1&2&3\\\\4&5&6\\\\7&8&8\\end{pmatrix}$ et $B=\\begin{pmatrix}8\\\\17\\\\24\\end{pmatrix}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = [[1, 2, 3], [4, 5, 6], [7, 8, 8]]\n", "B = [[8], [17], [24]]\n", "prnt(A)\n", "print()\n", "prnt(B)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Effectuons les opérations $L_1\\leftarrow L_1-4L_0$ et $L_2\\leftarrow L_2-7L_0$ sur les matrices $A$ et $B$. Rappelons que les indices commencent à 0 !" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "combiner_lignes(A, 1, 0, -4)\n", "combiner_lignes(B, 1, 0, -4)\n", "combiner_lignes(A, 2, 0, -7)\n", "combiner_lignes(B, 2, 0, -7)\n", "prnt(A)\n", "print()\n", "prnt(B)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous obtenons un nouveau système $A'X=B'$ où la matrice $A'$ possède des zéros sur la première colonne. La clé de tout est le fait suivant : __$A'$ est inversible et les systèmes $AX=B$ et $A'X=B'$ ont les mêmes solutions.__ \n", "\n", "Maintenant, effectuons l' opération $L_1 \\leftarrow -\\frac 1 3 L_1$ sur $A'$ et $B'$. Ceci place un 1 ligne 1, colonne 1." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "multiplier_ligne(A, 1, -1/3)\n", "multiplier_ligne(B, 1, -1/3)\n", "prnt(A)\n", "print()\n", "prnt(B)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Combinons maintenant les lignes 0 et 2 des matrices $A'$ et $B'$ avec leur ligne 1 pour faire apparaître des 0 sur la colonne 1 de la matrice $A'$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "combiner_lignes(A, 0, 1, -2)\n", "combiner_lignes(B, 0, 1, -2)\n", "combiner_lignes(A, 2, 1, 6)\n", "combiner_lignes(B, 2, 1, 6)\n", "prnt(A)\n", "print()\n", "prnt(B)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous obtenons un nouveau système $A''X=B''$ où $A''$ est inversible, et qui a toujours les mêmes solutions que le système d'origine. Un dernier effort : Opérons maintenant pour faire apparaître des zéros dans la colonne 2 de la matrice $A''$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "multiplier_ligne(A, 2, -1)\n", "multiplier_ligne(B, 2, -1)\n", "combiner_lignes(A, 0, 2, 1)\n", "combiner_lignes(B, 0, 2, 1)\n", "combiner_lignes(A, 1, 2, -2)\n", "combiner_lignes(B, 1, 2, -2)\n", "prnt(A)\n", "print()\n", "prnt(B)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous voici avec un système $A'''X=B'''$ ayant les mêmes solutions que le système d'origine. Comme $A'''=I_3$, ce dernier système est donc extrêmement simple : il s'agit de $X=B'''$. Nous avons résolu notre système d'origine.\n", "\n", "Ce que nous venons de faire sur un exemple peut être fait pour n'importe quel système linéaire $AX=B$ lorsque $A\\in GL_n(\\mathbb R)$ est une matrice inversible. L'algorithme est le suivant :\n", "\n", "Pour $k$ variant de 0 à $n-1$ :\n", "1. Trouver un indice $k\\le \\ell \\le n-1$ tel que $A_{\\ell k}\\ne 0$\n", "2. Si $\\ell\\ne k$ échanger les lignes $\\ell$ et $k$ des matrices $A$ et $B$. Maintenant, $A_{kk}\\ne 0$.\n", "3. Diviser la ligne $k$ de $A$ et $B$ par $\\frac 1 {A_{kk}}$. Maintenant, $A_{kk}=1$.\n", "4. Pour $i$ variant de 0 à $n-1$, $i\\ne k$, combiner la ligne $i$ de $A$ et $B$ avec leur ligne $k$ de sorte que $A_{ik}$ devienne nul. La combinaison est clairement $L_i\\leftarrow L_i - A_{ik}L_k$.\n", "\n", "__Problème__ : __Le__ point problématique est l'étape 1. Un tel indice $\\ell$ existe-t-il ? Nous verrons que oui, grâce à l'inversibilité de $A$. Admettons pour l'instant l'existence de cet indice $\\ell$ et lançons-nous dans l'écriture du code." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.2 Recherche du pivot" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $A\\in\\mathcal M_{n}(\\mathbb R)$. Soit $k\\in [0,n-1|]$. On cherche un indice $\\ell\\in[|k,n-1|]$ tel que $A_{\\ell k}\\ne 0$. Si on le trouve, on le renvoie. Sinon, on lève une exception." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def chercher_pivot(A, k):\n", " n = nb_lig(A)\n", " l = k\n", " while l < n and A[l][k] == 0: l = l + 1\n", " if l == n: raise Exception('Pivot non trouvé')\n", " return l" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = [[0, 2, 3, 4], [1, 5, 1, 1], [4, 2, 2, 0]]\n", "prnt(A)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "chercher_pivot(A, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Revenons à notre question existentielle : un tel indice $\\ell$ existe-t-il ? Nous allons montrer que si la matrice $A$ possède une certaine forme la réponse est oui. Plus loin, lors de la preuve de la correction de l'algorithme du pivot de Gauss, il nous suffira de montrer que $A$ possède effectivement la dite forme.\n", "\n", "__Proposition__ : Soit $A\\in GL_n(\\mathbb R)$. Soit $0\\le k\\le n-1$. On suppose que $A$ s'écrit par blocs\n", "\n", "$$A=\\begin{pmatrix}I_{k}&U\\\\0&V\\end{pmatrix}$$\n", "\n", "où $U\\in\\mathcal M_{k(n-k)}(\\mathbb R)$, $V\\in\\mathcal M_{(n-k)(n-k)}(\\mathbb R)$ et $0\\in\\mathcal M_{(n-k)k}(\\mathbb R)$. Alors il existe $\\ell\\in[0, n-1-k]$ tel que $V_{\\ell k}\\ne 0$.\n", "\n", "__Démonstration__ : Si l'on suppose que la colonne 0 de $V$ est nulle, alors la colonne $k$ de $A$ est combinaison linéaire des colonnes $0,1,\\ldots, k-1$. Ce n'est pas possible puisque $A$ est inversible.\n", "\n", "__Remarque__ : Au lecteur dérouté par le cas $k=0$, on convient qu'une matrice ayant 0 ligne ou 0 colonne est vide. On suppose donc dans ce cas que $A=V$. La colonne 0 de $A$ n'est pas la colonne nulle puisque $A$ est inversible." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.3 Pivoter" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction `pivoter` ci-dessous effectue les opérations suivantes :\n", "\n", "- On cherche un coefficient non nul sur la $k$ème colonne de $A$, à partir de la ligne $k$. Mettons que ce coefficient soit trouvé à la ligne $\\ell$. Un tel coefficient existe si la matrice est inversible.\n", "- On échange les lignes $k$ et $\\ell$ de $A$.\n", "- On divise la ligne $k$ de $A$ par $A_{kk}$ qui est maintenant non nul.\n", "- Pour $i\\ne k$, on soustrait à la $i$ème ligne de $A$ $A_{ik}$ fois la $k$ème ligne.\n", "\n", "- On effectue les mêmes opérations sur le second membre $B$\n", "\n", "Le résultat de cette fonction est d'annuler tous les coefficients de la colonne $k$ de $A$, sauf celui situé à la ligne $k$ qui devient égal à 1. \n", "\n", "La fonction `pivoter` prend également un paramètre $D$ qui est un flottant dont nous reparlerons plus loin." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def pivoter(A, B, k, D):\n", " n = nb_lig(A)\n", " l = chercher_pivot(A, k)\n", " if l != k:\n", " D = -D\n", " echanger_lignes(A, k, l)\n", " echanger_lignes(B, k, l)\n", " P = A[k][k]\n", " D = D * P\n", " multiplier_ligne(B, k, 1 / P)\n", " multiplier_ligne(A, k, 1 / P)\n", " for i in range(n):\n", " if i != k:\n", " Aik = A[i][k]\n", " combiner_lignes(A, i, k, -Aik)\n", " combiner_lignes(B, i, k, -Aik)\n", " return D" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour les amateurs de théorie, voici la\n", "\n", "__Proposition__ : Soit $A\\in GL_n(\\mathbb R)$. Soit $0\\le k\\le n-1$. On suppose que $A$ s'écrit par blocs\n", "\n", "$$A=\\begin{pmatrix}I_{k}&U\\\\0&V\\end{pmatrix}$$\n", "\n", "où $U\\in\\mathcal M_{k(n-k)}(\\mathbb R)$, $V\\in\\mathcal M_{(n-k)(n-k)}(\\mathbb R)$ et $0\\in\\mathcal M_{(n-k)k}(\\mathbb R)$. Soient $A'$ et $B'$ les valeurs des matrices $A$ et $B$ après un appel à `pivoter(A, B, k, D)`. \n", "- On a $A'=\\begin{pmatrix}I_{k+1}&U'\\\\0&V'\\end{pmatrix}$ où $U'\\in\\mathcal M_{(k+1)(n-k-1)}(\\mathbb R)$, $V'\\in\\mathcal M_{(n-k-1)(n-k-1)}(\\mathbb R)$ et $0\\in\\mathcal M_{(n-k-1)(k+1)}(\\mathbb R)$.\n", "- La matrice $A'$ est inversible.\n", "- On a $AX=B$ si et seulement si $A'X=B'$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Reprenons nos matrices d'exemple." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = [[1, 2, 3], [4, 5, 6], [7, 8, 8]]\n", "B = [[8], [17], [24]]\n", "prnt(A)\n", "print()\n", "prnt(B)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "D = 1\n", "D = pivoter(A, B, 0, D)\n", "print('D=', D)\n", "prnt(A)\n", "print()\n", "prnt(B)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "D = pivoter(A, B, 1, D)\n", "print('D=', D)\n", "prnt(A)\n", "print()\n", "prnt(B)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "D = pivoter(A, B, 2, D)\n", "print('D=', D)\n", "prnt(A)\n", "print()\n", "prnt(B)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous avons automatisé ce que nous avions fait manuellement au paragraphe 3.1. La dernière automatisation consiste en une simple boucle." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.4 Une simple boucle" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction `pivot_gauss` est maintenant immédiate. Elle prend en paramètres deux matrices $A$ et $B$. On suppose que\n", "- $A$ est inversible.\n", "- Le nombre de lignes de $B$ est le même que celui de $A$.\n", "\n", "Soit $n$ le nombre de lignes de $A$. Pour $k$ variant de 0 à $n - 1$, on applique la fonction `pivoter` à la matrice $A$ et à la matrice $B$.\n", "\n", "__Proposition__ : Soient $A\\in GL_n(\\mathbb R)$ et $B\\in\\mathcal M_{nq}(\\mathbb R)$. Un appel à `pivot_gauss(A, B)` renvoie un triplet $(A',X, D)$ où\n", "- $A'=I_n$.\n", "- $X\\in \\mathcal M_{nq}(\\mathbb R)$ est l'unique solution de l'équation $AX=B$.\n", "- $D=\\det A$ est le déterminant de la matrice $A$.\n", "\n", "__Démonstration__ : Tout a déjà été (presque) fait. Le lecteur consciencieux reprendra les propositions précédentes pour s'en convaincre." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def pivot_gauss(A, B):\n", " n = nb_lig(A)\n", " A = [A[i].copy() for i in range(n)]\n", " B = [B[i].copy() for i in range(n)]\n", " D = 1\n", " for i in range(n):\n", " D = pivoter(A, B, i, D)\n", " return (A, B, D)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Reprenons une dernière fois notre petit exemple" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = [[1, 2, 3], [4, 5, 6], [7, 8, 8]]\n", "B = [[8], [17], [24]]\n", "A1, B1, D = pivot_gauss(A, B)\n", "print(\"D =\", D)\n", "prnt(A1)\n", "print()\n", "prnt(B1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.5 Complexité" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Combien d'opérations sur des flottants la fonction `pivot_gauss` effectue-t-elle lorsqu'elle est appelée sur une matrice $A\\in GL_n(\\mathbb R)$ et une matrice $B\\in\\mathcal M_{n,q}(\\mathbb R)$ ?\n", "\n", "La fonction `pivoter` est appelée $n$ fois. À chaque appel de celle-ci :\n", "\n", "- On cherche la valeur du pivot par un appel à `chercher_pivot` : cette fonction ne fait aucune opération sur des flottants.\n", "- On échange éventuellement deux lignes de $A$ et $B$ : toujours aucune opération.\n", "- On multiplie une ligne de $A$ et une ligne de $B$ par un réel : $n + q$ opérations.\n", "- On appelle $n-1$ fois la fonction `combiner_lignes` sur $A$ et $B$ : $2(n+q)$ opérations.\n", "\n", "Notons $C_{n, q}$ le nombre d'opérations total. On a donc\n", "\n", "$$C_{n,q}=n(n+q + 2(n-1)(n+q)=n(n+q)(2n-1)$$\n", "\n", "Deux cas nous intéresseront dans la suite :\n", "\n", "- On résout un système $AX=B$ où $B\\in\\mathcal M_{n1}(\\mathbb R)$. Dans ce cas $q=1$ et\n", "\n", "$$C_{n1}=n(n+1)(2n-1)=2n^3+n^2-n\\sim2n^3$$\n", "\n", "- On inverse une matrice $A$. Nous verrons que dans ce cas $B=I_n\\in\\mathcal M_{n}(\\mathbb R)$. Ainsi $q=n$ et\n", "\n", "$$C_{nn}=n(2n)(2n-1)=4n^3-2n^2\\sim 4n^3$$\n", "\n", "__Remarque__ : Le nombre d'opérations nécessaire pour multiplier deux matrices $n\\times n$ avec notre fonction `prodmat` est $2n^3$. Ainsi, résoudre un système de Cramer coûte exactement la même chose (en termes d'opérations sur des flottants) que multiplier deux matrices. Et contrairement aux idées reçues, inverser une matrice ne coûte que deux fois plus cher ! " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tentons une expérience numérique. Résolvons un système de 200 équations à 200 inconnues. Quel est le temps nécessaire à sa résolution ? Comment ce temps ce compare-t-il au temps mis pour élever une matrice de taille 200 au carré ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n = 200\n", "A = randmat(n, n)\n", "B = randmat(n, 1)\n", "t1 = time.time()\n", "pivot_gauss(A, B)\n", "t2 = time.time()\n", "print('Temps pivot : %.1fs' % (t2 - t1))\n", "t3 = time.time()\n", "prodmat(A, A)\n", "t4 = time.time()\n", "print('Temps produit: %.1fs' % (t4 - t3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Les deux temps sont comparables, et nous sommes même plutôt agréablement surpris. Rappelons que notre complexité théorique ne compte que les opérations sur des flottants. Il s'avère que des instructions du type `A[i][j] = ...` n'ont pas un coût négligeable en Python. Nous ne creuserons pas plus avant pour ne pas rendre ce notebook incommensurablement long.\n", "\n", "Disons que notre complexité « pratique » est $C_n\\sim Cn^3$, exprimée en secondes. Que vaut $C$ ? Cela dépend de la machine utilisée, bien évidemment." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n = 200\n", "A = randmat(n, n)\n", "B = randmat(n, 1)\n", "t1 = time.time()\n", "pivot_gauss(A, B)\n", "t2 = time.time()\n", "C = (t2 - t1) / n ** 3\n", "print('C=%.2e' % C)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Traçons la courbe de la fonction $n\\longmapsto \\frac{t}{Cn^3}$ où $t$ est le temps mis pour résoudre un système $n\\times n$. Traçons également la courbe de la fonction $n\\longmapsto \\frac{t'}{Cn^3}$ où $t'$ est le temps mis pour multiplier deux matrices $n\\times n$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "s1 = []\n", "s2 = []\n", "rng = range(30, 201, 10)\n", "for n in rng:\n", " A = randmat(n, n)\n", " B = randmat(n, 1)\n", " t1 = time.time()\n", " pivot_gauss(A, B)\n", " t2 = time.time()\n", " s1.append((t2 - t1) / (C * n ** 3))\n", " t3 = time.time()\n", " prodmat(A, A)\n", " t4 = time.time()\n", " s2.append((t4 - t3) / (C * n ** 3))\n", "plt.ylim(0, 2)\n", "plt.plot(rng, s1, 'k')\n", "plt.plot(rng, s2, 'r')\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous obtenons ce que nous pensions. La courbe noire, celle qui correspond à l'algorithme du pivot, est quasi-constante, égale à 1. La courbe rouge (produit matriciel) est elle aussi-quasi constante, avec une valeur plus faible que 1, environ 0.8." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.6 Quelques justifications sur la valeur de $D$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $D'$ la valeur de $D$ après l'appel à `pivoter`. Soit $A'$ la valeur de la matrice $A$ après ce même appel.\n", "- Un éventuel échange de lignes change le signe de $\\det A$.\n", "- La multiplication d'une ligne de $A$ par $\\frac 1 P$ divise le déterminant de $A$ par $P$.\n", "- L'ajout à une ligne de $A$ d'un multiple d'une autre ligne de $A$ ne change pas son déterminant.\n", "\n", "Ainsi, en posant $K=\\pm P$,\n", "\n", "$$\\det A'=\\frac 1 K\\det A,\\quad D'=KD$$\n", "\n", "d'où \n", "\n", "$$D'\\det A'=D\\det A$$\n", "\n", "Ainsi, la quantité $D\\det A$ est un invariant de boucle. Après $n$ itérations de `pivoter`, la matrice $A$ devient la matrice identité. Soit $D$ la valeur finale. On a donc $D\\det I_n = 1\\det A$, ou encore\n", "\n", "$$D=\\det A$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Inversion d'une matrice" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "L'algorithme du pivot permet d'inverser une matrice $A\\in\\mathcal M_n(\\mathbb R)$ : il suffit de mettre pour second membre du système la matrice $I_n$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def inverse(A, dbg=False):\n", " p = nb_lig(A)\n", " q = nb_col(A)\n", " if p != q:\n", " raise Exception('inverser: matrice pas carrée')\n", " B = eye(p)\n", " A1, B1, D = pivot_gauss(A, B)\n", " if dbg: print('Déterminant: %.5e' % D)\n", " return B1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Inversons une matrice aléatoire $4\\times 4$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = randmat(4, 4)\n", "prnt(A)\n", "print()\n", "B = inverse(A)\n", "prnt(B)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Faisons le produit $AB$, cela devrait donner l'identité. En fin presque : je rappelle que nous faisons des calcul en flottants, et donc des calculs approchés." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "prnt(prodmat(A, B))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Faisons la même expérience qu'avec le pivot de Gauss : Comparons le temps mis pour inverser une matrice avec celui mis pour l'élever au carré." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n = 200\n", "A = randmat(n, n)\n", "t1 = time.time()\n", "B = inverse(A)\n", "t2 = time.time()\n", "print('Temps inversion : %.1fs' % (t2 - t1))\n", "t3 = time.time()\n", "B = prodmat(A, A)\n", "t4 = time.time()\n", "print('Temps produit : %.1fs' % (t4 - t3))\n", "print('Quotient des temps: %.1f' % ((t2 - t1) / (t4 - t3)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notre théorie prévoyait un temps d'inversion double du temps de multiplication. La pratique nous conforte dans l'idée que contrairement à ce que disent les mauvaises langues, la théorie, ça marche :-)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n = 200\n", "A = randmat(n, n)\n", "t1 = time.time()\n", "B = inverse(A)\n", "t2 = time.time()\n", "C = (t2 - t1) / n ** 3\n", "print('C = %.2e' % C)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "s1 = []\n", "s2 = []\n", "rng = range(10, 201, 10)\n", "for n in rng:\n", " A = randmat(n, n)\n", " B = randmat(n, 1)\n", " t1 = time.time()\n", " inverse(A)\n", " t2 = time.time()\n", " s1.append((t2 - t1) / (C * n ** 3))\n", " t3 = time.time()\n", " prodmat(A, A)\n", " t4 = time.time()\n", " s2.append((t4 - t3) / (C * n ** 3))\n", "plt.ylim(0, 2)\n", "plt.plot(rng, s1, 'k')\n", "plt.plot(rng, s2, 'r')\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5. Stabilité de l'algorithme du pivot" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Cette section mériterait un long développement ... que nous ne ferons pas." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 5.1 Norme de Frobenius" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous avons eu, lors de la vérification de nos algorithmes, à nous poser des questions du genre : A-t-on $AX=B$ ? Le souci est que, opérant sur des flottants, la vraie question serait plutôt : a-t-on $AX\\simeq B$ ? Ou encore : a-t-on $AX-B$ \"petit\" ? La notion de petitesse est très relative. Introduisons, pour toute matrice $A$ à coefficients réels, la quantité\n", "\n", "$$\\Vert A\\Vert = \\left(\\sum_{i,j}A_{ij}^2\\right)^{1/2}$$\n", "\n", "Le nombre $\\Vert A \\Vert$ est la __norme de Frobenius__ de la matrice $A$. On montre aisément que la fonction $A\\mapsto\\Vert A \\Vert$ est une norme sur l'espace vectoriel $\\mathcal M_{pq}(\\mathbb R)$. Une \"petite \" matrice sera donc, dans la discussion qui va suivre, une matrice de \"petite norme\"." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def norme_frob(A):\n", " s = 0\n", " p = nb_lig(A)\n", " q = nb_col(A)\n", " for i in range(p):\n", " for j in range(q):\n", " s += A[i][j] ** 2\n", " return math.sqrt(s)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(norme_frob([[1, 2], [3, 4]]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction `check_prod` prend en paramètres trois matrices $A, X, B$. Elle renvoie $\\Vert AX-B\\Vert$. Si $X$ est solution du système $AX=B$, le valeur renvoyée par la fonction `check_prod` est donc 0. Dans la pratique, à cause des erreurs d'arrondis, la valeur devrait être petite quoique non nulle." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def check_prod(A, X, B):\n", " return norme_frob(submat(prodmat(A, X), B))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = randmat(100, 100)\n", "B = randmat(100, 1)\n", "A1, X, D = pivot_gauss(A, B)\n", "print(check_prod(A, X, B))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Il semble bien que $AX\\simeq B$. Tentons autre chose. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 5.2 Les matrices de Hilbert" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Les matrices de Hilbert sont connues pour poser de gros problèmes aux algorithmes numériques. Pour tout $n\\ge 1$, posons\n", "\n", "$$H_n=\n", "\\begin{pmatrix}\n", "\\frac 1 1&\\frac 1 2&\\frac 1 3&\\ldots&\\frac 1 n\\\\\n", "\\frac 1 2&\\frac 1 3&\\frac 1 4&\\ldots&\\frac 1 {n+1}\\\\\n", "\\frac 1 3&\\frac 1 4&\\frac 1 5&\\ldots&\\frac 1 {n+2}\\\\\n", "\\vdots&\\vdots&\\vdots&&\\vdots\\\\\n", "\\frac 1 n&\\frac 1 {n+1}&\\frac 1 {n+2}&\\ldots&\\frac 1 {2n-1}\\\\\n", "\\end{pmatrix}$$\n", "\n", "On peut montrer que $H_n$ est inversible, quoique son déterminant soit très petit, vraiment TRÈS petit." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def hilbert(n):\n", " return [[1 / (i + j + 1) for i in range(n)] for j in range(n)]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "prnt(hilbert(8))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Prenons une matrice de Hilbert $A$ de taille 20, calculons son inverse $B$ avec notre fonction `inverse`, puis vérifions si $AB=I_{20}$. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n = 20\n", "A = hilbert(n)\n", "B = inverse(A, dbg=True)\n", "print(check_prod(A, B, eye(n)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Le résultat est catastrophique. Nous ne creuserons pas plus avant ce sujet, il nécessiterait un notebook à lui tout seul. Mais il est bon de garder à l'esprit que sur certaines matrices l'algorithme du pivot peut tout à fait donner des résultats __désastreux__." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6. Le pivot partiel" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous allons dans cette section améliorer la complexité de l'algorithme du pivot en étudiant la métode du __pivot partiel__. Nous verrons que par cette méthode les temps de calcul sont divisés par 3." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.1 Résolution d'un système triangulaire" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $A\\in\\mathcal M_n(\\mathbb R)$ une matrice triangulaire supérieure. Supposons que les coefficients diagonaux de $A$ sont non nuls. La matrice $A$ est alors inversible. Soit $B\\in\\mathcal M_{n1}(\\mathbb R)$. Il est alors aisé de résoudre le système $AX=B$, d'inconnue $X\\in\\mathcal M_{n1}(\\mathbb R)$. En effet, $AX=B$ si et seulement si \n", "\n", "$$A_{(n-1)(n-1)}x_{n-1}=b_{n-1}$$\n", "\n", "et pour tout $i\\in [0, n-2]$,\n", "\n", "$$A_{ii}x_i+\\sum_{j=i+1}^{n-1}A_{ij}x_j=b_i$$\n", "\n", "ou encore\n", "\n", "$$x_i= \\frac 1 {A_{ii}}\\left(b_i-\\sum_{j=i+1}^{n-1}A_{ij}x_j\\right)$$\n", "\n", "On obtient donc immédiatement $x_{n-1}$. Puis, supposant calculés $x_{n-1},\\ldots,x_{i+1}$, on en déduit $x_i$ en fonction de ceux-ci.\n", "\n", "La fonction `solve_triang` prend en paramètres une matrice triangulaire inversible $A$ et une matrice colonne $B$. Elle modifie $B$ sur place. En sortie de fonction, $B$ contient l'unique solution du système $AX=B$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def solve_triang(A, B):\n", " n = nb_lig(A)\n", " for i in range(n - 1, -1, -1):\n", " s = B[i][0]\n", " for j in range(i + 1, n):\n", " s -= A[i][j] * B[j][0]\n", " B[i][0] = s / A[i][i]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Quelle est la complexité $C_n$ de `solve_triang` en nombre d'opérations sur des flottants ? Pour chaque $i$ compris entre 0 et $n-1$, la fonction effectue $2(n-i-1)+1$ opérations. Ainsi,\n", "\n", "$$C_n=\\sum_{i=0}^{n-1}\\left(2(n-i-1)+1\\right)=\\sum_{k=0}^{n-1}\\left(2k+1\\right)=n^2$$\n", "\n", "Résoudre un système triangulaire a donc un coût négligeable par rapport à celui de la résolution d'un système quelconque." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.2 L'idée du pivot partiel" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Le faible coût de la résolution d'un système triangulaire suggère de modifier le pivot de Gauss de la façon suivante : au lieu d'opérer sur les lignes de la matrice $A$ pour obtenir la matrice identité, \n", "- On amène $A$ à une forme triangulaire supérieure.\n", "- On résout le système obtenu par l'algorithme du paragraphe précédent.\n", "Nous allons voir que le gain de temps est conséquent." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.3 Nouvelles opérations sur les lignes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Les nouvelles fonctions sur les lignes se passent de commentaires. \n", "- La fonction `echanger_lignes` est gardée telle quelle, son coût est négligeable par rapport au coût total.\n", "- La fonction `multiplier_ligne`, appelée sur la ligne $i$ d'une matrice, ne modifie les colonnes de la matrice qu'à partir de la colonne $i$. L'idée est que les colonnes précédentes contiendront des zéros lorsque cette fonction sera appelée.\n", "- Même remarque pour la fonction `combiner_lignes`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def multiplier_ligne_partiel(A, i, t):\n", " q = nb_col(A)\n", " for j in range(i, q): A[i][j] *= t" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def combiner_lignes_partiel(A, k, i, t):\n", " q = nb_col(A)\n", " for j in range(i, q): A[k][j] += t * A[i][j]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.4 La fonction de pivot" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La nouvelle fonction de pivot est un quasi copier-coller de l'ancienne. Les seules différences sont :\n", "- Des appels à des opérations partielles sur les lignes de la matrice $A$.\n", "- Et surtout une modification des lignes de $A$ uniquement à partir de la ligne $k+1$ : les lignes 0 à $k-1$ sont laissées inchangées, c'est ce point précis qui amène $A$ à une forme triangualire, et pas diagonale. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def pivoter_partiel(A, B, k, D):\n", " n = nb_lig(A)\n", " l = chercher_pivot(A, k)\n", " if l != k:\n", " D = -D\n", " echanger_lignes(A, k, l)\n", " echanger_lignes(B, k, l)\n", " P = A[k][k]\n", " D = D * P\n", " multiplier_ligne(B, k, 1 / P)\n", " multiplier_ligne_partiel(A, k, 1 / P)\n", " for i in range(k + 1, n):\n", " Aik = A[i][k]\n", " combiner_lignes(B, i, k, -Aik)\n", " combiner_lignes_partiel(A, i, k, -Aik)\n", " return D" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notre nouvelle fonction `pivot_gauss` appelle `pivoter_partiel` au lieu de `pivoter` et n'oublie pas de résoudre le système triangulaire obtenu." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def pivot_gauss_partiel(A, B):\n", " n = nb_lig(A)\n", " A = [A[i].copy() for i in range(n)]\n", " B = [B[i].copy() for i in range(n)]\n", " D = 1\n", " for k in range(n):\n", " D = pivoter_partiel(A, B, k, D)\n", " solve_triang(A, B)\n", " return (A, B, D)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Testons." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = randmat(100, 100)\n", "B = randmat(100, 1)\n", "A1, B1, D = pivot_gauss_partiel(A, B)\n", "print(check_prod(A, B1, B))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tout a l'air de fonctionner." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.5 Complexité" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Quelle est la complexité $C_n$ de `pivot_gauss_partiel` en termes d'opérations sur des flottants ? Supposons la matrice $A$ carrée de taille $n$, le second membre $B$ possédant $n$ lignes et 1 colonne.\n", "\n", "Pour $k=0,\\ldots,n-1$, on appelle `pivoter_partiel(A, B, k, D)`. Cette fonction effectue :\n", "- 1 multiplication sur un coefficient de $B$\n", "- $n-k$ multiplications sur la ligne $k$ de $A$.\n", "- Pour $i=k+1,\\ldots,n-1$, $2(n-i)$ opérations sur la ligne $i$ de $A$ et 2 opérations sur le $i$ième coefficient de $B$. \n", "\n", "Il faut ajouter à cela la résolution du système triangulaire, qui nécessite $n^2$ opérations. Ainsi,\n", "\n", "$$\\begin{array}{lll}\n", "C_n&=&\\sum_{i=0}^{n-1}\\left(1+(n-i)+\\sum_{k=i+1}^{n-1}(2(n-i)+2)\\right) + n^2\\\\\n", "&=&\\sum_{i=0}^{n-1}\\left(n-i+1+(n-1-i)(2(n-i)+2)\\right)+n^2\n", "\\end{array}$$\n", "\n", "Posons $k=n-i-1$. Il vient\n", "\n", "$$\\begin{array}{lll}\n", "C_n&=&\\sum_{k=0}^{n-1}\\left(k+2+k(2k+4)\\right)+n^2\\\\\n", "&=&\\sum_{k=0}^{n-1}\\left(2k^2+5k+2)\\right)+n^2\\\\\n", "&=&\\frac 1 3n(n-1)(2n-1)+\\frac 5 2n(n-1)+n^2+2n\\\\\n", "&=&\\frac 1 6n(n-1)(4n+13)+n^2+2n\\\\\n", "&=&\\frac 2 3 n^3+\\frac 5 2 n^2-\\frac 1 6 n\\\\\n", "&\\sim&\\frac 2 3 n^3\n", "\\end{array}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rappelons-nous que par la méthode du pivot total, la complexité était $2n^3$. Nous devrions donc obtenir un gain d'un facteur 3 avec l'algorithme du pivot partiel. Quelques vérifications s'imposent ..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.5 Quelques tests" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Commençons par vérifier que nos fonctions fonctionnent. Résolvons un système de 6 équtaions à 6 inconnues." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = randmat(6, 6)\n", "B = randmat(6, 1)\n", "prnt(A)\n", "print()\n", "A1, B1, D = pivot_gauss_partiel(A, B)\n", "prnt(A1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La matrice $A_1$ possède bien une forme triangulaire, ce qui est rassurant. Maintenant, la matrice $B_1$ devrait contenir la solution (unique) du système $AX=B$. Calculons donc $AB_1-B$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A1, X, D = pivot_gauss_partiel(A, B)\n", "print(check_prod(A, X, B))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tout a l'air de fonctionner. Passons maintenant à quelques tests sur la complexité de l'algorithme, comparée à celle du pivot total. Résolvons un système de 300 équations et comparons les temps pour le pivot partiel et le pivot total." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n = 300\n", "A = randmat(n, n)\n", "B = randmat(n, 1)\n", "t1 = time.time()\n", "A1, B1, D = pivot_gauss_partiel(A, B)\n", "t2 = time.time()\n", "print('Pivot partiel: %.2fs' % (t2 - t1))\n", "t3 = time.time()\n", "A1, B1, D = pivot_gauss(A, B)\n", "t4 = time.time()\n", "print('Pivot total: %.2fs' % (t4 - t3))\n", "print('Quotient: %.2f' % ((t4 - t3) / (t2 - t1)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Des calculs un peu plus poussés montrent que le quotient des complexités théoriques est \n", "\n", "$$3-\\frac{33}{4n}+\\frac{149}{16n^2}o\\left(\\frac 1 {n^2}\\right)$$\n", "\n", "Il suffit pour cela de reprendre leur valeur et de faire un développement asymptotique. Qu'obtenons-nous pour $n=200$ ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def quot(n):\n", " return 3 - 33 / (4 * n) + 149 / (16 * n ** 2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(quot(300))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La pratique et la théorie s'accordent vraiment bien. Cela justifie _a posteriori_ le choix de compter les opérations sur des flottants pour évaluer la complexité. \n", "\n", "Allons, un dernier test ! Résolvons des systèmes de $n$ équations à $n$ inconnues pour $n=30,40,50,\\ldots,200$. Traçons le quotient des temps pivot total / pivot partiel en fonction de $n$. Traçons également en rouge notre estimation théorique du quotient." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "s1 = []\n", "s2 = []\n", "rng = range(30, 201, 10)\n", "for n in rng:\n", " A = randmat(n, n)\n", " B = randmat(n, 1)\n", " t1 = time.time()\n", " A1, B1, D = pivot_gauss_partiel(A, B)\n", " t2 = time.time()\n", " t3 = time.time()\n", " A1, B1, D = pivot_gauss(A, B)\n", " t4 = time.time()\n", " s1.append((t4 - t3) / (t2 - t1))\n", " s2.append(quot(n))\n", "plt.ylim(1, 5)\n", "plt.plot(rng, s1, 'k')\n", "plt.plot(rng, s2, 'r')\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notre estimation théorique est tout à fait satisfaisante ... Si jamais il reste un lecteur pour lire ces lignes je lui dis bravo :-) ..." ] } ], "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.8" } }, "nbformat": 4, "nbformat_minor": 4 }