diff --git "a/Chapitre 9 - Produits scalaires et espaces euclidens/9.13-9.14 La factorisation QR - application \303\240 la r\303\251solution d'un syst\303\250me au sens des moindres carr\303\251es.ipynb" "b/Chapitre 9 - Produits scalaires et espaces euclidens/9.13-9.14 La factorisation QR - application \303\240 la r\303\251solution d'un syst\303\250me au sens des moindres carr\303\251es.ipynb" index 305a910..992f09b 100644 --- "a/Chapitre 9 - Produits scalaires et espaces euclidens/9.13-9.14 La factorisation QR - application \303\240 la r\303\251solution d'un syst\303\250me au sens des moindres carr\303\251es.ipynb" +++ "b/Chapitre 9 - Produits scalaires et espaces euclidens/9.13-9.14 La factorisation QR - application \303\240 la r\303\251solution d'un syst\303\250me au sens des moindres carr\303\251es.ipynb" @@ -1,371 +1,371 @@ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# **Concept(s)-clé(s) et théorie**\n", "\n", "## Définition 1\n", "Soit $A \\in \\mathcal{M}_{m \\times n}(\\mathbb{R})$ une matrice dont les colonnes sont linéairement indépendantes (vues comme vecteurs de $\\mathbb{R}^m$). Alors il existe une factorisation du type $A = QR$, où $Q$ est une matrice $m \\times n$ dont les colonnes forment une base orthonormée de l'espace colonnes de $A$ et $R$ est une matrice $n \\times n$ triangulaire supériore, inversible, dont les coefficients diagonaux sont strictement positifs.\n", "\n", "## Algorithme 1\n", "Soit $A \\in \\mathcal{M}_{m \\times n}(\\mathbb{R})$ une matrice dont les colonnes sont linéairement indépendantes (vues comme vecteurs de $\\mathbb{R}^m$). Afin de déterminer une factorisation $QR$ de $A$, on procède comme suit:\n", "1. Poser $\\mathcal{C} = \\left( c_1, \\dots, c_n \\right)$, qui ici forme une base de l'espace colonnes $W$ de $A$\n", "2. A l'aide du procède de Gram-Schmidt, calculer une base orthoormée $\\mathcal{B} = \\left( w_1, \\dots, w_n \\right)$ de $W$\n", "3. Définir $Q \\in \\mathcal{M}_{m \\times n}(\\mathbb{R})$ comme étant la matrice dont la $i$-ème colonne est $w_i$\n", "4. Pour tout $1 \\leq k \\leq n$, écrire $c_k = r_{1k}w_1 + r_{2k}w_2 + \\dots + r_{kk}w_k + 0w_{k+1} + \\dots + 0w_{n}$. On supposerà que $r_{ii} \\geq 0$, quitte à replacer $w_i$ par $-w_i$. Poser alors $r_k = \\begin{pmatrix} r_{1k} & r_{2k} & \\dots & r_{kk} & 0 & \\dots & 0 \\end{pmatrix}^T$\n", "5. Définir $R \\in \\mathcal{M}_{n \\times n}(\\mathbb{R})$ comme étant la matrice dont la $i$-ème colonne est $r_i$.\n", "\n", "Donc:\n", "\n", "$$ Q = \\left( \\begin{array}{@{}c|c|c|c@{}} w_1 & w_2 & \\dots & w_n \\end{array}\\right) \\qquad \n", "R = \\begin{pmatrix} \\langle c_1, w_1 \\rangle & \\langle c_2, w_1 \\rangle & \\langle c_3, w_1 \\rangle & \\dots & \\langle c_n, w_1 \\rangle \\\\\n", "0 & \\langle c_2, w_2 \\rangle & \\langle c_3, w_2 \\rangle & \\dots & \\langle c_n, w_2 \\rangle \\\\ \n", "0 & 0 & \\langle c_3, w_3 \\rangle & \\dots & \\langle c_n, w_2 \\rangle \\\\ \n", "\\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\ \n", "0 & 0 & 0 & \\dots & \\langle c_n, w_n \\rangle \\end{pmatrix}$$\n", "\n", "## Proposition 1\n", "Soit $Q \\in \\mathcal{M}_{m \\times n}(\\mathbb{R})$ une matrice dont les colonnes forment une base orthonormée de l'espace des colonnes $W$ de $Q$. Alors \n", "\n", "$$QQ^Tb = proj_Wb$$ \n", "\n", "pour tout $b \\in \\mathcal{M}_{m \\times 1}(\\mathbb{R})$.\n", "\n", "## Proposition 2\n", "Soit $A \\in \\mathcal{M}_{m \\times n}(\\mathbb{R})$ une matrice dont les colonnes sont linéairement indépendantes et soit $A = QR$ une factorisation $QR$ de $A$. Alors pour tout $b \\in \\mathcal{M}_{m \\times 1}(\\mathbb{R})$, l'equation $Ax=b$ admet une unique solution au sens du moindres carrées, donnée par la formule\n", "\n", "$$ \\hat{x} = R^{-1} Q^T b$$\n", "\n", "Également, la solution du système au sens du moindres carrée peut être calculée en résolvant le système linéaire traingulaire supérieur suivant\n", "\n", "$$ R\\hat{x} = Q^T b$$\n", "\n", "\n", "\n", "En utilisant la méthode de substitution vers l'arrière, il n'est alors pas nécessaire de dériver explicitement l'espression de l'inverse de R." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercises et Examples" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import Librairie.AL_Fct as al\n", "import Corrections.corrections as corrections\n", "import numpy as np\n", "import sympy as sp\n", "import time\n", "from scipy.linalg import solve_triangular, lu" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice 1\n", "Considérez les matrices suivantes. En tirant parti des cellules interactives suivantes, appliquez l'algorithme de Gram-Schmidt et dérivez leur factorisation QR.\n", "\n", "1. $\\qquad A_1 = \\begin{pmatrix} 1 & 2 & 0 \\\\ 0 & 1 & 1 \\\\ 1 & 0 & 1 \\end{pmatrix}$\n", "2. $\\qquad A_2 = \\begin{pmatrix} 1 & 1 \\\\ 0 & 2 \\\\ 0 & -2 \\\\ 0 & -1 \\end{pmatrix}$\n", "3. $\\qquad A_3 = \\begin{pmatrix} -3 & 1 & 0 \\\\ 1 & 1 & 1 \\\\ 5 & 0 & 2 \\\\ -1 & -1 & 0 \\end{pmatrix}$ \n", "4. $\\qquad A_4 = \\begin{pmatrix} -4 & 20 & 35 & 5 \\\\ -4 & -30 & -15 & 55 \\\\ -8 & 40 & -80 & -65 \\\\ 23 & -15 & 30 & 15 \\end{pmatrix}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Instructions\n", "\n", "Pour utiliser la méthode interactive pour la factorisation QR, procédez comme suit:\n", "\n", "1. Insérez le numéro de cas souhaité dans la cellule appelèe \"SÉLECTION DU CAS\". Exécutez le cellules appelées \"SÉLECTION DU CAS\" et \"INITIALISATION DES VARIABLES\"\n", "2. Exécutez la cellule appelée \"SÉLECTION DES PARAMÈTRES\" pour insérer les coefficients necessaires\n", "3. Exécutez la cellule appelée \"EXÉCUTER L'ÉTAPE DE L'ALGORITHME\" pour exécuter l'étape de l'algorithme avec les paramètres précédemment sélectionnés\n", "4. Répétez les étapes 2 et 3 jusqu'à ce que la factorisation soit terminée" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SÉLECTION DU CAS\n", "case_number = 1 # CHOISISSEZ LE NUMÉRO DE CAS ICI ET EXECUTEZ LA CELLULE SUIVANTE!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# INITIALISATION DES VARIABLES\n", "if case_number == 1: \n", " A = np.array([[1,2,0], [0,1,1], [1,0,1]])\n", "elif case_number == 2:\n", " A = np.array([[1,1], [0,2], [0,-2], [0,-1]])\n", "elif case_number == 3:\n", " A = np.array([[-3,1,0], [1,1,1], [5,0,2], [-1,-1,0]])\n", "elif case_number == 4:\n", " A = np.array([[-4,20,35,5], [-4,-30,-15,55], [-8,40,-80,-65], [23,-15,30,15]])\n", "else:\n", " print(f\"{case_number} n'est pas un numéro de cas valide!\" \n", " f\"Numéros de cas disponibles: [1,2,3,4]\")\n", " \n", "Q = np.full(A.shape, np.nan)\n", "R = np.full((A.shape[1],A.shape[1]), np.nan)\n", "QList = [Q]\n", "RList = [R]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SÉLECTION DES PARAMÈTRES\n", "norm_coeff, proj_coeffs = al.manual_QR(QList[-1], RList[-1], A)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# EXÉCUTER L'ÉTAPE DE L'ALGORITHME\n", "al.interactive_QR(norm_coeff, proj_coeffs, QList, RList, A)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "8/25" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice 2\n", "\n", "Considérez les systèmes linéaires suivants, avec forme $Ax=b$. Calculez leur solution au sens des moindres carrés en tirant parti de la factorisation QR de A. \n", "\n", "1. $\\quad A = \\begin{pmatrix}1 & 1 \\\\ 1 & -1\\end{pmatrix} \\qquad \\qquad \\ \\ b = \\begin{pmatrix}1 \\\\ -1\\end{pmatrix}$\n", "2. $\\quad A = \\begin{pmatrix}-2 & 1 \\\\ 1 & 3 \\\\ 2 & 0\\end{pmatrix} \\qquad \\qquad \\quad b = \\begin{pmatrix}1 \\\\ 2 \\\\ 0\\end{pmatrix}$\n", "3. $\\quad A = \\begin{pmatrix}1 & 0 & 0\\\\ 2 & 0 & 1 \\\\ 0 & 4 & 0 \\\\ 0 & 3 & 2 \\\\ 2 & 0 & 2 \\end{pmatrix} \\qquad \\qquad \\ b = \\begin{pmatrix}0 \\\\ 1 \\\\ 0 \\\\ -1 \\\\ 1 \\end{pmatrix}$\n", "4. $\\quad A = \\begin{pmatrix}5 & 0 & 0\\\\ -1 & -2 & -2 \\\\ 1 & 1 & 2 \\\\ 3 & 2 & 0 \\\\ 0 & 0 & 1\\end{pmatrix} \\qquad \\ \\ b = \\begin{pmatrix}0 \\\\ 2 \\\\ 0 \\\\ 1 \\\\ 4\\end{pmatrix}$\n", "\n", "### !! ATTENTION !!\n", "Si $A$ ne correspond pas aux exigences de la *Définition 1*, insérez **None** pour les valeurs des $Q$, $R$ et de la solution $\\hat{x}$ au système linéaire donnée au sens des moindres carrés!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SÉLECTION DU CAS\n", "case_number = 1 # CHOISISSEZ LE NUMÉRO DE CAS ICI ET EXECUTEZ LA CELLULE SUIVANTE!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# INITIALISATION DES VARIABLES\n", "if case_number == 1:\n", " A_cols = [[1,1], [1,-1]]\n", " A = np.array(A_cols).T\n", " b = [1,-1]\n", " dim=2\n", "elif case_number == 2:\n", " A_cols = [[-2,1,2], [1,3,0]]\n", " A = np.array(A_cols).T\n", " b = [1,2,0]\n", " dim=2\n", "elif case_number == 3:\n", " A_cols = [[1,2,0,0,2], [0,0,4,3,0], [0,1,0,2,2]]\n", " A = np.array(A_cols).T\n", " b = [0,1,0, -1, 1]\n", " dim=3\n", "elif case_number == 4:\n", " A_cols = [[5,-1,1,3,0], [0,-2,1,2,0], [0,-2,2,0,1]]\n", " A = np.array(A_cols).T\n", " b = [0,2,0,1,4]\n", " dim=3\n", "else:\n", " print(f\"{case_number} n'est pas un numéro de cas valide!\" \n", " f\"Numéros de cas disponibles: [1,2,3,4]\")\n", "\n", "step = 0\n", "VectorsList = [A_cols]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Aide\n", "\n", "Pour calculer la factorisation QR de $A$, vous pouvez exécutez la cellule suivant. Si $A$ ne correspond pas aux exigences de la *Définition 1*, supprimez le résultat de cette cellule et insérez **None** comme réponse dans la suivante." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Q1, R1 = al.unique_qr(A)\n", "al.printA(Q1, name=\"Q\")\n", "al.printA(R1, name=\"R\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Aide\n", "\n", "Pour calculer la solution du système linéaire $R\\hat{x}=Q^Tb$, plutôt que de dériver explicitement l'expression pour l'inverse de $R$, utilisez la **méthode de substitution vers l'arrière**!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# INSÉREZ LES VALEURS DE Q,R ET DE LA SOLUTION AU SYSTÈME AU SENS DES MOINS CARREÉS\n", "Q = Q1 # None\n", "R = R1 # None\n", "sol = sp.sets.FiniteSet((0, 0, 0)) # None" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "corrections.Ex2Chapitre9_13_14(Q, R, sol, case_nb=case_number)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exemple 1\n", "\n", "Dans le chapitre 2, tout en considérant la décomposition LU, nous avons signalé qu'un cas d'utilisation très important est celui où plusieurs systèmes linéaires, présentant tous la même matrice de gauche mais différents vecteurs de droite, doivent être résolus.\n", "\n", "De la même manière, la factorisation QR est extrêmement utile lorsque plusieurs systèmes linéaires doivent être résolus au sens des moindres carrés. En effet, une fois la factorisation effectuée, il suffit de résoudre un seul système linéaire triangulaire supérieur (i.e. $R\\hat{x} = Q^Tb$) pour dériver la solution souhaitée.\n", "\n", "Pour des systèmes suffisamment grands, cette approche s'avère plus rapide à la fois que la solution directe des différents systèmes (c'est-à-dire résoudre à chaque étape $A^TAx = A^Tb$, avec $A^TA$ étant calculé une fois pour toutes au début) et que l'application de la décomposition LU (dans ce cas, utilisée avec pivotement partiel pour gagner en performances supplémentaires!) à la matrice $A^TA$.\n", "\n", "**Exécutez la cellule suivante et évaluez les différences de performances ... cela peut prendre quelques secondes**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "M, N = 1000, 500 # dimension des systèmes linéaires à résoudre au sens des moindres carrés\n", "Nt = 10000 # nombre de systèmes linéaires à résoudre au sens des moindres carrés\n", "A = np.random.rand(M, N);\n", "A_ls = A.T @ A\n", "\n", "print(\"----------------------------Factorisation QR----------------------------------\")\n", "start = time.time()\n", - "Q, R = al.unique_qr(A)\n", + "Q, R = np.linalg.qr(A, mode=\"reduced\")\n", "time_qr = time.time() - start\n", "print(\"Temps d'exécution: % f s\" %(time_qr))\n", "\n", - "print(\"----------------------------Décomposition LU--------------------------------\")\n", + "print(\"----------------------------Factorisation LU--------------------------------\")\n", "start = time.time()\n", "P, L, U = lu(A_ls)\n", "time_lu = time.time() - start\n", "print(\"Temps d'exécution: % f s\" %(time_lu))\n", "\n", "print(\"\\n-------Résolution du systèmes linéaires au sens des moindres carreés--------\")\n", "\n", "# résoudre sans utiliser les factorisations QR et LU\n", "start = time.time()\n", "for cnt in range(Nt):\n", " b = np.random.rand(M,1)\n", " x = np.linalg.solve(A_ls, A.T@b)\n", "time_no_fact = time.time() - start\n", "print(\"Sans factorisation: temps d'exécution: % f s\" %time_no_fact)\n", "\n", "# résoudre en utilisant la factorisation LU\n", "start = time.time()\n", "for cnt in range(Nt):\n", " b = np.random.rand(M,1)\n", " y = solve_triangular(L, P@(A.T@b))\n", " x = solve_triangular(U, y)\n", "time_lu += time.time() - start\n", "\n", - "print(\"Avec décomposition LU: temps d'exécution: % f s\" % time_lu)\n", + "print(\"Avec factorisation LU: temps d'exécution: % f s\" % time_lu)\n", "\n", "# résoudre en utilisant la factorisation QR\n", "start = time.time()\n", "for cnt in range(Nt):\n", " b = np.random.rand(M,1)\n", " x = solve_triangular(R, Q.T@b)\n", "time_qr += time.time() - start\n", "\n", "print(\"Avec factorisation QR: temps d'exécution: % f s\" %time_qr)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Vous pouvez comparer les temps d'exécution pour différentes tailles de matrice (c'est-à-dire changer les paramètres M et N) et pour un nombre différent de systèmes lineaires (c'est-à-dire changer le paramètre N_t)**" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.7.4" } }, "nbformat": 4, "nbformat_minor": 4 }