{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# La décomposition QR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Marc Lorenzi\n",
    "\n",
    "14 février 2020"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Théorème__ : Toute matrice inversible à coefficients réels $A\\in GL_n(\\mathbb R)$ s'écrit de façon unique sous la forme \n",
    "\n",
    "$$A=QR$$\n",
    "\n",
    "où $Q\\in\\mathcal O_n(\\mathbb R)$ est une matrice orthogonale et $R\\in\\mathcal T_n(\\mathbb R)$ est une matrice triangulaire supérieure à  coefficients strictement positifs.\n",
    "\n",
    "Il existe différents algorithmes permettant d'obtenir la décomposition QR : algorithme de Gram-Schmidt, rotations \"élémentaires\" (matrices de Givens), symétries \"élémentaires\" (matrices de Householder). J'ai choisi ici d'utiliser des rotations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import math \n",
    "import random\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (12, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Quelques fonctions auxiliaires"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Facile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous représenterons les matrices en Python par des listes de listes. Si $A$ est une matrice, la liste $A[i]$ est la $i$ème ligne de $A$. Et $A[i][j]$ est le coefficient ligne $i$, colonne $j$ de $A$. Les fonctions `nb_lig` et `nb_col` renvoient le nombre de lignes et le nombre de colonnes de la matrice $A$.\n",
    "\n",
    "__Avertissement__ : Dans toutes nos discussions mathématiques, les indices de ligne ou de colonne des matrices commenceront à 0, afin de pouvoir facilement transposer nos calculs en Python."
   ]
  },
  {
   "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": [
    "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": [
    "print(eye(4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Opérations matricielles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `prodmat` prend en paramètres deux matrices $A$ et $B$. Elle renvoie leur produit $AB$. Pour deux matrices $A\\in\\mathcal M_{pq}(\\mathbb R)$ et $B\\in\\mathcal M_{qr}(\\mathbb R)$, cette fonction effectue $2pqr$ opérations. En particulier, si $p=q=r=n$, le nombre d'opérations est $2n^3$. Retenez ce nombre, il va bientôt réapparaître."
   ]
  },
  {
   "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",
    "            for k in range(q):\n",
    "                C[i][j] += A[i][k] * B[k][j]\n",
    "    return C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = [[1, 2, 3], [4, 5, 6]]\n",
    "B = [[5, 6], [7, 8], [9,10]]\n",
    "print(prodmat(A, B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `submat` prend en paramètres deux matrices $A$ et $B$. Elle renvoie leur différence $A-B$. Elle effectue $pq$ opérations sur des flottants si $A$ et $B$ sont 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": [
    "A = [[1, 2, 3], [4, 5, 6]]\n",
    "B = [[6, 5, 4], [3, 2, 1]]\n",
    "print(submat(A, B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `transp` renvoie la transposée $A^T$ de la matrice $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def transp(A):\n",
    "    p = nb_lig(A)\n",
    "    q = nb_col(A)\n",
    "    B = [p * [0] for i in range(q)]\n",
    "    for i in range(q):\n",
    "        for j in range(p):\n",
    "            B[i][j] = A[j][i]\n",
    "    return B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = [[1, 2], [3, 4], [5, 6]]\n",
    "print(transp(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Affichage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `prnt` affiche proprement une matrice de flottants. Si le choix de 5 chiffres après la virgule ne vous satisfait pas, changez-le."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prnt(A):\n",
    "    p = nb_lig(A)\n",
    "    q = nb_col(A)\n",
    "    for i in range(p):\n",
    "        s = ''\n",
    "        for j in range(q):\n",
    "            s += '%+10.5f ' % A[i][j]\n",
    "        print(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]\n",
    "prnt(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4 Matrices aléatoires"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici enfin une fonction `randmat` qui nous sera bien utile pour nos tests. Elle renvoie une matrice \"aléatoire\" ayant $p$ lignes et $q$ colonnes. Chacun des coefficients de la matrice est un réel choisi uniformément dans l'intervalle $[-r,r]$. Le réel $r$ est par défaut égal à 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": [
    "A = randmat(4, 7, r=3)\n",
    "prnt(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. L'algorithme QR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Quelques rappels sur les matrices orthogonales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $Q\\in\\mathcal M_n(\\mathbb R)$. La matrice $Q$ est dite __orthogonale__ lorsque \n",
    "\n",
    "$$Q^TQ=I_n$$\n",
    "\n",
    "Cette égalité est équivalente à $QQ^T=I_n$. Cela revient encore à dire que $Q$ est inversible et $Q^{-1}=Q^T$. L'ensemble $O_n(\\mathbb R)$ des matrices orthogonales est un groupe multiplicatif, le __groupe orthogonal d'ordre $n$__.\n",
    "\n",
    "Il est facile de voir que les matrices orthogonales ont un déterminant égal à $\\pm 1$. Celles dont le déterminant est égal à 1 sont appelées __matrices de rotation__. Nous noterons $SO_n(\\mathbb R)$ l'ensemble des matrices de rotation. L'ensemble $SO_n(\\mathbb R)$ est un sous-groupe de $O_n(\\mathbb R)$, le __groupe spécial-orthogonal d'ordre $n$__."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Rotations élémentaires"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parmi les éléments de $SO_n(\\mathbb R)$, certains sont plus \"simples\" que les autres en ce sens qu'ils possèdent beaucoup de coefficients nuls. \n",
    "\n",
    "Soit $\\theta\\in\\mathbb R$. Notons $c=\\cos\\theta$ et $s=\\sin\\theta$. Donnons nous $0\\le j < i\\le n-1$. Définissons la matrice $Q_{ij}=Q_{ij}(\\theta)$ comme suit :\n",
    "\n",
    "- Pour tout $k\\ne i,j$, $(Q_{ij})_{kk}=1$.\n",
    "- $(Q_{ij})_{ii}=(Q_{ij})_{jj}=c$.\n",
    "- $(Q_{ij})_{ij}=s$ et $(Q_{ij})_{ji}=-s$.\n",
    "- Tous les autres coefficients de $Q_{ij}$ sont nuls.\n",
    "\n",
    "Les matrices $Q_{ij}(\\theta)$ sont des matrices de rotations, que nous qualifierons d'__élémentaires__. On les appelle aussi __rotations de Givens__. Voici $Q_{ij}(\\theta)$ ci-dessous : le $s$ se trouve ligne $i$ colonne $j$, et le $-s$ se trouve ligne $j$ colonne $i$. Tous les coefficients qui ne sont pas explicités sont nuls."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$Q_{ij}(\\theta)=\\begin{pmatrix}\n",
    "1\\\\\n",
    "&\\ddots\\\\\n",
    "&&1\\\\\n",
    "&&&c&&&&-s\\\\\n",
    "&&&&1&&&\\\\\n",
    "&&&&&\\ddots&&\\\\\n",
    "&&&&&&1&\\\\\n",
    "&&&s&&&&c\\\\\n",
    "&&&&&&&&1\\\\\n",
    "&&&&&&&&&\\ddots\\\\\n",
    "&&&&&&&&&&1\n",
    "\\end{pmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercice__ : Vérifier que $Q_{ij}(\\theta)\\in SO_n(\\mathbb R)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Faire tourner pour annuler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $A\\in\\mathcal M_n(\\mathbb R)$. Soient $0\\le j < i\\le n-1$. Est-il possible de trouver un réel $\\theta$ tel que $A'=Q_{ij}(\\theta)A$ ait son coefficient ligne $i$, colonne $j$ égal à 0 ? La réponse est oui.\n",
    "\n",
    "Par souci de légèreté, notons $Q=Q_{ij}(\\theta)$. Calculons $A'_{kl}$ pour tous $k,l\\in[0,n-1]$. On a\n",
    "\n",
    "$$A'_{kl}=\\sum_{m=0}^{n-1}Q_{km}A_{ml}$$\n",
    "\n",
    "Plusieurs cas se présentent :\n",
    "\n",
    "- Si $k\\ne i,j$, alors $A'_{kl}=A_{kl}$\n",
    "- Si $k=i$, $A'_{il}=Q_{ii}A_{il}+Q_{ij}A_{jl}=cA_{il}+sA_{jl}$.\n",
    "- Si $k=j$, $A'_{jl}=Q_{ji}A_{il}+Q_{jj}A_{jl}=-sA_{il}+cA_{jl}$.\n",
    "\n",
    "En particulier, $A'_{ij}=cA_{ij}+sA_{jj}$. Si nous voulons annuler ce coefficient, il suffit donc de choisir $\\theta$ tel que\n",
    "\n",
    "$$cA_{ij}+sA_{jj}=0$$\n",
    "\n",
    "Posons $x=A_{ij}$ et $y=A_{jj}$.\n",
    "\n",
    "- Si $x=y=0$, $\\theta=0$ fait l'affaire.\n",
    "- Sinon, soit $r=\\sqrt{x^2+y^2}$. Il existe un réel $\\varphi$ tel que $\\cos \\varphi=\\frac y r$ et $\\sin\\varphi =-\\frac x r$. L'équation qui nous intéresse devient \n",
    "\n",
    "$$\\sin(\\varphi-\\theta)=0$$\n",
    "\n",
    "Ainsi, $\\theta=\\varphi$ convient. Remarquons que nous n'avons absolument pas besoin de calculer $\\theta$. Seuls son sinus et son cosinus nous intéressent.\n",
    "\n",
    "La fonction `angle` prend en paramètres deux réels $x$ et $y$. Elle renvoie un couple $(c,s)$ où $c=\\cos\\theta$ et $s=\\sin\\theta$, de sorte que $cx-sy=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def angle(x, y):\n",
    "    if x == 0 and y == 0: return (1, 0)\n",
    "    else:\n",
    "        r = math.sqrt(x ** 2 + y ** 2)\n",
    "        return (y / r, -x / r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = 2, 7\n",
    "c, s = angle(x, y)\n",
    "print(c, s)\n",
    "print(c ** 2 + s ** 2)\n",
    "print(c * x + s * y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il est maintenant aisé, étant donnée une matrice $A$, d'annuler un coefficient non diagonal de $A$ en multipliant $A$ à gauche par une matrice de rotation élémentaire. La fonction `rotation` ci-dessous fait le travail. Elle prend en paramètres la matrice $A$ ainsi que deux indices $i$ et $j$, supposés distincts. Elle prend également en paramètre une matrice $Q$ supposée être une matrice de rotation.\n",
    "\n",
    "La fonction calcule un réel $\\theta$ judicieux, puis effectue le produit $Q_{ij}(\\theta)A$. Le coefficient ligne $i$ colonne $j$ de $A$ est annulé. Elle effectue également le produit $QQ_{ij}(\\theta)^T=QQ_{ij}(-\\theta)$. Comment la matrice $Q$ est-elle modifiée ? Notons $Q'=QQ_{ij}(\\theta)^T$. Calculons $Q'_{kl}$ pour tous $k,l\\in[0,n-1]$. On a\n",
    "\n",
    "$$Q'_{kl}=\\sum_mQ_{km}(Q_{ij}(-\\theta))_{ml}$$\n",
    "\n",
    "Plusieurs cas se présentent :\n",
    "\n",
    "- Si $l\\ne i,j$, alors $Q'_{kl}=Q_{kl}$\n",
    "- Si $l=i$, $Q'_{ki}=cQ_{ki}+sQ_{kj}$.\n",
    "- Si $l=j$, $Q'_{kj}=-sQ_{ki}+cQ_{kj}$.\n",
    "\n",
    "Le calcul de ce produit modifie donc les __colonnes__ $i$ et $j$ de la matrice $Q$. Nous verrons plus loin le pourquoi de cete opération.\n",
    "\n",
    "Remarquons pour terminer que les matrices $A$ et $Q$ sont modifiées sur place."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rotation(A, i, j, Q):\n",
    "    n = len(A)\n",
    "    c, s = angle(A[i][j], A[j][j])\n",
    "    for l in range(n):\n",
    "        Ail, Ajl = A[i][l], A[j][l]\n",
    "        A[i][l] = c * Ail + s * Ajl\n",
    "        A[j][l] = -s * Ail + c * Ajl\n",
    "    for k in range(n):\n",
    "        Qki, Qkj = Q[k][i], Q[k][j]\n",
    "        Q[k][i] = c * Qki + s * Qkj\n",
    "        Q[k][j] = -s * Qki + c * Qkj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = randmat(4,4)\n",
    "Q = eye(4)\n",
    "prnt(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rotation(A, 1, 0, Q)\n",
    "prnt(A)\n",
    "print()\n",
    "prnt(Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercice__ : Réévaluez la cellule ci-dessus en prenant $(i,j)$ successivement égaux à $(2,0),(3,0),(2,1),(3,1)$ et $(3, 2)$. Regardez comment les coefficients de $A$ au-dessous de la diagonale sont l'un après l'autre mis à 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Itérations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $A$ une matrice $n\\times n$. En itérant la fonction `rotation` sur tous les couples $(i, j)$ tels que $0\\le j < i\\le n-1$, les coefficients de $A$ au-dessous de la diagonale sont successivement annulés. \n",
    "\n",
    "Nous n'allons pas faire cela au hasard. En effet, il ne s'agit pas qu'une rotation \"désannule\" un coefficient précédemment annulé ! Nous allons opérer colonne par colonne. On calcule donc $R=Q_{n(n-1)}\\ldots Q_{31}Q_{21}A$, où les $Q_{ij}$ sont des matrices de rotation élémentaires avec des angles judicieux. \n",
    "\n",
    "__Proposition__ : La matrice $R$ est triangulaire. \n",
    "\n",
    "__Démonstration__ : On munit l'ensemble $E=\\{(i,j), 0\\le j < i \\le n-1\\}$ de l'ordre lexicographique :\n",
    "\n",
    "$$(i',j')\\preceq (i,j)\\Leftrightarrow j'<j {\\rm\\ ou\\ }(j'=j {\\rm\\ et\\ }i'\\le i)$$\n",
    "\n",
    "On montre ensuite par induction sur le couple $(i,j)$ que pour tout $(i,j)\\in E$, après l'itération numéro $(i,j)$, on a \n",
    "\n",
    "$$\\forall(i',j')\\preceq(i,j), A_{i'j'}=0$$\n",
    "\n",
    "Récurrence bien posée étant à moitié démontrée, je ne donnerai pas ici les détails de la preuve.\n",
    "\n",
    "__Qu'en est-il de la matrice $Q$ ?__\n",
    "\n",
    "Si, lors de l'appel à `QR`, la matrice $Q$ est égale à $I_n$, la matrice identité d'ordre $n$, à la fin de la dernière itération elle est égale à $Q=Q_{21}^TQ_{31}^T\\ldots Q_{n(n-1)}^T$. On a ainsi $R=Q^TA$. Ainsi, $A=QR$ puisque $Q^{-1}=Q^T$.\n",
    "\n",
    "Ceci explique pourquoi, dans la fonction `QR`, nous avons modifié à chaque itération la matrice $Q$ en la multipliant __à droite__ par $Q_{ij}^T$. Cela évite, en fin de boucle, de transposer la matrice obtenue : nous avons déjà la bonne matrice.\n",
    "\n",
    "__Un dernier petit détail à régler__ : si $A$ est inversible, les coefficients diagonaux de $R$ sont non nuls puisque $R$, tout comme $A$, est inversible. Cependant, certains coefficients diagonaux peuvent être négatifs. Si jamais $R_{ii}<0$, il suffit de changer le signe des coefficients de la ligne $i$ de $R$ et de la colonne $i$ de $Q$. Le produit $QR$ reste égal à $A$ (exercice). De plus, la matrice $Q$ demeure orthogonale (ce n'est en revanche plus forcément une matrice de rotation) et la matrice $R$ reste triangulaire supérieure. La fonction `ajuster_signes` fait cette modification : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ajuster_signes(Q, R):\n",
    "    n = nb_lig(R)\n",
    "    for i in range(n):\n",
    "        if R[i][i] < 0:\n",
    "            for j in range(n):\n",
    "                R[i][j] = - R[i][j]\n",
    "                Q[j][i] = - Q[j][i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercice__ : Soient $Q, R$ deux matrices carrées $n\\times n$. Soit $i\\in[0,n-1]$. On multiplie la colonne $i$ de $Q$ et la ligne $i$ de $R$ par $-1$. Montrer que le produit $QR$ n'est pas modifié."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Avant d'écrire notre fonction `QR` améliorons un peu `rotation`. Comme nous procédons colonne par colonne, les coefficients dans chaque colonne de $A$ s'annulent au fur et à mesure. Lors d'un appel à `rotation` ligne $i$, colonne $j$, inutile de calculer ce qui se passe à gauche de la colonne $j$ ... car il ne s'y passe rien. Voici donc notre nouvelle fonction `rotation2`.\n",
    "\n",
    "Cela n'a l'air de rien mais lors de l'exécution du futur algorithme QR, le nombre d'opérations sur la matrice $A$ est divisé ... par 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rotation2(A, i, j, Q):\n",
    "    n = len(A)\n",
    "    c, s = angle(A[i][j], A[j][j])\n",
    "    for l in range(j, n):\n",
    "        Ail, Ajl = A[i][l], A[j][l]\n",
    "        A[i][l] =  c * Ail + s * Ajl\n",
    "        A[j][l] = -s * Ail + c * Ajl\n",
    "    for k in range(n):\n",
    "        Qki, Qkj = Q[k][i], Q[k][j]\n",
    "        Q[k][i] =  c * Qki + s * Qkj\n",
    "        Q[k][j] = -s * Qki + c * Qkj"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combien `rotation2` effectue-t-elle d'opérations sur des flottants ? Précisément $6(n-j)+6n$. Retenons ce résultat.\n",
    "\n",
    "Nous y voilà. L'algorithme QR est maintenant clair."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def QR(A):\n",
    "    n = nb_lig(A)\n",
    "    R = [A[i].copy() for i in range(n)]\n",
    "    Q = eye(n)\n",
    "    for j in range(n):\n",
    "        for i in range(j + 1, n):\n",
    "            rotation2(R, i, j, Q)\n",
    "    ajuster_signes(Q, R)\n",
    "    return (Q, R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.5 Complexité"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quelle est la complexité $C_n$ de `QR` en termes d'opérations sur des flottants ? Ne comptons pas `ajuster_signes` qui effectue $O(n^2)$ changements de signes. On a\n",
    "\n",
    "$$C_n=\\sum_{j=0}^{n-1}(n-j-1)(6(n-j)+6n)$$\n",
    "\n",
    "Posons $k=n-j-1$.\n",
    "\n",
    "$$\\begin{eqnarray*}\n",
    "C_n&=&\\sum_{k=0}^{n-1}k(6k+6n+6)\\\\\n",
    "&=&6\\sum_{k=0}^{n-1}k^2+6n\\sum_{k=0}^{n-1}k+6n\\\\\n",
    "&=&6\\frac{(n-1)n(2n-1)}{6}+6n\\frac{(n-1)n}2+6n\\\\\n",
    "&\\sim&5n^3\n",
    "\\end{eqnarray*}$$\n",
    "\n",
    "À titre de comparaison, la complexité d'un produit matriciel est $2n^3$. Notre fonction `QR` effectue environ $2.5$ fois plus d'opérations pour décomposer une matrice $A$ que ce qu'il n'en faut pour calculer $A^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lançons `QR` sur une matrice de grande taille. Regardons le temps de calcul et comparons avec le temps de calcul de $A^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = randmat(200, 200)\n",
    "t1 = time.time()\n",
    "Q, R = QR(A)\n",
    "t2 = time.time()\n",
    "print('Temps: %.1fs' % (t2 - t1))\n",
    "t1 = time.time()\n",
    "B = prodmat(A, A)\n",
    "t2 = time.time()\n",
    "print('Temps: %.1fs' % (t2 - t1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La différence entre les deux temps est plus faible que prévu, ce qui est une bonne nouvelle pour QR. __Trop__ bonne nouvelle ? Clairement, compter les opérations sur les flottants n'est pas suffisant. En Python, d'autres opérations sont importantes, je pense en particulier à l'obtention de la valeur d'un coefficient de matrice, ou le changement de la valeur d'un tel coefficient."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tentons une expérience graphique : affichons $\\frac{t'}t$ en fonction de $n$, où $t$ est le temps mis pour calculer le carré d'une matrice $n\\times n$ et $t'$ est le temps mis pour calculer la décomposition QR de cette même matrice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_temps():\n",
    "    s = []\n",
    "    rg = range(20, 101, 2)\n",
    "    for n in rg:\n",
    "        A = randmat(n, n)\n",
    "        t1 = time.time()\n",
    "        Q, R = QR(A)\n",
    "        t2 = time.time()\n",
    "        B = prodmat(A, A)\n",
    "        t3 = time.time()\n",
    "        s.append((t2 - t1) / (t3 - t2))\n",
    "    plt.ylim(1, 2)\n",
    "    plt.grid()\n",
    "    plt.plot(rg, s, 'k')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_temps()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Après quelques égarements, le quotient des deux temps se stabilise autour d'une valeur environ égale à $1.3$. Notre fonction `QR` est donc (avec Python) à peine plus coûteuse en temps que la fonction `prodmat`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.6 Le théorème"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Venons-en au résultat annoncé au début du notebook.\n",
    "\n",
    "__Proposition__ : Soit $A\\in GL_n(\\mathbb R)$ une matrice inversible. Il existe une matrice $Q\\in\\mathcal O_n(\\mathbb R)$ orthogonale et une matrice $R\\in\\mathcal T_n(\\mathbb R)$ triangulaire supérieure dont les coefficients diagonaux sont strictement positifs telles que $A=QR$.\n",
    "\n",
    "__Démonstration__ : Nous venons de la faire. Nous avons même un algorithme pour calculer $Q$ et $R$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Proposition__ : Les matrices $Q$ et $R$ sont définies de façon unique.\n",
    "\n",
    "__Démonstration__ : Supposons $A=QR=Q'R'$. On a alors $Q'^{-1}Q=R'R^{-1}$. La matrice $R''=R'R^{-1}$ est une matrice triangulaire supérieure à coefficients strictement positifs. C'est de plus une matrice orthogonale puisqu'elle est égale à $Q'^{-1}Q$. Il est alors facile de voir que $R''=I_n$. Ainsi, $R'=R$ et $Q'=Q$.\n",
    "\n",
    "__Exercice__ : Montrer qu'une matrice à la fois triangulaire et orthogonale est symétrique et n'a que des $\\pm 1$ sur sa diagonale."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.6 Vérifications"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rappelons-le, la fonction `QR` prend une matrice $A$ en paramètre et renvoie un couple $(Q,R)$ de matrices tel que\n",
    "\n",
    "- $A=QR$\n",
    "- $Q$ est orthogonale.\n",
    "- $R$ est triangulaire supérieure\n",
    "- Les coefficients diagonaux de $R$ sont strictement positifs \n",
    "\n",
    "Ne rêvons pas ! Bien entendu, les calculs en virgule flottante ne donnent que des résultats approchés. La matrice $Q$ est presque orthogonale, $R$ est quasi-triangulaire, et $QR$ est vaguement égal à $A$. __Qu'entendons-nous par \"presque\", \"quasiment\" et \"vaguement\" ?__\n",
    "\n",
    "La fonction `norme2` renvoie la __norme de Frobenius__ d'une matrice $A$ :\n",
    "\n",
    "$$\\Vert A\\Vert=\\left(\\sum_{i,j}A_{ij}^2\\right)^{1/2}$$\n",
    "\n",
    "La fonction $A\\mapsto \\Vert A\\Vert$ est une norme sur l'espace vectoriel $\\mathcal M_n(\\mathbb R)$. La quantité $\\Vert A\\Vert$ est donc un bon indicateur de la \"petitesse\" d'une matrice $A$ et la quantité $\\Vert A-B\\Vert$ est un bon indicateur de la \"proximité\" de deux matrices $A$ et $B$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def norme2(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": "markdown",
   "metadata": {},
   "source": [
    "La fonction `verif_orthog` prend une matrice carrée $Q$ en paramètre et renvoie $\\Vert Q^TQ-I\\Vert$. Elle devrait renvoyer presque 0 lorsque $Q$ est quasi-orthogonale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def verif_orthog(Q):\n",
    "    n = nb_lig(Q)\n",
    "    return norme2(submat(prodmat(Q, transp(Q)), eye(n)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `verif_triang` prend une matrice $R$ en paramètre et renvoie \n",
    "\n",
    "$$\\Gamma(R) = \\left(\\sum_{i>j}R_{ij}^2\\right)^{1/2}$$\n",
    "\n",
    "Cette quantité devrait être quasi-nulle pour une matrice $R$ presque triangulaire supérieure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def verif_triang(R):\n",
    "    n = nb_col(R)\n",
    "    s = 0\n",
    "    for j in range(n):\n",
    "        for i in range(j + 1, n):\n",
    "            s += R[i][j] ** 2\n",
    "    return math.sqrt(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `verif_diag_pos` renvoie `True` si tous les coefficients diagonaux de $R$ sont strictement positifs, et `False` sinon. Ici, c'est tout ou rien. Un coefficient presque positif ne nous satisferait point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def verif_diag_pos(R):\n",
    "    n = nb_col(R)\n",
    "    for i in range(n):\n",
    "        if R[i][i] <= 0: return False\n",
    "    return True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `verif_QR`\n",
    "\n",
    "- Appelle `verif_orthog`\n",
    "- Appelle `verif_triang`\n",
    "- Appelle `verif_diag_pos`\n",
    "- Calcule $\\Vert QR-A\\Vert$\n",
    "\n",
    "et affiche les résultats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def verif_QR(A, Q, R):\n",
    "    print('||Q^TQ-I||: %e' % (verif_orthog(Q)))\n",
    "    print('G(R)      : %e' % (verif_triang(R)))\n",
    "    if verif_diag_pos(R): print('Diag pos  : OK')\n",
    "    else: print('Diag pos : PAS OK')\n",
    "    print('||QR-A||  : %e' % (norme2(submat(prodmat(Q, R), A))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Allons-y."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = randmat(100, 100)\n",
    "Q, R = QR(A)\n",
    "verif_QR(A, Q, R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pas trop mal. Sachant qu'une matrice $100\\times 100$ a $10^4$ coefficients, la moyenne quadratique des erreurs est de l'ordre de $10^{-18}$, c'est à dire moins que la précision de la machine. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.7 Un test avec des matrices de mauvaise réputation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La matrice de Hilbert $(H_{ij})_{0\\le i,j\\le n-1}$ est définie par $H_{ij}=\\frac 1 {i+j+1}$. Les matrices de Hilbert sont connues pour jouer des tours à beaucoup d'algorithmes numériques. Elles sont inversibles mais très mal \"conditionnées\" (je n'entrerai pas dans les détails). Comment se comporte `QR` sur ces matrices ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hilbert(n):\n",
    "    H = [[1 / (i + j + 1) for j in range(n)] for i in range(n)]\n",
    "    return H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prnt(hilbert(10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = hilbert(100)\n",
    "Q, R = QR(H)\n",
    "verif_QR(H, Q, R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous n'avons plus qu'à féliciter notre fonction `QR`. Rien ne lui résiste :-)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 Approximation au sens des moindres carrés"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Généralisation de la décomposition QR à des matrices non carrées"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Considérons un système linéaire $Ax=b$ de $p$ équations à $q$ inconnues, où $p\\ge q$. Supposons que la matrice $A\\in\\mathcal M_{pq}(\\mathbb R)$ est de rang maximal $q$. \n",
    "\n",
    "- Lorsque $p=q$ on a un système de Cramer. Ce système possède une unique solution. \n",
    "- En revanche, lorsque $p>q$, ce système peut fort bien ne pas avoir de solution. Peut-on cependant trouver $x\\in\\mathbb R^q$ tel que $Ax$ soit _le plus près possible_ de $b$ ? Plus précisément, est-il possible de minimiser $||Ax-b||$, où $||u||$ dénote la norme euclidienne du vecteur $u$ ? \n",
    "\n",
    "Complétons la matrice $A$ en une matrice $A'$ de taille $p\\times p$ inversible, en rajoutant $p-q$ colonnes à $A$. Cela est possible car $A$ est de rang $q$. Ses colonnes sont libres dans $\\mathbb R^p$ et peuvent donc être complétées en une base de $\\mathbb R^p$. Écrivons ensuite $A'=QR'$ où $Q\\in\\mathcal O_p(\\mathbb R)$ est orthogonale et $R'\\in\\mathcal T_p(\\mathbb R)$ est triangulaire à coefficients diagonaux strictement positifs.\n",
    "\n",
    "Tronquons les colonnes de $A'$ et $R'$ à partir de la colonne $q+1$. $A'$ redevient $A$ et $R'$ se transforme en une matrice $R$. On a alors $A=QR$ où $Q$ est orthogonale et $R=\\begin{pmatrix}\\hat R\\\\0\\end{pmatrix}$ où la matrice $\\hat R\\in\\mathcal T_q(\\mathbb R)$ est triangulaire à coefficients diagonaux strictement positifs. On a donc le théorème suivant.\n",
    "\n",
    "__Proposition__ : Soient $1\\le q\\le p$ deux entiers. Soit $A\\in\\mathcal M_{pq}(\\mathbb R)$ une matrice que rang $q$. Il existe une matrice orthogonale $Q\\in\\mathcal O_p(\\mathbb R)$ et une matrice triangulaire supérieure à coefficients strictement positifs $\\hat R\\in\\mathcal T_q(\\mathbb R)$ telles que\n",
    "\n",
    "$$A=Q\\begin{pmatrix}\\hat R\\\\0\\end{pmatrix}$$\n",
    "\n",
    "où $0$ désigne la matrice nulle de taille $(p-q)\\times q$.\n",
    "\n",
    "__Remarque__ : Si $p>q$, il n'y a plus unicité : la matrice $A$ peut-être complétée d'une infinité de façons en une matrice inversible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 L'algorithme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `QR1` utilise la fonction `QR`. Il suffit de\n",
    "- Rajouter $p-q$ colonnes à $A$ (on prend ici des coefficients aléatoires) pour obtenir une matrice carrée $A_1$.\n",
    "- Appliquer l'algorithme QR à la matrice $A_1$.\n",
    "- Tronquer la matrice $R$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def QR1(A):\n",
    "    p = nb_lig(A)\n",
    "    q = nb_col(A)\n",
    "    A1 = [A[i].copy() for i in range(p)]\n",
    "    for i in range(p):\n",
    "        A1[i] = A1[i] + [random.uniform(-1, 1) for k in range(p - q)]\n",
    "    Q, R = QR(A1)\n",
    "    for i in range(p): R[i] = R[i][:q]\n",
    "    return (Q, R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prenons un exemple."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = randmat(8, 6)\n",
    "prnt(A)\n",
    "Q, R = QR1(A)\n",
    "print()\n",
    "prnt(Q)\n",
    "print()\n",
    "prnt(R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le concepteur de `verif_QR` a pensé à tout, cette fonction marche encore pour des matrices non carrées."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "verif_QR(A, Q, R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Application à l'approximation au sens des moindres carrés"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dans ce paragraphe, nous confondons allègrement matrices à $q$ lignes et une colonne avec des vecteurs de $\\mathbb R^q$.\n",
    "\n",
    "Revenons à notre problème. Pour tout $x\\in\\mathbb R^q$, on a\n",
    "\n",
    "$$\\Vert Ax-b\\Vert=\\Vert QRx-b\\Vert=\\Vert Q(Rx-Q^Tb)\\Vert$$\n",
    "\n",
    "Comme $Q$ est orthogonale, les produits par $Q$ conservent la norme euclidienne. Ainsi,\n",
    "\n",
    "$$\\Vert Ax-b\\Vert=\\Vert Rx-Q^Tb\\Vert$$\n",
    "\n",
    "Posons $Q^Tb=\\begin{pmatrix}c\\\\d\\end{pmatrix}$ où $c$ possède $q$ lignes et $d$ possède $p-q$ lignes. On a alors\n",
    "\n",
    "$$\\Vert Ax-b\\Vert^2=\\left\\Vert\\begin{pmatrix}\\hat R\\\\0\\end{pmatrix}x-\\begin{pmatrix}c\\\\d\\end{pmatrix}\\right\\Vert^2=\\Vert\\hat Rx-c\\Vert^2+\\Vert d\\Vert^2\\ge\\Vert d\\Vert^2$$\n",
    "\n",
    "Et on a égalité si et seulement si $\\hat Rx=c$. La matrice $\\hat R$ étant inversible, il existe un unique $x\\in\\mathbb R^q$ vérifiant l'égalité.\n",
    "\n",
    "La fonction `moindres_carres` prend en paramètres une matrice $A\\in\\mathcal M_{pq}(\\mathbb R)$ de rang $q$ telle que $1\\le q\\le p$ et une matrice colonne $b\\in\\mathcal M_{p1}(\\mathbb R)$. Elle renvoie l'unique $x\\in\\mathcal M_{q1}(\\mathbb R)$ minimisant $\\Vert Ax-b\\Vert^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def moindres_carres(A, b):\n",
    "    p = nb_lig(A)\n",
    "    q = nb_col(A)\n",
    "    Q, R = QR1(A)\n",
    "    b1 = prodmat(transp(Q), b)\n",
    "    R = R[:q]\n",
    "    b1 = b1[:q]\n",
    "    x = solve_triang_sup(R, b1) \n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il reste à expliquer la ligne `solve_triang_sup(R, b1)` ci-dessus. Elle résout le système $Rx=b_1$. Cela dit, résoudre un système dont la matrice est triangulaire n'est pas bien difficile ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4 Résolution d'un système triangulaire"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit à résoudre le système $Rx=b$ où $R\\in\\mathcal T_n(\\mathbb R)$ est triangulaire supérieure à coefficients diagonaux non nuls. Ce système est de Cramer et possède donc une unique solution. On a pour tout $i\\in[0,n-1]$, \n",
    "\n",
    "$$\\sum_{j=0}^{n-1}R_{ij}x_j=b_i$$\n",
    "\n",
    "On en déduit donc \n",
    "\n",
    "$$x_{n-1}=\\frac 1 {R_{(n-1)(n-1)}}$$\n",
    "\n",
    "et pour $i=n-2,n-3,\\ldots, 0$,\n",
    "\n",
    "$$R_{ii}x_i + \\sum_{j=i+1}^{n-1}R_{ij}x_j=b_i$$\n",
    "\n",
    "et donc\n",
    "\n",
    "$$x_i=\\frac 1 {R_{ii}}\\left(b_i-\\sum_{j=i+1}^{n-1} R_{ij}x_j\\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_triang_sup(R, b):\n",
    "    n = nb_lig(R)\n",
    "    x = [b[i] for i in range(n)]\n",
    "    for i in range(n - 1, -1, -1):\n",
    "        for j in range(i + 1, n):\n",
    "            x[i][0] = x[i][0] - R[i][j] * x[j][0]\n",
    "        x[i][0] = x[i][0] / R[i][i]\n",
    "    return x "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.5 Quelques tests"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testons tout d'abord la qualité de l'approximation sur la résolution d'un système \"aléatoire\" de $n+1$ équations à $n$ inconnues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def norme_euclidienne(u):\n",
    "    s = 0\n",
    "    for i in range(len(u)):\n",
    "        s += u[i][0] ** 2\n",
    "    return math.sqrt(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "norme_euclidienne([[1], [2], [3]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dans la cellule ci-dessous on minimise $\\Vert Ax-b\\Vert$ où $A$ possède $n+1$ lignes et $n$ colonnes. On renvoie $\\Vert Ax-b\\Vert$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 50\n",
    "A = randmat(n + 1, n)\n",
    "b = randmat(n + 1, 1)\n",
    "x = moindres_carres(A, b)\n",
    "E = submat(prodmat(A, x), b)\n",
    "print(norme_euclidienne(E))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Second test, on se donne une matrice $A$ carrée de taille $n$ et une matrice $b\\in\\mathcal M_{n1}(\\mathbb R)$. Pour $k=1,\\ldots, n$, on minimise successivement les quantités $\\Vert A_kx-b\\Vert$ où $A-k$ est la matrice formée des $k$ premières colonnes de $A$. Puis on trace le graphe de l'erreur commise en fonction de $k$. Cette erreur devrait décroître, et devenir nulle lorsque $k=n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = []\n",
    "n = 50\n",
    "A = randmat(n, n)\n",
    "b = randmat(n, 1)\n",
    "for k in range(1, n + 1):\n",
    "    A1 = [[A[i][j] for j in range(k)] for i in range(n)]\n",
    "    x = moindres_carres(A1, b)\n",
    "    t = norme_euclidienne(submat(prodmat(A1, x), b))\n",
    "    s.append(t)\n",
    "plt.grid()\n",
    "plt.plot(range(1, n + 1), s, 'k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4 Approximation polynomiale au sens des moindres carrés"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Passons pour finir à une application intéressante de ce que nous venons d'étudier : l'approximation polynomiale."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 Approximation versus interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soient $x_0<x_1<\\ldots < x_{p-1}$ $p$ réels distincts. Soient $y_0,\\ldots,y_{p-1}$ $p$ réels. Soit $q\\le p$. Existe-t-il un polynôme $P\\in\\mathbb R_{q-1}[X]$ tel que $P(x_0)=y_0,\\ldots,P(x_{p-1})=y_{p-1}$ ?\n",
    "\n",
    "- Si $p=q$, un tel polynôme existe et est unique. Il s'agit du _polynôme d'interpolation de Lagrange_.  \n",
    "\n",
    "- Si, en revanche, $p>q$ l'existence d'un tel polynôme n'est plus assurée. On peut néanmoins chercher un polynôme $P$ tel que la quantité\n",
    "\n",
    "$$\\sum_{k=0}^{p-1} (P(x_k)-y_k)^2$$\n",
    "\n",
    "soit minimale. Au lieu de faire de __l'interpolation__, nous allons faire de __l'approximation__. En fait, nous savons déjà résoudre ce problème. Posons $P=\\sum_{k=0}^{q-1}a_kX^k$. On a pour tout $i\\in[0,p-1]$, \n",
    "\n",
    "$$P(x_i)= \\sum_{k=0}^{q-1}a_kx_i^k$$\n",
    "\n",
    "ce qui peut encore s'écrire matriciellement \n",
    "\n",
    "$$\\begin{pmatrix}P(x_0)\\\\ \\vdots\\\\P(x_{p-1})\\end{pmatrix}=A\\begin{pmatrix}y_0\\\\ \\vdots\\\\y_{p-1}\\end{pmatrix}$$\n",
    "\n",
    "où\n",
    "\n",
    "$$A=\\begin{pmatrix}\n",
    "x_0^0&x_0^1&\\ldots&x_0^{q-1}\\\\\n",
    "x_1^0&x_1^1&\\ldots&x_1^{q-1}\\\\\n",
    "\\vdots&\\vdots&&\\vdots\\\\\n",
    "x_{p-1}^0&x_{p-1}^1&\\ldots&x_{p-1}^{q-1}\\\\\n",
    "\\end{pmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La matrice $A$ est une __matrice de Vandermonde__ de taille $p\\times q$ et de rang $q$ car les $x_i$ sont distincts deux à deux. On peut donc en faire la décomposition $A=QR$ pour obtenir $x=\\begin{pmatrix}a_0\\\\ \\vdots\\\\a_{q-1}\\end{pmatrix}$ tel que $||Ax-b||$ soit minimal, où $b=\\begin{pmatrix}y_0\\\\ \\vdots\\\\y_{p-1}\\end{pmatrix}$. Mais précisément,\n",
    "\n",
    "$$||Ax-b||^2=\\sum_{k=0}^{p-1} (P(x_k)-y_k)^2$$\n",
    "\n",
    "et les coordonnées de $x$ sont les coefficients d'un polynôme $P$ réalisant le minimum cherché.\n",
    "\n",
    "Créons tout d'abord une fonction `vandermonde`. Celle-ci prend en paramètres une liste $xs=[x_0,\\ldots,x_{p-1}]$ de $p$ réels distincts et un entier $q$. Elle renvoie la matrice de Vandermonde de taille $p\\times q$ définie ci-dessus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vandermonde(s, q):\n",
    "    p = len(s)\n",
    "    A = [[s[i] ** j for j in range(q)] for i in range(p)]\n",
    "    return A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prnt(vandermonde([1, 2, 3, 4, 5], 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `approximer` ci-dessous se passe de commentaires. Elle prend en paramètres une liste $xs=[x_0,\\ldots,x_{p-1}]$ de $p$ réels distincts, une liste $b=[y_0,\\ldots,y_{p-1}]$ de $p$ réels et un entier $q$. Elle renvoie la liste des coefficients d'un polynôme $P$ de degré inférieur ou égal à $q-1$ tel que $\\sum_{k=0}^{p-1} (P(x_k)-y_k)^2$ soit minimal. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def approximer(xs, b, q):\n",
    "    A = vandermonde(xs, q)\n",
    "    b = [[x] for x in b]\n",
    "    P = moindres_carres(A, b)\n",
    "    P = [x for [x] in P]\n",
    "    return P"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 Tests"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `eval_poly` prend en paramètre la liste des coefficients d'un polynôme $P$ et un réel $x$. Elle renvoie $P(x)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def eval_poly(P, x):\n",
    "    y, u = 0, 1\n",
    "    for a in P:\n",
    "        y += a * u\n",
    "        u *= x\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prenons $1=1+2X+3X^2$. Que vaut $P(2)$ ? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = [-1, 7, 2]\n",
    "eval_poly(P, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prenons une liste de points situés sur une sinusoïde déformée aléatoirement. Cherchons ensuite un polynôme d'approximation et traçons le tout. En noir, les points sur la sinusoïde déformée. En rouge, le polynôme d'approximation.\n",
    "\n",
    "On a pris sur l'exemple 30 points approximés par un polynôme de degré $<10$. Changez à votre guise ces valeurs. Voyez en particulier ce qui se passe si le degré du polynôme est trop petit, s'il est trop grand. Changez aussi la fonction, les perturbations apportées à celle-ci, etc. Bref, testez.\n",
    "\n",
    "__Avertissement__ : Pour un nombre de points trop élevé, le polynôme d'approximation possède des coefficients extrêmement grands et aussi des coefficients extrêmement petits. En conséquence, les erreurs d'arrondi deviennent prépondérantes dans la fonction `eval_poly` et le résultat renvoyé n'a plus de sens. Par exemple, en prenant $p=q=100$ ci-dessous on devrait obtenir le polynôme d'interpolation : la courbe rouge _devrait_ passer par les points noirs. Je vous laisse tester, je détecte des problèmes dès que $p$ dépasse environ 35 ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_approx_poly(p, q):\n",
    "    a = -3\n",
    "    b = 3\n",
    "    d = (b - a) / p\n",
    "    xs = [a + k * d for k in range(p)]\n",
    "    f = lambda x: math.sin(x)\n",
    "    B = [f(x) for x in xs]\n",
    "    C = [t + 0.2 * random.uniform(-1, 1) for t in B]\n",
    "    P = approximer(xs, C, q)\n",
    "\n",
    "    n1 = 10 * p\n",
    "    d1 = (b - a) / float(n1)\n",
    "    xs1 = [a + k * d1 for k in range(n1)]\n",
    "    ys = [eval_poly(P, t) for t in xs1]\n",
    "    plt.grid()\n",
    "    plt.ylim(-2, 2)\n",
    "    plt.plot(xs, C, '.k', markersize=8)\n",
    "    plt.plot(xs1, ys, 'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_approx_poly(30, 10)"
   ]
  }
 ],
 "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.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}