{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# L'algorithme de Gram-Schmidt\n",
    "\n",
    "Marc Lorenzi\n",
    "\n",
    "20 juin 2021"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Le théorème de Gram-Schmidt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $E$ un $\\mathbb R$-espace vectoriel de dimension finie $n$ muni d'un produit scalaire $<\\bullet ,\\bullet>$. Soit $\\mathcal B=(e_0,e_1,\\ldots,e_{n-1})$ une base de $E$. On fabrique par récurrence forte sur $k$ une famille $\\mathcal B'=(u_0,u_1,\\ldots,u_{n-1})$ de vecteurs de $E$ en posant :\n",
    "\n",
    "- $u_0=e_0$.\n",
    "- Pour $1\\le k\\le n-1$, \n",
    "\n",
    "$$u_k=e_k-\\sum_{j=0}^{k-1}\\frac{<e_k,u_j>}{<u_j, u_j>}u_j$$\n",
    "\n",
    "__Théorème__ : $\\mathcal B'$ est une base orthogonale de $E$, vérifiant de plus $\\forall k\\in[0,n-1], Vect(u_0,\\ldots,u_k)=Vect(e_0,\\ldots,e_k)$.\n",
    "\n",
    "__Démonstration__ : elle se trouve dans tous les bons cours de mathématiques :-)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. L'algorithme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 La fonction `gram_schmidt`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Écrivons la fonction `gram_schmidt`. Elle prend en paramètre une base $\\mathcal B$ de notre espace vectoriel $E$. Oui, mais c'est qui, $E$ ? Il faut assi passer $E$ en paramètre, d'une façon ou d'une autre. Si on regarde bien les formules de Gram-Schmidt, on voit qu'on a besoin\n",
    "\n",
    "- du produit scalaire\n",
    "- de la soustraction dans $E$\n",
    "- de la multiplication d'un vecteur de $E$ par un réel\n",
    "\n",
    "Nous allons donc représenter $E$ par un triplet `(ps, sub, mul)` de trois fonctions. L'écriture de l'algorithme est alors immédiate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gram_schmidt(B, E):\n",
    "    ps, sub, mul = E\n",
    "    n = len(B)\n",
    "    B1 = [B[0]]\n",
    "    for k in range(1, n):\n",
    "        u = B[k]\n",
    "        for j in range(k):\n",
    "            mu = ps(B[k], B1[j]) / ps(B1[j], B1[j])\n",
    "            u = sub(u, mul(mu, B1[j]))\n",
    "        B1.append(u)\n",
    "    return B1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et voici une fonction qui renvoie `True` lorsque la base $\\mathcal B$ est orthogonale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def est_orthogonale(B, E):\n",
    "    n = len(B)\n",
    "    ps = E[0]\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            if j != i and ps(B[i], B[j]) != 0: return False\n",
    "    return True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Orthonormalisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si l'on veut une base orthonormée, il suffit évidemment d'appliquer Gram-Schmidt puis de diviser les vecteurs de la base obtenue par leur norme.\n",
    "\n",
    "La fonction `normaliser` prend en paramètre un vecteur $u$ et renvoie $u/||u||$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normaliser(u, E):\n",
    "    ps, sub, mul = E\n",
    "    return mul(1 / sqrt(ps(u, u)), u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `orthonormer` prend en paramètre une base orthogonale $\\mathcal B$. Elle renvoie une base orthonormée. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def orthonormer(B, E):\n",
    "    return [normaliser(u, E) for u in B]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Complexité"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous appelons dans ce qui suit « opération » toute opération sur des réels (addition, multiplication, etc.). Notons \n",
    "\n",
    "- $C_p(E)$ le nombre d'opérations nécéssaires pour effectuer le produit scalaire de deux vecteurs de $E$.\n",
    "- $C_s(E)$ le nombre d'opérations nécéssaires pour effectuer le produit d'un vecteur  de $E$ par un réel.\n",
    "- $C_m(E)$ le nombre d'opérations nécéssaires pour effectuer la soustraction de deux vecteurs de $E$.\n",
    "\n",
    "La fonction `gram_schmidt` présente deux boucles imbriquées. La boucle indexée par $k$ effectue $n-1$ itérations, où $n$ est la dimension de $E$. Pour chaque $k$, la boucle indexée par $j$ effectue $k$ itérations. Le nombre total d'itérations est donc\n",
    "$$\\sum_{k=1}^{n-1}k=\\frac 1 2 n(n-1)$$\n",
    "\n",
    "À chaque itération la fonction effectue :\n",
    "\n",
    "- 2 produits scalaires (on pourrait faire mieux en mémorisant certains résultats)\n",
    "- 1 soustraction\n",
    "- 1 produit par un scalaire\n",
    "- 1 division\n",
    "\n",
    "Tout cela ne dépend ni de $j$ ni de $k$. Si l'on note, avec un peu de redondance, $C_n(E)$ la complexité de l'algorithme de Schmidt, on a donc\n",
    "$$C_n(E)=\\frac 1 2 n(n-1) (2C_p(E) + C_s(E) + C_m(E) + 1)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On ne peut pas aller plus loin dans le cas général : tout dépend de l'espace $E$ sur lequel on travaille. Regardons deux exemples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Exemples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 $\\mathbb R^n$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Munissons $\\mathbb R^n$ du produit scalaire canonique. Comme `sympy` sait faire des produits scalaires, soustraire des vecteurs et les multiplier par un réel, il n'y a rien à définir.\n",
    "\n",
    "Voici $\\mathbb R^n$ :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E1 = (\n",
    "    lambda u, v: u.dot(v),\n",
    "    lambda u, v: u - v,\n",
    "    lambda t, u: t * u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Faisons quelques tests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ps1 = E1[0]\n",
    "ps1(Matrix([1,2,3]), Matrix([4,5,6]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sub1 = E1[1]\n",
    "sub1(Matrix([1, 2, 3]), Matrix([2, 1, 4]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Considérons la base $((1, 2, 3), (4, 5,6), (7, 8, 8))$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "B1 = [Matrix([1, 2, 3]), Matrix([4, 5, 6]), Matrix([7, 8, 8])]\n",
    "B1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Orthogonalisons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "B2 = gram_schmidt(B1, E1)\n",
    "B2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La base est-elle bien orthogonale ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "est_orthogonale(B2, E1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tant que nous y sommes, orthonormalisons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "orthonormer(B2, E1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quelle est la complexité de l'algorithme ?\n",
    "\n",
    "- Pour effectuer un produit scalaire dans $\\mathbb R^n$, on fait $n$ multiplications de réels et $n-1$ additions de réels. Donc, $C_p=2n-1$.\n",
    "- $C_s=n$ et $C_m=n$, de façon évidente. \n",
    "\n",
    "Ainsi, \n",
    "\n",
    "$$2C_p+C_s+C_m+1=4n-2+n+n+1=6n-1$$\n",
    "\n",
    "et \n",
    "\n",
    "$$C(n)=\\frac 1 2 n(n-1)(6n-1)$$\n",
    "\n",
    "Ainsi,\n",
    "\n",
    "$$C(n)\\sim 3n^3$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 $\\mathbb R_n[X]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il existe beaucoup de produits scalaires intéressants sur l'espace $\\mathbb R[X]$. Nous allons choisir celui défini par \n",
    "$$<P, Q>=\\int_{-1}^1P(t)Q(t)\\,dt$$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = Symbol('X')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E2 = (\n",
    "    lambda P, Q: integrate(P * Q, (X, -1, 1)),\n",
    "    lambda P, Q: P - Q,\n",
    "    lambda t, P: t * P\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prenons la base canonique de $\\mathbb R_4[X]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "B3 = [1, X, X ** 2, X ** 3, X ** 4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "B4 = gram_schmidt(B3, E2)\n",
    "B4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "est_orthogonale(B4, E2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et orthonormalisons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "orthonormer(B4, E2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(map(simplify, _))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les polynômes que nous obtenons font partie d'une famille importante, ils sont appelés les __polynômes de Legendre__."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quelle est ici la complexité de `gram_schmidt` ? Plaçons nous dans $\\mathbb R_n[X]$. On a $C_s=C_m=n+1$ puisque soustraction et produit par un réel se font par des soustractions ou multiplications coefficient par coefficient.\n",
    "\n",
    "Concernant $C_p$, je n'ai évidemment pas la moindre idée de la complexité de la fonction d'intégration de `sympy`. Cela dit, si l'on code soi-même les calculs d'intégrales, il faut, pour calculer $\\int_{-1}^1 PQ$ :\n",
    "\n",
    "- Multiplier $P$ et $Q$ : cela donne, si on utilise la définition du produit de deux polynômes, $O(n^2)$ opérations élémentaires. On peut calculer précisément cette complexité, mais n'entrons pas dans les détails.\n",
    "- Intégrer terme à terme : $O(n)$ opérations.\n",
    "\n",
    "Ainsi, $C_p=O(n^2)$, donc \n",
    "\n",
    "$$2C_p+C_s+C_m+1=O(n^2)$$\n",
    "\n",
    "et\n",
    "\n",
    "$$C(n)=O(n^4)$$"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}