{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": 203,
      "metadata": {
        "id": "VwavQQWX2Fhk"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import time\n",
        "import scipy\n",
        "from scipy.linalg import circulant\n",
        "import scipy.io\n",
        "import matplotlib.pyplot as plt"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "On veut résoudre:\n",
        "$$\\min_x \\frac{1}{2} ||Ax-b||^2_2$$\n",
        "Le problème peut s'écrire comme:\n",
        "$$\\min_x (a_{11}x_1+a_{12}x_2-b_1)^2 + ...$$\n",
        "$$\\frac{\\delta f}{\\delta x_2}=0 \\Leftrightarrow 2a_{12}(a_{11}x_1+a_{12}x_2-b_1) + ...$$\n",
        "$$(a_{12}^2+a_{22}^2+a_{32}^2)x_2 = -(A(:,2)^TA(:,1))x_1 + ...$$\n",
        "$$(A(:,2)^TA(:,2))x_2 = -(A^TAx)_2 + (A^Tb)_2 + A(:,2)^TA(:,2)x_2$$\n",
        "$$x_2=x_2-\\frac{(A^TAx-A^Tb)_2}{(A^TA)_{2,2}}$$\n",
        "$$x_i^+ = x_i-\\frac{[\\Delta f(x)]_i}{[A^TA]_{i,i}}$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 204,
      "metadata": {
        "id": "ch8nmbMD2Psl"
      },
      "outputs": [],
      "source": [
        "def CDLS(A,b,x0,maxiter=10**4,maxtime=2):\n",
        "  t0 = time.time()\n",
        "\n",
        "  t = np.zeros(maxiter)\n",
        "  fobj = np.zeros(maxiter)\n",
        "\n",
        "  x = x0.copy()\n",
        "\n",
        "  # O(mr^2)\n",
        "  AtA = A.T@A\n",
        "  # O(mr)\n",
        "  Atb = A.T@b\n",
        "  # O(m)\n",
        "  btb = b@b / 2\n",
        "\n",
        "  g = AtA@x-Atb\n",
        "\n",
        "  fobj[0] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2\n",
        "  k = 0\n",
        "  while k < maxiter-1 and t[k] < maxtime:\n",
        "    # On mets à jour le vecteur x entrée par entrée\n",
        "    for i in range(r):\n",
        "      xinew = x[i] - g[i] / AtA[i,i]\n",
        "      g = g - AtA[:,i]*x[i]+AtA[:,i]*xinew # O(r)\n",
        "      x[i] = xinew\n",
        "    fobj[k+1] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2\n",
        "    t[k+1] = time.time() - t0\n",
        "    k += 1\n",
        "  return t[:k+2], fobj[:k+2], x"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 205,
      "metadata": {
        "id": "mx47170C2T1N"
      },
      "outputs": [],
      "source": [
        "def GRLS(A,b,x0,maxiter=10**4,maxtime=2):\n",
        "  t0 = time.time()\n",
        "\n",
        "  t = np.zeros(maxiter)\n",
        "  fobj = np.zeros(maxiter)\n",
        "\n",
        "  x = x0.copy()\n",
        "\n",
        "  # O(mr^2)\n",
        "  AtA = A.T@A\n",
        "  # O(mr)\n",
        "  Atb = A.T@b\n",
        "  # O(m)\n",
        "  btb = b@b / 2\n",
        "\n",
        "  g = AtA@x-Atb\n",
        "\n",
        "  L = max(np.linalg.eigvals(AtA))\n",
        "\n",
        "  alpha = 1/L\n",
        "\n",
        "  fobj[0] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2\n",
        "  k = 0\n",
        "  while k < maxiter-1 and t[k] < maxtime:\n",
        "    d = -g\n",
        "    x = x + alpha*d\n",
        "    g = AtA@x - Atb\n",
        "    fobj[k+1] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2\n",
        "    t[k+1] = time.time() - t0\n",
        "    k += 1\n",
        "  return t[:k+1], fobj[:k+1], x"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 206,
      "metadata": {
        "id": "ZYuPRORf2WVl"
      },
      "outputs": [],
      "source": [
        "def GALS(A,b,x0,maxiter=10**4,maxtime=2):\n",
        "  t0 = time.time()\n",
        "\n",
        "  t = np.zeros(maxiter)\n",
        "  fobj = np.zeros(maxiter)\n",
        "\n",
        "  x = x0.copy()\n",
        "\n",
        "  # O(mr^2)\n",
        "  AtA = A.T@A\n",
        "  # O(mr)\n",
        "  Atb = A.T@b\n",
        "  # O(m)\n",
        "  btb = b@b / 2\n",
        "\n",
        "  g = AtA@x-Atb\n",
        "\n",
        "  L = max(np.linalg.eigvals(AtA))\n",
        "\n",
        "  alpha = 1/L\n",
        "\n",
        "  fobj[0] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2\n",
        "  k = 0\n",
        "\n",
        "  xprec = x.copy()\n",
        "  while k < maxiter-1 and t[k] < maxtime:\n",
        "    y = x + ((k-1)/(k+2))*(x-xprec)\n",
        "    xprec = x\n",
        "    \n",
        "    g = AtA@y - Atb\n",
        "    d = -g\n",
        "    x = y + alpha*d\n",
        "\n",
        "    fobj[k+1] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2\n",
        "    t[k+1] = time.time() - t0\n",
        "    k += 1\n",
        "  return t[:k+1], fobj[:k+1], x"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 207,
      "metadata": {
        "id": "wGTAgRsp2b4r"
      },
      "outputs": [],
      "source": [
        "#Pour la Question 3\n",
        "m, r = 1000, 100\n",
        "A  = np.random.rand(m,r)\n",
        "b  = np.random.rand(m)\n",
        "x0 = np.ones(r)\n",
        "maxiter, maxtime = 10**4, 2\n",
        "t_cd, f_cd, x_cd = CDLS(A,b,x0,maxiter=maxiter,maxtime=maxtime)\n",
        "t_gr, f_gr, x_gr = GRLS(A,b,x0,maxiter=maxiter,maxtime=maxtime)\n",
        "t_ga, f_ga, x_ga = GALS(A,b,x0,maxiter=maxiter,maxtime=maxtime)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 208,
      "metadata": {
        "id": "zroLtOO42jpW"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "[]"
            ]
          },
          "execution_count": 208,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "[encoded data, 41124 characters]",
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "image/png": "[encoded data, 30128 characters]",
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "#Pour la Question 3\n",
        "xopt = np.linalg.solve(A.T@A,A.T@b)\n",
        "fopt = 0.5*xopt@((A.T@A)@xopt)-((A.T@b)@xopt)+b@b/2\n",
        "\n",
        "#Affichage première figure\n",
        "plt.figure()\n",
        "plt.semilogy(f_cd-fopt, label=\"Coord. Descent\")\n",
        "plt.semilogy(f_gr-fopt, label=\"Meth.Gradient\")\n",
        "plt.semilogy(f_ga-fopt, label=\"Meth.Gradient ACC\")\n",
        "plt.legend()\n",
        "plt.plot()\n",
        "\n",
        "\n",
        "#Affichage deuxième figure\n",
        "plt.figure()\n",
        "plt.semilogy(t_cd, f_cd-fopt, label=\"Coord. Descent\")\n",
        "plt.semilogy(t_gr, f_gr-fopt, label=\"Meth.Gradient\")\n",
        "plt.semilogy(t_ga, f_ga-fopt, label=\"Meth.Gradient ACC\")\n",
        "plt.legend()\n",
        "plt.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 209,
      "metadata": {
        "id": "T9duI0AS2vKH"
      },
      "outputs": [],
      "source": [
        "def CDNNLS(A,b,x0,maxiter=10**4,maxtime=2):\n",
        "  t0 = time.time()\n",
        "\n",
        "  t = np.zeros(maxiter)\n",
        "  fobj = np.zeros(maxiter)\n",
        "\n",
        "  x = x0.copy()\n",
        "\n",
        "  # O(mr^2)\n",
        "  AtA = A.T@A\n",
        "  # O(mr)\n",
        "  Atb = A.T@b\n",
        "  # O(m)\n",
        "  btb = b@b / 2\n",
        "\n",
        "  g = AtA@x-Atb\n",
        "\n",
        "  fobj[0] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2\n",
        "  k = 0\n",
        "  while k < maxiter-1 and t[k] < maxtime:\n",
        "    # On mets à jour le vecteur x entrée par entrée\n",
        "    for i in range(r):\n",
        "      xinew = max(0, x[i] - g[i] / AtA[i,i])\n",
        "      g = g - AtA[:,i]*x[i]+AtA[:,i]*xinew # O(r)\n",
        "      x[i] = xinew\n",
        "    fobj[k+1] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2\n",
        "    t[k+1] = time.time() - t0\n",
        "    k += 1\n",
        "  return t[:k+2], fobj[:k+2], x"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 210,
      "metadata": {
        "id": "k1Y3ofiJ2zNc"
      },
      "outputs": [],
      "source": [
        "def GRNNLS(A,b,x0,maxiter=10**4,maxtime=2):\n",
        "  t0 = time.time()\n",
        "\n",
        "  t = np.zeros(maxiter)\n",
        "  fobj = np.zeros(maxiter)\n",
        "\n",
        "  x = x0.copy()\n",
        "\n",
        "  # O(mr^2)\n",
        "  AtA = A.T@A\n",
        "  # O(mr)\n",
        "  Atb = A.T@b\n",
        "  # O(m)\n",
        "  btb = b@b / 2\n",
        "\n",
        "  g = AtA@x-Atb\n",
        "\n",
        "  L = max(np.linalg.eigvals(AtA))\n",
        "\n",
        "  alpha = 1/L\n",
        "\n",
        "  fobj[0] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2\n",
        "  k = 0\n",
        "  while k < maxiter-1 and t[k] < maxtime:\n",
        "    d = -g\n",
        "    x = np.maximum(0, x + alpha*d)\n",
        "    g = AtA@x - Atb\n",
        "    fobj[k+1] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2\n",
        "    t[k+1] = time.time() - t0\n",
        "    k += 1\n",
        "  return t[:k+1], fobj[:k+1], x"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 211,
      "metadata": {
        "id": "acjPuR7x2z-7"
      },
      "outputs": [],
      "source": [
        "def GANNLS(A,b,x0,maxiter=10**4,maxtime=2):\n",
        "  t0 = time.time()\n",
        "\n",
        "  t = np.zeros(maxiter)\n",
        "  fobj = np.zeros(maxiter)\n",
        "\n",
        "  x = x0.copy()\n",
        "\n",
        "  # O(mr^2)\n",
        "  AtA = A.T@A\n",
        "  # O(mr)\n",
        "  Atb = A.T@b\n",
        "  # O(m)\n",
        "  btb = b@b / 2\n",
        "\n",
        "  g = AtA@x-Atb\n",
        "\n",
        "  L = max(np.linalg.eigvals(AtA))\n",
        "\n",
        "  alpha = 1/L\n",
        "\n",
        "  fobj[0] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2\n",
        "  k = 0\n",
        "\n",
        "  xprec = x.copy()\n",
        "  while k < maxiter-1 and t[k] < maxtime:\n",
        "    y = x + ((k-1)/(k+2))*(x-xprec)\n",
        "    xprec = x\n",
        "    \n",
        "    g = AtA@y - Atb\n",
        "    d = -g\n",
        "    x = np.maximum(0, y + alpha*d)\n",
        "\n",
        "    fobj[k+1] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2\n",
        "    t[k+1] = time.time() - t0\n",
        "    k += 1\n",
        "  return t[:k+1], fobj[:k+1], x"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 212,
      "metadata": {
        "id": "P_pEc_un23n3"
      },
      "outputs": [],
      "source": [
        "#Pour la Question 4\n",
        "m, r = 10000, 1000\n",
        "A = np.random.rand(m,r)\n",
        "b = np.random.rand(m)\n",
        "\n",
        "x0 = np.ones(r)\n",
        "maxiter, maxtime = 10**4, 2\n",
        "t_cd, f_cd, x_cd = CDNNLS(A,b,x0,maxiter=maxiter,maxtime=maxtime)\n",
        "t_gr, f_gr, x_gr = GRNNLS(A,b,x0,maxiter=maxiter,maxtime=maxtime)\n",
        "t_ga, f_ga, x_ga = GANNLS(A,b,x0,maxiter=maxiter,maxtime=maxtime)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 213,
      "metadata": {
        "id": "A7zXKfw32-C5"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "[]"
            ]
          },
          "execution_count": 213,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "[encoded data, 29312 characters]",
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "image/png": "[encoded data, 41164 characters]",
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "#Pour la Question 4\n",
        "xopt, fopt = scipy.optimize.nnls(A,b)\n",
        "fopt = (fopt**2)/2\n",
        "\n",
        "#Affichage première figure\n",
        "plt.figure()\n",
        "plt.semilogy(f_cd-fopt, label=\"Coord. Descent\")\n",
        "plt.semilogy(f_gr-fopt, label=\"Meth.Gradient\")\n",
        "plt.semilogy(f_ga-fopt, label=\"Meth.Gradient ACC\")\n",
        "plt.legend()\n",
        "plt.plot()\n",
        "\n",
        "\n",
        "#Affichage deuxième figure\n",
        "plt.figure()\n",
        "plt.semilogy(t_cd, f_cd-fopt, label=\"Coord. Descent\")\n",
        "plt.semilogy(t_gr, f_gr-fopt, label=\"Meth.Gradient\")\n",
        "plt.semilogy(t_ga, f_ga-fopt, label=\"Meth.Gradient ACC\")\n",
        "plt.legend()\n",
        "plt.plot()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "$$\\min_{H\\geq0, W\\geq0} ||X-WH||_F^2$$\n",
        "Si on considère $H$ fixe -> convexe \\\n",
        "Si on considère $W$ fixe-> convexe\n",
        "\n",
        "Pour l'algo on alterne entre les deux: \\\n",
        "for i = maxiter \\\n",
        "W <- Update(X,H) \\\n",
        "H <- Update(X,H)\n",
        "\n",
        "$||X-WH||_F^2=||X^T-H^TW^T||_F^2$ \\\n",
        "$H^* \\leftarrow UpdateH(X,W)$ \\\n",
        "$W^* \\leftarrow UpdateH(X^T,W^T)^T$\n",
        "\n",
        "$||X-WH||_F^2 = ||x(:,1)-WH(:,1)||_2^2 + ||x(:,2)-WH(:,2)||_2^2 + \\ldots$\n",
        "\n",
        "Prenons un sous-problème: \\\n",
        "$||x(:,2)-WH(:,2)||_2^2$ \\\n",
        "$\\min_{x\\geq0} ||Ax-b||_2^2$\n",
        "=> NNLS"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 224,
      "metadata": {
        "id": "6jAZ1LdT3FIC"
      },
      "outputs": [],
      "source": [
        "def initialisation(m,n,r):\n",
        "  W = np.random.random((m,r))\n",
        "  H = np.random.random((r,n))\n",
        "  return W, H"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 225,
      "metadata": {
        "id": "UM8Q2o_a3Q4L"
      },
      "outputs": [],
      "source": [
        "def NNLS(AtA,Atb,x,maxiter=10**2):\n",
        "  r = len(x)\n",
        "  g = AtA@x - Atb\n",
        "\n",
        "  for iter in range(maxiter):\n",
        "    for i in range(r):\n",
        "      xinew = np.maximum(0, x[i] - g[i] / AtA[i,i])\n",
        "      g = g - AtA[:,i]*x[i]+AtA[:,i]*xinew # O(r)\n",
        "      x[i] = xinew\n",
        "\n",
        "  return x"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 230,
      "metadata": {
        "id": "3bnnXhgm3TEP"
      },
      "outputs": [],
      "source": [
        "def optimizeH(X,W,H):\n",
        "  # Pour X et W donnés,\n",
        "  # ce code renvoie la matrice H minimisant ||X-WH||_F^T avec H>=0\n",
        "  # en résolvant un NNLS pour chaque colonne de H\n",
        "  n = X.shape[1]\n",
        "  AtA = W.T@W # O(mr^2)\n",
        "  AtB = W.T@X # O(mrn) <- le + cher\n",
        "\n",
        "  for j in range(n):\n",
        "    H[:,j] = NNLS(AtA, AtB[:,j], H[:,j])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 236,
      "metadata": {
        "id": "9lAMyXi23cVA"
      },
      "outputs": [],
      "source": [
        "def alternatealgo_l2nmf(X,r,maxiter):\n",
        "  m, n = X.shape\n",
        "  W, H = initialisation(m,n,r)\n",
        "  fobj = np.zeros(maxiter)\n",
        "\n",
        "  for iter in range(maxiter):\n",
        "    optimizeH(X,W,H) # optimisation de H étant donnés X et W\n",
        "    optimizeH(X.T, H.T, W.T) # optimisation de W étant donnés X^T et W^T\n",
        "\n",
        "  return W, H"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 237,
      "metadata": {
        "id": "Go3S1Mik3kSm"
      },
      "outputs": [],
      "source": [
        "#Pour la Question 5\n",
        "X = np.array([[0 , 1 , 2 , 2 , 1 , 0],\n",
        "              [0 , 0 , 1 , 2 , 2 , 1],\n",
        "              [1 , 0 , 0 , 1 , 2 , 2],\n",
        "              [2 , 1 , 0 , 0 , 1 , 2],\n",
        "              [2 , 2 , 1 , 0 , 0 , 1],\n",
        "              [1 , 2 , 2 , 1 , 0 , 0]])\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 238,
      "metadata": {
        "id": "Ruodf3-13-UT"
      },
      "outputs": [
        {
          "data": {
            "image/png": "[encoded data, 180444 characters]",
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "#Pour la Question 5\n",
        "def load_matrix(name):\n",
        "  S = scipy.io.loadmat(name)\n",
        "  X = S[\"M\"]\n",
        "  return X\n",
        "def monaffichage(X,nbimages,nbrow,nbcol,m_image,n_image):\n",
        "  if nbrow*nbcol==nbimages and m_image*n_image==X.shape[0]:\n",
        "    Xtot = np.zeros((nbrow*m_image,nbcol*n_image))\n",
        "    i, j = 0, 0\n",
        "    for k in range(nbimages):\n",
        "      reshaped_image = np.reshape(X[:,k], (m_image,n_image)) / max(X[:,k])\n",
        "      Xtot[(i)*m_image:(i+1)*m_image,(j)*n_image:(j+1)*n_image] = reshaped_image.T\n",
        "      if j>=nbcol-1:\n",
        "        i+=1\n",
        "        j = 0\n",
        "      else: j+=1\n",
        "    plt.figure()\n",
        "    plt.imshow(1 - Xtot, cmap='gray')\n",
        "\n",
        "\n",
        "nbvisages = 200;\n",
        "X         = load_matrix('exemplenmf2.mat')\n",
        "Xtemp     = X[:,:nbvisages]\n",
        "monaffichage(Xtemp,200,10,20,19,19)\n",
        "\n",
        "\n",
        "#Lancement de l'algorithme sur la matrice Xtemp avec r=5\n",
        "r = 10\n",
        "W, H = alternatealgo_l2nmf(Xtemp,r,maxiter=5)\n",
        "#plt.figure()\n",
        "#plt.semilogy(f)\n",
        "#plt.xlabel(\"Itérations\")\n",
        "#plt.ylabel(\"Erreur ||X-WH||\")\n",
        "#plt.title(\"Evolution de l'erreur\")\n",
        "#plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 239,
      "metadata": {
        "id": "i8zEO66q41mB"
      },
      "outputs": [
        {
          "data": {
            "image/png": "[encoded data, 18780 characters]",
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "#Affichage de W\n",
        "monaffichage(W,10,2,5,19,19)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 240,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "image/png": "[encoded data, 166824 characters]",
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "#Affichage de W\n",
        "monaffichage(W@H,200,10,20,19,19)"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "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": 0
}
