basis_pursuit_subgrad.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "unique-opinion",
   "metadata": {},
   "source": [
    "# Subgradient projection method for basis pursuit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "floating-camera",
   "metadata": {},
   "source": [
    "Consider the basis pursuit problem:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\text{minimize } &&\\|x\\|_1 = \\sum_{i=1}^n |x_i|,\\\\\n",
    "&\\text{subject to } && Ax=b.\n",
    "\\end{aligned}\n",
    "$$\n",
    "We assume that $A\\in\\mathbb{R}^{m\\times n}$ has full row rank, $m<n$, $b\\in \\text{range }A \\subseteq \\mathbb{R}^m$.  A code snippet below generates a random instance of such a problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "nearby-daniel",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def random_instance(m,n):\n",
    "    A = np.random.rand(m,n)\n",
    "    b = A.sum(axis=1)\n",
    "    return (A,b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "academic-albuquerque",
   "metadata": {},
   "source": [
    "We now provide a prototypical implementation of the subgradient projection algorithm for this problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bound-lightning",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import cho_factor, cho_solve\n",
    "\n",
    "def subgradient(x):\n",
    "    # subgradient of 1-norm\n",
    "    return np.sign(x)\n",
    "\n",
    "def project(y,A,b,AAT_fact):\n",
    "    # Projection of y onto the affine manifold {x: Ax=b}\n",
    "    return y-A.T@(cho_solve(AAT_fact, A@y-b))\n",
    "\n",
    "def subgradient_method(x0,A,b):\n",
    "    xk       = x0.copy()\n",
    "    x_best   = xk.copy()\n",
    "    f_best   = float('inf')\n",
    "    f_history= []\n",
    "    # we prefactorize AA^T, since we are solving this system\n",
    "    # at every iteration\n",
    "    AAT_fact = cho_factor(A@A.T)\n",
    "    k_max    = 1000000\n",
    "    k_print  = 10000\n",
    "    x_diff   = 1.0E-06\n",
    "    for k in range(k_max):\n",
    "        # get a subgradient\n",
    "        gk = subgradient(xk)\n",
    "        # Since subradients are not necessarily descent directions,\n",
    "        # we cannot use linesearch. Instead, we use the rule\n",
    "        # where alpha_k->0 but sum(alpha_k)=infinity.\n",
    "        alphak = 1.0/(1+k)\n",
    "        # Project the step onto the feasible set\n",
    "        xk1 = project(xk-alphak*gk,A,b,AAT_fact)\n",
    "        # Compute the new objective and possibly update\n",
    "        # the best so far\n",
    "        fk1 = np.linalg.norm(xk1,ord=1)\n",
    "        f_history.append(fk1)\n",
    "        if fk1 < f_best:\n",
    "            x_best = xk1.copy()\n",
    "            f_best = fk1\n",
    "        # Print messages every k_print iterations\n",
    "        if k%k_print == 0:\n",
    "            print(\"Iter: %6d, fk: %e\" % (k,fk1))\n",
    "        # We stop when the method stops making progress.\n",
    "        # This is not a good practice; one should instead \n",
    "        # check the optimality conditions in general.\n",
    "        if np.linalg.norm(xk-xk1,np.inf) < x_diff:\n",
    "            break\n",
    "        xk = xk1\n",
    "    return x_best,k,f_history"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "moderate-advancement",
   "metadata": {},
   "source": [
    "Let us test it on some small problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "legitimate-happening",
   "metadata": {},
   "outputs": [],
   "source": [
    "A,b=random_instance(5,20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "municipal-square",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter:      0, fk: 1.903090e+01\n",
      "Iter:  10000, fk: 1.646753e+01\n",
      "Iter:  20000, fk: 1.642056e+01\n",
      "Iter:  30000, fk: 1.639310e+01\n",
      "Iter:  40000, fk: 1.637594e+01\n",
      "Iter:  50000, fk: 1.636577e+01\n",
      "Iter:  60000, fk: 1.635744e+01\n",
      "Iter:  70000, fk: 1.635043e+01\n",
      "Iter:  80000, fk: 1.634433e+01\n",
      "Iter:  90000, fk: 1.633896e+01\n",
      "Iter: 100000, fk: 1.633415e+01\n",
      "Iter: 110000, fk: 1.632980e+01\n",
      "Iter: 120000, fk: 1.632583e+01\n",
      "Iter: 130000, fk: 1.632218e+01\n",
      "Iter: 140000, fk: 1.631880e+01\n",
      "Iter: 150000, fk: 1.631566e+01\n",
      "Iter: 160000, fk: 1.631271e+01\n",
      "Iter: 170000, fk: 1.630995e+01\n",
      "Iter: 180000, fk: 1.630733e+01\n",
      "Iter: 190000, fk: 1.630487e+01\n",
      "Iter: 200000, fk: 1.630253e+01\n",
      "Iter: 210000, fk: 1.630031e+01\n",
      "Iter: 220000, fk: 1.629819e+01\n",
      "Iter: 230000, fk: 1.629616e+01\n",
      "Iter: 240000, fk: 1.629422e+01\n",
      "Iter: 250000, fk: 1.629235e+01\n",
      "Iter: 260000, fk: 1.629057e+01\n",
      "Iter: 270000, fk: 1.628885e+01\n",
      "Iter: 280000, fk: 1.628719e+01\n",
      "Iter: 290000, fk: 1.628559e+01\n",
      "Iter: 300000, fk: 1.628404e+01\n",
      "Iter: 310000, fk: 1.628255e+01\n",
      "Iter: 320000, fk: 1.628110e+01\n",
      "Iter: 330000, fk: 1.627970e+01\n",
      "Iter: 340000, fk: 1.627834e+01\n",
      "Iter: 350000, fk: 1.627701e+01\n",
      "Iter: 360000, fk: 1.627573e+01\n",
      "Iter: 370000, fk: 1.627448e+01\n",
      "Iter: 380000, fk: 1.627326e+01\n",
      "Iter: 390000, fk: 1.627208e+01\n",
      "Iter: 400000, fk: 1.627092e+01\n",
      "Iter: 410000, fk: 1.626980e+01\n",
      "Iter: 420000, fk: 1.626870e+01\n",
      "Iter: 430000, fk: 1.626763e+01\n",
      "Iter: 437311, f : 1.626686e+01\n"
     ]
    }
   ],
   "source": [
    "x,k,f_history =subgradient_method(np.zeros(A.shape[1]),A,b)\n",
    "f   = np.linalg.norm(x,ord=1)\n",
    "print(\"Iter: %6d, f : %e\" % (k,f))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "certain-woman",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEGCAYAAAB2EqL0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAZ60lEQVR4nO3df5Dcd33f8ed7b++HLMmyjezDsp3IDsbBdmMVn00CTTgZ4siCCZ0kJVanUwIuamlgCGnaipIxpdPOEIfMAAltrIDrMC2WIdiFgmphTC6G8Mt2sbGMrUgYOVxlLIyM5JN1urvdT//Y793t3n5X+up8t7unez5mbvb7/Xy/3933vRn80vfHfTZSSkiSdDKlThcgSVoaDAxJUiEGhiSpEANDklSIgSFJKqTc6QIW09q1a9P69evndezRo0dZuXLlwha0xNmTfPYln31pthR68tBDDz2bUjo3b9tpHRjr16/nwQcfnNexIyMjDA8PL2xBS5w9yWdf8tmXZkuhJxHxVKttXpKSJBViYEiSCjEwJEmFGBiSpEIMDElSIQaGJKkQA0OSVEjbAiMibouIgxGxu27sqoj4RkQ8GhH/OyLObHHspojYExH7ImLbYtf6p/ft5dEfTy32x0jSktLOM4zbgU1zxj4ObEsp/QPgbuDfzj0oInqAjwE3AJcDWyLi8sUs9L+OfJ/HflJdzI+QpCWnbYGRUrofODRn+DLg/mz5XuA3cw69FtiXUnoypTQB7ADetGiFSpJydXpqkN3ArwOfA/4JcFHOPhcAP6xbHwVe1eoNI2IrsBVgcHCQkZGRUy6qUq0wOTk1r2NPZ2NjY/Ykh33JZ1+aLfWedDow3gZ8NCJuBj4PTOTsEzljLb9XNqW0HdgOMDQ0lOYzb0vPfffQ21vq+jlf2m0pzIPTCfYln31pttR70tHASCk9AVwPEBEvB96Qs9sojWceFwIH2lDd4n+EJC0hHX2sNiLOy15LwB8Cf56z2wPApRFxcUT0ATdSOxtZxLoW890laWlq52O1dwDfAC6LiNGIuInaE09/BzxB7azhv2f7rouInQAppSngncAu4HHg0ymlx9pVtySppm2XpFJKW1ps+kjOvgeAzXXrO4Gdi1SaJKkA/9JbklSIgSFJKsTAaMFnpCSpkYGRw4ekJKmZgSFJKsTAkCQVYmBIkgoxMFrxrrckNTAwcoRzg0hSEwNDklSIgSFJKsTAkCQVYmBIkgoxMFrwISlJamRg5PAZKUlqZmBIkgoxMCRJhRgYkqRCDIwWvOktSY0MjDze9ZakJgaGJKkQA0OSVIiBIUkqxMCQJBViYEiSCjEwcviQlCQ1MzAkSYW0LTAi4raIOBgRu+vGNkTENyPi4Yh4MCKubXHs/oh4dHq/dtUsSZrVzjOM24FNc8ZuAT6QUtoA3Jytt7IxpbQhpTS0OOVJkk6kbYGRUrofODR3GDgzW14DHGhXPSeTnBtEkhpEauN/GSNiPfCFlNKV2forgF3U7jOXgFenlJ7KOe4HwHPUAubWlNL2E3zGVmArwODg4NU7duw45Tp/976jDK1NvPWqVad87OlsbGyMVavsyVz2JZ99abYUerJx48aHWl3JKbe7mDneAbwnpfTZiHgz8Ang9Tn7vSaldCAizgPujYgnsjOWJlmYbAcYGhpKw8PDp1xU+W++RG9vYj7Hns5GRkbsSQ77ks++NFvqPen0U1JvAe7Klj8D5N70TikdyF4PAne32k+StHg6HRgHgNdmy9cBe+fuEBErI2L19DJwPbB77n6SpMXVtktSEXEHMAysjYhR4P3A24GPREQZGCe79xAR64CPp5Q2A4PA3RExXe+nUkr3tKtuSVJN2wIjpbSlxaarc/Y9AGzOlp8ErlrE0nL5kJQkNer0JamuFM4NIklNDAxJUiEGhiSpEANDklSIgSFJKsTAyOE9b0lqZmBIkgoxMCRJhRgYkqRCDAxJUiEGRgtODSJJjQyMHOHcIJLUxMCQJBViYEiSCjEwJEmFGBiteNdbkhoYGDm85S1JzQwMSVIhBoYkqRADQ5JUiIEhSSrEwGjBh6QkqZGBkcOZQSSpmYEhSSrEwJAkFWJgSJIKMTBa8Ka3JDVqW2BExG0RcTAidteNbYiIb0bEwxHxYERc2+LYTRGxJyL2RcS2NlS7+B8hSUtMO88wbgc2zRm7BfhASmkDcHO23iAieoCPATcAlwNbIuLyRa1UktSkbYGRUrofODR3GDgzW14DHMg59FpgX0rpyZTSBLADeNOiFSpJylXu8Of/HrArIj5ELbxenbPPBcAP69ZHgVe1esOI2ApsBRgcHGRkZOSUi5qYmGBysjqvY09nY2Nj9iSHfclnX5ot9Z50OjDeAbwnpfTZiHgz8Ang9XP2ybuh0PKedEppO7AdYGhoKA0PD59yUX1f+zK9vRXmc+zpbGRkxJ7ksC/57Euzpd6TTj8l9Rbgrmz5M9QuP801ClxUt34h+ZeuFpaPSUlSg04HxgHgtdnydcDenH0eAC6NiIsjog+4Efj8Yhbl1CCS1Kxtl6Qi4g5gGFgbEaPA+4G3Ax+JiDIwTnbvISLWAR9PKW1OKU1FxDuBXUAPcFtK6bF21S1JqmlbYKSUtrTYdHXOvgeAzXXrO4Gdi1SaJKmATl+SkiQtEQZGC97zlqRGBkYO73lLUjMDQ5JUiIEhSSrEwJAkFWJgSJIKMTBa8CkpSWpkYORwahBJamZgSJIKedGBERH/fiEKkSR1t1OeSyoiPl2/CmwA/mihCpIkdaf5TD54JKX0L6ZXIuK/LWA9kqQuNZ9LUv9lzvr7FqKQbhJODiJJTU45MFJKP5izfmjhypEkdauTXpKKiJ8p+F4/TSkdeZH1SJK6VJF7GH9J7e/YTnSdJgG3A59cgJokSV3opIGRUto4dywi+lNKxxenJElSN5rv32H8x+mFiHjNwpTSXZJzg0hSg/l+p/euiNgKHAVeAfztwpXUeU4NIknN5vOHe58ADlP7g71vppT+w0IXJUnqPqccGCmlmyJiBfBK4JqIuDWl9C8XvjRJUjcpHBgR8WHgPanmGLXLUKfVpShJUmunctN7DPh8RKwEiIjrI8LAkKRlovAZRkrpDyPinwIjEXGc2g3vbYtWWQd5z1uSmp3KJanXAW+nFhTnAzellPYsVmGSpO5yKpek3gfcnFIaBn4LuDMirluUqiRJXedULkldV7f8aETcAHwWeHWR4yPiNuCNwMGU0pXZ2J3AZdkuZ1Gbj2pDzrH7geeBCjCVUhoqWrckaWG82MkHb6rbfrLJB28H/oy6+aZSSr9d9zl/Qu3vO1rZmFJ69mT1SpIWx4udfHB6/KSTD6aU7o+I9XnbIiKANwNdc4nLmUEkqdG8Jh9cBL8MPJNS2tuqDOBLEZGAW1NK21u9UTZlyVaAwcFBRkZGTrmY48ePM9lfndexp7OxsTF7ksO+5LMvzZZ6T+Y7l9RC2wLccYLtr0kpHYiI84B7I+KJlNL9eTtmYbIdYGhoKA0PD59yMQPf/Aq9vZPM59jT2cjIiD3JYV/y2ZdmS70n852tdsFERBn4DeDOVvuklA5krweBu4Fr21OdJGlaxwMDeD3wREppNG9jRKyMiNXTy8D1wO421idJoo2BERF3AN8ALouI0Yi4Kdt0I3MuR0XEuojYma0OAl+LiEeAbwNfTCnds9j1+n0YktSobfcwUkpbWoz/Ts7YAWBztvwkcNWiFidJOqluuCQlSVoCDAxJUiEGhiSpEANDklSIgSFJKsTAyBF+g5IkNTEwJEmFGBiSpEIMDElSIQZGC84MIkmNDIwc3vSWpGYGhiSpEANDklSIgSFJKsTAkCQVYmC0kHxOSpIaGBg5Ah+TkqS5DAxJUiEGhiSpEANDklSIgdGK97wlqYGBkcOpQSSpmYEhSSrEwJAkFWJgSJIKMTAkSYW0LTAi4raIOBgRu+vG7oyIh7Of/RHxcItjN0XEnojYFxHb2lGvD0lJUqN2nmHcDmyqH0gp/XZKaUNKaQPwWeCuuQdFRA/wMeAG4HJgS0RcvpiF+pCUJDVrW2CklO4HDuVti4gA3gzckbP5WmBfSunJlNIEsAN406IVKknKVe50AZlfBp5JKe3N2XYB8MO69VHgVa3eKCK2AlsBBgcHGRkZOeVijh07xlSpOq9jT2djY2P2JId9yWdfmi31nnRLYGwh/+wC8q8QtbzFkFLaDmwHGBoaSsPDw6dczIoH/ppy+TjzOfZ0NjIyYk9y2Jd89qXZUu9JxwMjIsrAbwBXt9hlFLiobv1C4MBi1yVJatQNj9W+HngipTTaYvsDwKURcXFE9AE3Ap9fzILCuUEkqUk7H6u9A/gGcFlEjEbETdmmG5lzOSoi1kXEToCU0hTwTmAX8Djw6ZTSY+2qW5JU07ZLUimlLS3Gfydn7ACwuW59J7Bz0YqTJJ1UN1ySkiQtAQaGJKkQA6OF5NwgktTAwMjhM1KS1MzAkCQVYmBIkgoxMHJMVRPVThchSV2m41ODdKO/P/QCf9/pIiSpy3iG0cI5A976lqR6nmHkuPz8M+mbOtrpMiSpq3iGkaO3J6j4dxiS1MAzjByPjB72bzEkaQ7PMFrwBEOSGhkYOV778nO5aLWtkaR6/lcxx8r+HqaqnmNIUj0DI8ePnz/O00cNDEmqZ2DkeGD/c50uQZK6joGR4xcvOQeA5BznkjTDwMjxzScPAbUpQiRJNQbGCTz0lJemJGmagZHj3a+7FICzz+jrcCWS1D0MjBzXrK/dw/je00c6XIkkdQ8DI8fLzlsFwB/v2tPhSiSpexgYOV66ZqDTJUhS1zEwTqLiX3xLEmBgnNRH79vb6RIkqSsYGC18dOMZAHzEwJAkoI2BERG3RcTBiNg9Z/xdEbEnIh6LiFtaHLs/Ih6NiIcj4sF21Htm/+w3Ynzhuwfa8ZGS1NXaeYZxO7CpfiAiNgJvAn4hpXQF8KETHL8xpbQhpTS0eCU2+qt/9UsAvPNT33GaEEnLXtsCI6V0P3BozvA7gA+mlI5n+xxsVz1FDGV/jwFw8Xt3drASSeq8aOe/nCNiPfCFlNKV2frDwOeonXmMA3+QUnog57gfAM9R+yK8W1NK20/wGVuBrQCDg4NX79ixY161jo2NsWrVKqop8bZds3NK3b5p5bze73Qw3RM1si/57EuzpdCTjRs3PtTqSk6nA2M38BXg3cA1wJ3AJWlOURGxLqV0ICLOA+4F3pWdsZzQ0NBQevDB+d3yGBkZYXh4GIBnx44z9J+/PLPtkfdfz5oVvfN636WsvieaZV/y2ZdmS6EnEdEyMDr9lNQocFeq+TZQBdbO3SmldCB7PQjcDVzbziLXrurne//p12bWr/rAl7ji5nvaWYIkdVynA+N/AdcBRMTLgT7g2fodImJlRKyeXgauB3bTZmf0ldn/wTfMrB+dqLB+2xf5y6/vb3cpktQR7Xys9g7gG8BlETEaETcBtwGXZJemdgBvSSmliFgXEdN3mQeBr0XEI8C3gS+mlDr2z/v9H3wDX/79X5lZf//nH2P9ti/yVw+NcviFyU6VJUmLrtyuD0opbWmx6Z/l7HsA2JwtPwlctYilnbKXnbea/R98A/ufPcrwh0YA+IPPPEK5FLz6ZWu54cqXcv3lg7xkVX9nC5WkBdS2wDgdrV+7kv0ffAMpJb47epidu5/mnt0/4r13Pcr77n6Unzt3FXsPjvEX/3yI1/38eZRKcfI3laQuZWAsgIjgqovO4qqLzmLbpp/ne08f4Z7dP+JPv7IPgLd/8sG6feFd113KjddcxLqzVnSqZEk6ZQbGAosIrli3hivWreHfXH8Zjz99hBs+8tWZ7SnVJjT86H17WdVfZvDMfgbPHOClZw4wuCZ7nR5bM8C5q/op93T62QRJMjAW3SvOP7Ph6aqDz4/zoV17KEUw0NvDM0fGeebION/6wSGeOTLO1Jzp1CNqj/XWgmSAl67pZ3D1AC9Z1c+KvhJ9PT1cePYKzj9rgLNW9NFXNlwkLQ4Do83OWz3ALb+Vfw+/Wk0cemGCHx0ez4LkOD86Ms4zh8f50ZFxRp97gYeeOsRzJ3gaa2VfD6sHelk1UGb1QJlV/bXX1f05YwO9dctlVvX3snqgzBl9PUR4v0VSIwOji5RKwdpV/axd1c+VF6xpud/xqQqHjk5w6OgEX9/3E1b09ZBS4rkXJjl8bJLnxycZOz7F8+O1n6cPj9fGxqc4OlE5eR1BFiSzgbJqoMyxw+PsOvQoZ/T1cEZfDyv6eljZV2ZFtl77KbOit7at/nWgt4ceb/pLS5qBsQT1l3s4f80Kzl+zgivWtQ6WPJVqysKkFipjWag8Pz2Wrc8GTm2/n4xNcPBIle9/7xlemJji2GSFU51Vpq9cYqBcYkVfLUAGyj0M9PUwUC7V1ntLM+P9M8sl+nt76K9/LZfoz/ZpWC+X6MuW+2aWS5RL4RmTtAAMjGWmpxSsWdE7r7mw6ufBSSkxPlnlhYkpXpiocGyywtHjUxybqMysH5usMD5Z4Vj9+kSF8ckq41O18fGpKuOTFX76wgTHJiscz9brX1/sdGcR1MKkp0RfXbD09gT95R56eyILmB76suXentpPX3Zcb9347Fjt5/ujkxx++P/NrPf2xAmXy9NjpRK95aBcqm0z1NTtDAzNS0TULjf19fCSRfyclBITlSrjk1WOT1U4Plnl+FS2PFVloi5YJrKf2nKFiUqV45PV2mvDtiqTldp71JYTE1NVDh+bZHKqtv9kpZotJyamKkxWEpOVatNDCTN2P/yif9dyKWbDpKd2ZjQdMI3LtYApl2bDZ3p7T2l2ffq4ct179JRK9E6Pzdmnp1S3/8zy7DE9devl6fW6/afXa68ljk0lxicrM+MG4tJnYKirRUR2uakH6PwMwZVqLTgmskCZrCS++rdf55XXXJuFTGKyOrttcjp8Kompai2spqppZmyyUmWqfrk6PZYajpvK2T5VrXJsMs3UNFmpZsu1bTPLlSqT1dp+lVaBt1i+PDuLTymgXJoNptKc17mB05Mt1/+U69ejeVspZ5+mz5kzNvd9StG8bWb/mX2ZXW7xfnn7PnusyoGfHmvcHkGpRNNYBF0XsgaGdApq/wGo3YOZdu4ZJX7u3O7+joNpKSWmqmkmcKZDZTp0KtXa+PTYVDVRqVuf3j7VsH12v+n3rlQTe/bu42cvvphKJVFJaeazqml6n+rM8ZXq7D6196hSqTKzz/Qxk5VaSFbrPrvh2Mrs+vQ+M/umDoXmXH/zlcK7liI/SOrDbfq1VGImvNau7OfT2TeGLiQDQ1pGIiK7lwLQc7LdX5SRqacYHn7Zon7GfM2ESV1IVecETaVuezVN78PMcv32+uCqNmyjYex7jz/OpS+/rGGskpj57Ep9HXXjs+/J7Gdm49OfUz++un9x/tNuYEhadkqloK8Dj3m/5Pl9DF/7M23/3IXinwVLkgoxMCRJhRgYkqRCDAxJUiEGhiSpEANDklSIgSFJKsTAkCQVEunFTgXaxSLix8BT8zx8LfDsApZzOrAn+exLPvvSbCn05GdTSufmbTitA+PFiIgHU0pDna6jm9iTfPYln31pttR74iUpSVIhBoYkqRADo7XtnS6gC9mTfPYln31ptqR74j0MSVIhnmFIkgoxMCRJhRgYc0TEpojYExH7ImJbp+tZCBFxW0QcjIjddWPnRMS9EbE3ez27btt7s99/T0T8Wt341RHxaLbto5F94XBE9EfEndn4tyJifd0xb8k+Y29EvKVNv3IhEXFRRPx1RDweEY9FxLuz8WXbm4gYiIhvR8QjWU8+kI0v257Ui4ieiPhORHwhW19efUkp+ZP9UPvOyu8DlwB9wCPA5Z2uawF+r18BXgnsrhu7BdiWLW8D/ihbvjz7vfuBi7N+9GTbvg38EhDA/wFuyMb/NfDn2fKNwJ3Z8jnAk9nr2dny2Z3uR10PzgdemS2vBv4u+/2XbW+y+ldly73At4BfXM49mdOf3wc+BXxhOf7/qOP/A3TTT/Y/4q669fcC7+10XQv0u62nMTD2AOdny+cDe/J+Z2BX1pfzgSfqxrcAt9bvky2Xqf0la9Tvk227FdjS6V6coEefA37V3szUdAbwf4FX2ZMEcCFwH3Ads4GxrPriJalGFwA/rFsfzcZOR4MppacBstfzsvFWPbggW5473nBMSmkKOAy85ATv1XWy0/9/SO1f1Mu6N9lll4eBg8C9KaVl35PMh4F/B1TrxpZVXwyMRnnfCr/cnjtu1YMT9WY+x3SNiFgFfBb4vZTSkRPtmjN22vUmpVRJKW2g9i/qayPiyhPsvix6EhFvBA6mlB4qekjO2JLvi4HRaBS4qG79QuBAh2pZbM9ExPkA2evBbLxVD0az5bnjDcdERBlYAxw6wXt1jYjopRYW/zOldFc2bG+AlNJPgRFgE/bkNcCvR8R+YAdwXUT8D5ZbXzp9XbCbfqhdN3yS2k2q6ZveV3S6rgX63dbTeA/jj2m8WXdLtnwFjTfrnmT2Zt0D1G6ATt+s25yN/y6NN+s+nS2fA/yA2o26s7Plczrdi7oeBPBJ4MNzxpdtb4BzgbOy5RXAV4E3Luee5PRomNl7GMuqLx1vfrf9AJupPS3zfeB9na5ngX6nO4CngUlq/1q5idq10fuAvdnrOXX7vy/7/feQPcGRjQ8Bu7Ntf8bsTAEDwGeAfdSeALmk7pi3ZeP7gLd2uhdz+vKPqJ3afxd4OPvZvJx7A/wC8J2sJ7uBm7PxZduTnB4NMxsYy6ovTg0iSSrEexiSpEIMDElSIQaGJKkQA0OSVIiBIUkqxMCQ2igi1kfdrMHSUmJgSJIKMTCkDomIS7LvVrim07VIRRgYUgdExGXU5rB6a0rpgU7XIxVR7nQB0jJ0LrXv3vjNlNJjnS5GKsozDKn9DlP7foPXdLoQ6VR4hiG13wTwj4FdETGWUvpUh+uRCjEwpA5IKR3NvpTn3og4mlL6XKdrkk7G2WolSYV4D0OSVIiBIUkqxMCQJBViYEiSCjEwJEmFGBiSpEIMDElSIf8fUhmFCRGrmpkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(f_history)\n",
    "plt.xlabel('k')\n",
    "plt.ylabel(r'$\\|x_k\\|_1$')\n",
    "plt.grid('on')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "flying-commission",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}