{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8e81b5c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4301d681",
   "metadata": {},
   "source": [
    "# Question 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c3e61f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradientquad(Q,c,p,x0,maxiter=100,choix='exact'):\n",
    "    fobj = np.zeros(maxiter)\n",
    "    x = x0.copy()\n",
    "\n",
    "    valp, vecp = np.linalg.eig(Q)\n",
    "    L = np.max(valp)\n",
    "\n",
    "    # Si le choix est \"constant\", alpha ne dépend pas de l'itéré, on le calcule une seule fois.\n",
    "    alpha = 1/L\n",
    "\n",
    "    for k in range(maxiter):\n",
    "        # Gradient de fonction quadratique\n",
    "        d = -(Q@x-c)\n",
    "        if choix == \"exact\":\n",
    "            # Dépend de l'itéré, on le calcule à chaque itéré\n",
    "            alpha = (d@d)/(d@Q@d)\n",
    "        x = x + alpha*d\n",
    "        fobj[k] = 0.5*x@Q@x - c@x + p\n",
    "        if np.linalg.norm(d) <= 1e-6:\n",
    "            break\n",
    "    return x, fobj, k"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c418d62",
   "metadata": {},
   "source": [
    "## Expérience 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7f35335f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x240ff222660>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 27116 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Q = np.array([[10.0,-6.0],[-6.0,10.0]])\n",
    "c = np.array([4.0,-3.0])\n",
    "p = 0\n",
    "x0 = np.array([1.0,1.0])\n",
    "\n",
    "# Méthode exacte\n",
    "xopt, fobj, nbiter = gradientquad(Q,c,p,x0,maxiter=25,choix='exact')\n",
    "plt.plot(fobj, 'b', label=\"Gradient avec pas exact\")\n",
    "\n",
    "# Méthode constante\n",
    "xopt, fobj, nbiter = gradientquad(Q,c,p,x0,maxiter=25,choix='constant')\n",
    "plt.plot(fobj, 'r', label=\"Gradient avec pas constant\")\n",
    "\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c266cf47",
   "metadata": {},
   "source": [
    "## Expérience 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "90781b5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import hilbert"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "4b09233b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\conta\\AppData\\Local\\Temp\\ipykernel_10208\\2863094032.py:18: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  fobj[k] = 0.5*x@Q@x - c@x + p\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, \"Nombre d'itérations\")"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 37224 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nbiter_exact = []\n",
    "nbiter_constant = []\n",
    "\n",
    "for n in range(3,101):\n",
    "    Q = hilbert(n)\n",
    "    c = np.zeros(n)\n",
    "    p = 0\n",
    "    x0 = np.ones(n)\n",
    "    # Méthode exacte\n",
    "    xopt, fobj, nbiter = gradientquad(Q,c,p,x0,maxiter=10**6,choix='exact')\n",
    "    nbiter_exact.append(nbiter)\n",
    "    # Méthode constante\n",
    "    xopt, fobj, nbiter = gradientquad(Q,c,p,x0,maxiter=10**6,choix='constant')\n",
    "    nbiter_constant.append(nbiter)\n",
    "\n",
    "plt.plot(nbiter_exact, label=\"Pas exact\")\n",
    "plt.plot(nbiter_constant, label=\"Pas constant\")\n",
    "plt.xlabel(\"Dimension n\")\n",
    "plt.ylabel(\"Nombre d'itérations\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c28bbfd",
   "metadata": {},
   "source": [
    "# Question 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "795525ae",
   "metadata": {},
   "source": [
    "## Implémenter les fonctions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e8da3a46",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_grad_conv(x, c=10.0):\n",
    "    f = 2*np.log(np.cosh((c/2)*(x[0]+2*x[1]))) + 2*np.log(np.cosh((c/2)*(2*x[0]-x[1])))\n",
    "    grad = np.array([\n",
    "        c*np.tanh((c/2)*(x[0]+2*x[1])) + 2*c*np.tanh((c/2)*(2*x[0]-x[1])),\n",
    "        2*c*np.tanh((c/2)*(x[0]+2*x[1])) - c*np.tanh((c/2)*(2*x[0]-x[1]))\n",
    "        ])\n",
    "    L = 2.5*c**2\n",
    "    return f, grad, L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "8fececb8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_grad_strong(x, c=10.0, mu=10**2):\n",
    "    f, grad, L = f_grad_conv(x, c)\n",
    "    f += (mu/2)*(x[0]**2 + x[1]**2)\n",
    "    grad += mu*x\n",
    "    L += mu\n",
    "    return f, grad, L"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e106212",
   "metadata": {},
   "source": [
    "## Implémenter la fonction methgradient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "aad0cb25",
   "metadata": {},
   "outputs": [],
   "source": [
    "def methgradient(fct, x0, maxiter=100, tol=1e-9):\n",
    "    x = x0.copy()\n",
    "    fobj = np.zeros(maxiter)\n",
    "\n",
    "    for k in range(maxiter):\n",
    "        f, grad, L = fct(x)\n",
    "        d = grad\n",
    "        alpha = 1/L\n",
    "        x = x - alpha*d\n",
    "        f, grad, L = fct(x)\n",
    "        fobj[k] = f\n",
    "        if np.linalg.norm(grad) <= tol:\n",
    "            break\n",
    "    return x, fobj[:k]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b2a063e",
   "metadata": {},
   "source": [
    "## Expérience"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "5a665f93",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x2409f1e4550>"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 26912 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x0 = np.array([1.0,1.0])\n",
    "\n",
    "xopt, fobj = methgradient(f_grad_conv, x0, maxiter=70)\n",
    "plt.semilogy(fobj, 'b', label=\"Fonction convexe\")\n",
    "\n",
    "xopt, fobj = methgradient(f_grad_strong, x0, maxiter=70)\n",
    "plt.semilogy(fobj, 'r', label=\"Fonction fortement convexe\")\n",
    "\n",
    "plt.xticks(np.arange(0,18,2))\n",
    "plt.legend(loc='upper right')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6e7014d",
   "metadata": {},
   "source": [
    "# Question 3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6adbe45c",
   "metadata": {},
   "source": [
    "## Implémenter la fonction\n",
    "Note: La fonction n'est pas L-smooth parce qu'on peut pas trouver une fonction quadratique avec un seul L tel que cette fonction est au dessus de $x^4$.\n",
    "La courbature de $x^4$ en $\\infty$ est aussi $\\infty$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "ca042253",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fct(x):\n",
    "    f = (x[0]+x[1])**4 - 2*(x[0]+x[1])**2 + 1\n",
    "    g1 = 4*(x[0]+x[1])**3 - 4*(x[0]+x[1])\n",
    "    g = np.array([g1, g1])\n",
    "    return f, g"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e8069ca",
   "metadata": {},
   "source": [
    "## Implémenter les deux conditions de Wolfe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "id": "baf4f554",
   "metadata": {},
   "outputs": [],
   "source": [
    "def w1(fct, x, d, alpha, beta1):\n",
    "    f0, g0 = fct(x)\n",
    "    f1, _ = fct(x+alpha*d)\n",
    "    return f1 <= f0 + beta1*alpha*(g0@d)\n",
    "\n",
    "def w2(fct, x, d, alpha, beta2):\n",
    "    _, g0 = fct(x)\n",
    "    _, g1 = fct(x+alpha*d)\n",
    "    return d@g1 >= beta2*(d@g0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b935acb",
   "metadata": {},
   "source": [
    "## Implémenter la bissection de Wolfe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "id": "253f97be",
   "metadata": {},
   "outputs": [],
   "source": [
    "def wolfebissection(fct, x, d, alpha0=1.0, beta1=1e-4, beta2=0.9):\n",
    "    aleft = 0\n",
    "    aright = np.inf\n",
    "    alpha = alpha0\n",
    "\n",
    "    while True:\n",
    "        if w1(fct, x, d, alpha, beta1) and w2(fct, x, d, alpha, beta2):\n",
    "            return alpha\n",
    "        if not w1(fct, x, d, alpha, beta1):\n",
    "            # On bouge la borne de droite vers la gauche\n",
    "            aright = alpha\n",
    "            alpha = (aleft + aright)/2\n",
    "        elif not w2(fct, x, d, alpha, beta2):\n",
    "            # On bouge la borne de gauche vers la droite\n",
    "            aleft = alpha\n",
    "            if aright < np.inf:\n",
    "                alpha = (aleft + aright)/2\n",
    "            else:\n",
    "                alpha *= 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dc9c161",
   "metadata": {},
   "source": [
    "## Implémenter la méthode du gradient avec recherche de Wolfe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "id": "acd3c8f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def methgradient(fct, x0, maxiter=100, tol=1e-6, alpha0=1e-6):\n",
    "    x = x0.copy()\n",
    "    fobj = np.zeros(maxiter)\n",
    "    alpha = alpha0\n",
    "\n",
    "    for k in range(maxiter):\n",
    "        f, g = fct(x)\n",
    "        d = -g\n",
    "        alpha = wolfebissection(fct, x, d, alpha)\n",
    "        x = x + alpha*d\n",
    "        fobj[k] = f\n",
    "        if np.linalg.norm(g) <= tol:\n",
    "            break\n",
    "    return x, fobj[:k]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6f5ebc5",
   "metadata": {},
   "source": [
    "## Exemple numérique"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "id": "fcef62aa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.50000003 -0.50000003]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 13112 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x0 = np.array([1.0,1.0])\n",
    "maxiter = 30\n",
    "# Importance de alpha0!\n",
    "xopt, fobj = methgradient(fct, x0, maxiter=maxiter, alpha0=0.05)\n",
    "plt.plot(fobj, 'b', label=\"Gradient avec Wolfe\")\n",
    "plt.legend(loc=\"upper right\")\n",
    "print(xopt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c36d5972",
   "metadata": {},
   "source": [
    "## Bonus: Méthode gradient accéléré"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "id": "94113279",
   "metadata": {},
   "outputs": [],
   "source": [
    "def methgradientaccel(fct, x0, maxiter=100, tol=1e-6, alpha0=1e-6):\n",
    "    x = x0.copy()\n",
    "    fobj = np.zeros(maxiter)\n",
    "\n",
    "    xprev = x\n",
    "    alpha = alpha0\n",
    "\n",
    "    for k in range(maxiter):\n",
    "        beta = k/(k+3) # Paul Tseng (version plus simple)\n",
    "        y = x + beta*(x - xprev)\n",
    "        f, g = fct(y)\n",
    "        d = -g\n",
    "        alpha = wolfebissection(fct, y, d, alpha)\n",
    "        xprev = x.copy()\n",
    "        x = y + alpha*d\n",
    "        f, g = fct(x)\n",
    "        fobj[k] = f\n",
    "        if np.linalg.norm(g) <= tol:\n",
    "            break\n",
    "    return x, fobj[:k]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8b110f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = np.array([1.0,1.0])\n",
    "maxiter = 30\n",
    "# Importance de alpha0!\n",
    "xopt, fobj = methgradientaccel(fct, x0, maxiter=maxiter, alpha0=0.05)\n",
    "plt.plot(fobj, 'b', label=\"Gradient avec Wolfe\")\n",
    "plt.legend(loc=\"upper right\")\n",
    "print(xopt)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "optinonlin",
   "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.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
