{ "cells": [ { "cell_type": "markdown", "source": [ "**Matrix completion using trace norm regularization**\n", "\n", "This shows the use of Douglas-Rachford's algorithm to perform matrix completion with trace norm (nuclear norm) regularization." ], "metadata": { "id": "2PL4yuRZK1tO" }, "id": "2PL4yuRZK1tO" }, { "cell_type": "code", "execution_count": null, "id": "obvious-parker", "metadata": { "id": "obvious-parker" }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "silver-facing", "metadata": { "id": "silver-facing" }, "source": [ "We consider $n \\times n$ matrices." ] }, { "cell_type": "code", "execution_count": null, "id": "silver-boost", "metadata": { "id": "silver-boost" }, "outputs": [], "source": [ "n = 100" ] }, { "cell_type": "markdown", "id": "defensive-butler", "metadata": { "id": "defensive-butler" }, "source": [ "Generate a rank $r$ random matrix." ] }, { "cell_type": "code", "execution_count": null, "id": "upset-government", "metadata": { "id": "upset-government" }, "outputs": [], "source": [ "r = 10\n", "x0 = np.dot(np.random.randn(n,r),np.random.randn(r,n))" ] }, { "cell_type": "markdown", "id": "unavailable-knitting", "metadata": { "id": "unavailable-knitting" }, "source": [ "Display the singular values." ] }, { "cell_type": "code", "execution_count": null, "id": "constitutional-realtor", "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 282 }, "id": "constitutional-realtor", "outputId": "7e202837-25d2-41a4-fcc1-22781455b987" }, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "[]" ] }, "metadata": {}, "execution_count": 5 }, { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAY5UlEQVR4nO3dfZBd9X3f8ffn7t29KwmttJIWPXclgZAQJC7s1l7XLiGAG2F7LGbqerA9ttIKNNPQxondMdiettM/PBNPMnbsxqEjA7HckTGUuIU4jmMqyyZJLcLKmCdJoEWOhIRkLUIPPFharfbbP+7R5las0N6nveee+3nNaPY83Lv3ezjw4aff73d/RxGBmZllS67RBZiZWe053M3MMsjhbmaWQQ53M7MMcribmWVQvtEFAMybNy+WLVvW6DLMzJrKjh07XomInonOpSLcly1bxuDgYKPLMDNrKpL2Xeicu2XMzDLI4W5mlkEOdzOzDHK4m5ll0EXDXdJ9ko5IenaCc5+RFJLmJfuS9DVJQ5KelnRtPYo2M7O3N5mW+zeBtecflLQU+JfA/pLDNwMrkz8bgburL9HMzMp10XCPiMeAVyc49RXgs0DpspLrgG9F0XZgtqSFNal0Ajv2HePr24bYse9YvT7CzKwpVTTPXdI64GBEPCWp9NRi4KWS/QPJsUMVV3gBO/Yd42Pf2M7I6BiFfI4ttw/Q19td648xM2tKZQ+oSpoOfB74z9V8sKSNkgYlDQ4PD5f9/u17jzIyOkYAI2fH2L73aDXlmJllSiWzZS4DlgNPSfoHYAnwM0kLgIPA0pLXLkmOvUVEbIqI/ojo7+mZ8Nuzb2tgxVw68sXy23JiYMXcsn+HmVlWlR3uEfFMRFwaEcsiYhnFrpdrI+Iw8AjwyWTWzABwIiJq3iUD0NfbzbdvexcdbTluunK+u2TMzEpMZirk/cBPgVWSDkja8DYv/z6wFxgCvgH8Tk2qvIC+ZXP4tSWzOPr6SD0/xsys6Vx0QDUiPnqR88tKtgO4o/qyJm/1gpn8xVMvExGcN7hrZtaymv4bqqsXzOTkqVEOnTjV6FLMzFKj+cN9YRcAzx9+rcGVmJmlR9OH+xXzZwKw6/DJBldiZpYeTR/us6a1s3j2NLfczcxKNH24A6xaMJPdhxzuZmbnZCLcVy+YyYvDrzMyOtboUszMUiET4b5qwUxGx4IXh19vdClmZqmQiXC/Mpkxs9uDqmZmQEbCffm8GbS3id0eVDUzAzIS7u1tOS6/1IOqZmbnZCLcAa5cMNPTIc3MEpkJ9xmFPIdPnuInzx9pdClmZg2XiXDfse8YDzxRfADU7d/a4cfumVnLy0S4b997lNGx4hz3M34qk5lZNsK99KlMOfmpTGZmmQj3vt5uttw2wLJ5M1g6Z5qfymRmLS8T4Q7FgF971QIOHv+VlyEws5aXmXAHWLOoizNng6EjXobAzFpbtsI9WYZg5yEvQ2BmrS1T4b583gw623PsfNnhbmat7aLhLuk+SUckPVty7A8l7Zb0tKT/JWl2ybnPSRqS9Lyk36pX4RNpy4nVC7rYeejEVH6smVnqTKbl/k1g7XnHHgWujohfB14APgcgaQ1wK3BV8p4/ldRWs2on4apFXex8+SQRMZUfa2aWKhcN94h4DHj1vGM/jIjRZHc7sCTZXgd8JyJOR8QvgCHgnTWs96LWLOri5KlRDh7/1VR+rJlZqtSiz/3fAn+VbC8GXio5dyA59haSNkoalDQ4PDxcgzKKxgdV3e9uZi2sqnCX9AVgFNhS7nsjYlNE9EdEf09PTzVl/H9WL+giJ8+YMbPWlq/0jZJ+G/ggcGP8Ywf3QWBpycuWJMemzLSONpbPm+GWu5m1tIpa7pLWAp8FPhQRb5acegS4VVJB0nJgJfD31ZdZnjWLZrnlbmYtbTJTIe8HfgqsknRA0gbgT4CZwKOSfi7pvwNExHPAg8BO4AfAHRFxtm7VX8CahV0cOPYrvvzD5738r5m1pIt2y0TERyc4fO/bvP6LwBerKapaHXkB8Cfbhtj0N3vZctuAFxMzs5aSqW+onvPqGyMAjAWcGfX67mbWejIZ7jesno+S7fZ8zuu7m1nLyWS49/V28y9WzmNGRxtbNrzLXTJm1nIyGe4AN145nzdGzrJw9rRGl2JmNuUyG+7X/JPiWmZP7j/e4ErMzKZeZsN99YIuCvkcT+73VEgzaz2ZDfeOfI5fWzyLJ19yy93MWk9mwx2KXTPPHDzhZ6qaWcvJeLh3MzI6xi4vRWBmLSbj4X5uUNX97mbWWjId7gtnTWNBV6f73c2s5WQ63KHYevd0SDNrNS0R7vtffZM//OvdXiHSzFpG5sN9ekdx4cs/3fYiH79nuwPezFpC5sP96OunAQi8QqSZtY7Mh/t7V/aQS5aI9AqRZtYqMh/ufb3dfPLdvQB89dZrvEKkmbWEzIc7wMfeVQz3Y8lDPMzMsq4lwn3lpZewoKuTx/YMN7oUM7MpMZkHZN8n6YikZ0uOzZH0qKQ9yc/u5LgkfU3SkKSnJV1bz+InSxLXXTGPv93zCqNnvc6MmWXfZFru3wTWnnfsLmBrRKwEtib7ADcDK5M/G4G7a1Nm9a67ooeTp0Z56sCJRpdiZlZ3Fw33iHgMePW8w+uAzcn2ZuCWkuPfiqLtwGxJC2tVbDXee/k8JHjsBXfNmFn2VdrnPj8iDiXbh4H5yfZi4KWS1x1Ijr2FpI2SBiUNDg/XP3BnT+/g15fM5m/c725mLaDqAdWICIrfESr3fZsioj8i+nt6eqotY1J+Y+U8ntx/nC//8Hl/U9XMMq3ScP/lue6W5OeR5PhBYGnJ65Ykx1JhflcnAfy3Hw15KQIzy7RKw/0RYH2yvR54uOT4J5NZMwPAiZLum4Z79c3iPHcvRWBmWTeZqZD3Az8FVkk6IGkD8AfA+yTtAW5K9gG+D+wFhoBvAL9Tl6or9M8vKw6qgpciMLNsy1/sBRHx0QucunGC1wZwR7VF1Utfbzc3XTmfx14YZsuGd3kpAjPLrJb4hmqp61f1cHp0jEu7OhtdiplZ3bRcuF+1aBYAz73sh2abWXa1XLivmj+TnGDny/6mqpllV8uF+7SONi7ruYSdh9xyN7PsarlwB7hqUZe7Zcws01oy3Ncs6uLQiVO86vXdzSyjWjLc/3FQ1f3uZpZNLRruXQDsdNeMmWVUS4b77OkdLJ49zf3uZpZZLRnuAFcu7HK3jJllVsuG+1WLutj7yhu8OTLa6FLMzGqupcM9Ar74l7u89K+ZZU7LhvtYFJ8v8u3H93ttdzPLnJYN9xeHXwe8truZZVPLhvvAinm0txUXd89JXtvdzDKlZcO9r7eb79w+wNLuaXTkcyyfN6PRJZmZ1UzLhjtA37I53Pvb/4zTo2Pc+dBTfH3bkPvezSwTLvokpqy7Yv5Mbr56AX/x9CG27j5CRz7HltsG/JQmM2tqLd1yP2dFT7FLZiw8uGpm2VBVuEv6fUnPSXpW0v2SOiUtl/S4pCFJD0jqqFWx9XLdFZfSlisOrvrB2WaWBRWHu6TFwO8C/RFxNdAG3Ap8CfhKRFwOHAM21KLQeurr7eYz77sCgM+//0p3yZhZ06u2WyYPTJOUB6YDh4AbgIeS85uBW6r8jCnxiXf30t4mDh7/VaNLMTOrWsXhHhEHgT8C9lMM9RPADuB4RJxbsOUAsLjaIqfCzM52+nvn8JPnhxtdiplZ1arplukG1gHLgUXADGBtGe/fKGlQ0uDwcDoC9fpVPew+/BqHTrj1bmbNrZpumZuAX0TEcEScAb4LvAeYnXTTACwBDk705ojYFBH9EdHf09NTRRm185urLwXgx269m1mTqybc9wMDkqZLEnAjsBPYBnw4ec164OHqSpw6Ky+9hEWzOvnx80caXYqZWVWq6XN/nOLA6c+AZ5LftQm4E/i0pCFgLnBvDeqcEpK4fvWl/N3QUUZGxxpdjplZxaqaLRMR/yUiVkfE1RHxiYg4HRF7I+KdEXF5RPzriDhdq2KnwvVX9PD66VH+08PPeCkCM2ta/obqeaYXisMFDz5xwOu8m1nTcrif56mXjgPFdd5HvBSBmTUph/t5BlbMpTNf/McyFrCgq7PBFZmZlc/hfp6+3m623D7Av/uNy5gzvYMv/uUuvvRXu909Y2ZNxeE+gb7ebu68eTVf+MCVvPrmCHf/5EX3v5tZU3G4v43DJ0+hZNv972bWTBzub2NgxVw6kv53P2fVzJqJw/1t9PV28+3bB7isZwYzO/O8Y8msRpdkZjYpDveL6Ovt5s61qzn25hmvOWNmTcPhPgm/ufpS5l1S4IHBlxpdipnZpDjcJ6G9LceH+5bwo91HOPLaqUaXY2Z2UQ73SfpI/xLOjgX/8cGnPSXSzFLP4T5Jx948gwSP7Rn2nHczSz2H+yRt33uUiOL2Gc95N7OUc7hP0sCKubTlil9pas/nPOfdzFLN4T5Jfb3d/NZV8ynkc2y5bYC+3u5Gl2RmdkEO9zIs7Z6OhIPdzFLP4V6GQnsbp0fHiHOd72ZmKeVwL0MhnyMCRs76+apmlm4O9zIUkkXETvvh2WaWclWFu6TZkh6StFvSLknvljRH0qOS9iQ/M9NBXWhvA+D0GYe7maVbtS33rwI/iIjVwDuAXcBdwNaIWAlsTfYzoXO85X62wZWYmb29isNd0izgOuBegIgYiYjjwDpgc/KyzcAt1RaZFuda7qfccjezlKum5b4cGAb+TNKTku6RNAOYHxGHktccBuZP9GZJGyUNShocHm6OpXQLbrmbWZOoJtzzwLXA3RFxDfAG53XBRHHO4ITzBiNiU0T0R0R/T09PFWVMHQ+omlmzqCbcDwAHIuLxZP8himH/S0kLAZKfR6orMT06PaBqZk2i4nCPiMPAS5JWJYduBHYCjwDrk2PrgYerqjBFzrXcT7lbxsxSLl/l+/8DsEVSB7AX+DcU/4fxoKQNwD7gI1V+RmoU8m65m1lzqCrcI+LnQP8Ep26s5vemVaHdA6pm1hz8DdUyjPe5e0DVzFLO4V6G8dkyZ9xyN7N0c7iXwVMhzaxZONzLMD6g6nA3s5RzuJehvU3k5G4ZM0s/h3sZJFHIt3HKLXczSzmHe5kK7Tm33M0s9RzuZSrkc+5zN7PUc7iXqTN5jqqZWZo53MtUyOc45W4ZM0s5h3uZCnm33M0s/RzuZepsz3ltGTNLPYd7mQr5Nq8KaWap53AvUyGf83ruZpZ6DvcyFee5u+VuZunmcC9TpwdUzawJONzLVPCAqpk1AYd7mQr5Nk65W8bMUs7hXqbi8gNuuZtZulUd7pLaJD0p6XvJ/nJJj0sakvRA8vDszCgkyw9ERKNLMTO7oFq03D8F7CrZ/xLwlYi4HDgGbKjBZ6RGIZ8jAkbOumvGzNKrqnCXtAT4AHBPsi/gBuCh5CWbgVuq+Yy08aP2zKwZVNty/2Pgs8C5pJsLHI+I0WT/ALB4ojdK2ihpUNLg8PBwlWVMnUJ78qg9D6qaWYpVHO6SPggciYgdlbw/IjZFRH9E9Pf09FRaxpTrHG+5e1DVzNIrX8V73wN8SNL7gU6gC/gqMFtSPmm9LwEOVl9mepxruXs6pJmlWcUt94j4XEQsiYhlwK3AjyLi48A24MPJy9YDD1ddZYoU3HI3syZQj3nudwKfljREsQ/+3jp8RsN4QNXMmkE13TLjIuLHwI+T7b3AO2vxe9Oo0wOqZtYE/A3VMp1ruXvZXzNLM4d7mQp5t9zNLP0c7mUqtHtA1czSz+FepvE+dw+omlmKOdzLND5b5oxb7maWXg73MnkqpJk1A4d7mcYHVB3uZpZiDvcytbeJnNwtY2bp5nAvk6Tio/bccjezFHO4V6DQnnPL3cxSzeFegeJzVN1yN7P0crhXoDN5jqqZWVo53CtQyOc45W4ZM0sxh3sFCnm33M0s3RzuFSj2ubvlbmbp5XCvQGd7m1eFNLNUc7hXoJDPeT13M0s1h3sFivPc3XI3s/RyuFeg0wOqZpZyFYe7pKWStknaKek5SZ9Kjs+R9KikPcnP7tqVmw6Fdg+omlm6VdNyHwU+ExFrgAHgDklrgLuArRGxEtia7GdKId/GKXfLmFmKVRzuEXEoIn6WbL8G7AIWA+uAzcnLNgO3VFtk2ngqpJmlXU363CUtA64BHgfmR8Sh5NRhYH4tPiNNCsnyAxHR6FLMzCZUdbhLugT4c+D3IuJk6bkopt+ECShpo6RBSYPDw8PVljGlCvkcEXDmrMPdzNKpqnCX1E4x2LdExHeTw7+UtDA5vxA4MtF7I2JTRPRHRH9PT081ZUy5c4/a81x3M0urambLCLgX2BURXy459QiwPtleDzxceXnpVGhPHrXnQVUzS6l8Fe99D/AJ4BlJP0+OfR74A+BBSRuAfcBHqisxfTrHH5LtlruZpVPF4R4RfwvoAqdvrPT3NoPxlru/yGRmKeVvqFZgvM/da7qbWUo53CtQGO+WccvdzNLJ4V6BTg+omlnKOdwrUPCAqpmlnMO9AoV8seXu9WXMLK0c7hUotLvlbmbp5nCvQKenQppZyjncKzDe5+6pkGaWUg73CngqpJmlncO9AucGVB3uZpZWDvcKtLeJnNwtY2bp5XCvgKTio/bccjezlHK4V6jQnnPL3cxSy+FeoeJzVN1yN7N0crhXqDN5jqqZWRo53CtUyOe85K+ZpZbDvUKFvFvuZpZeDvcKFfvc3XI3s3RyuFeos73N67mbWWo53CtUyOc45Za7maVU3cJd0lpJz0saknRXvT6nUYrz3N1yN7N0ytfjl0pqA74OvA84ADwh6ZGI2FmPz2uE10+PcuS10+zYdwyA7XuPMrBi7qS2+3q72bHvWFnvmcrttNfXTLWmvb5mqjXt9VVba19vN7VUl3AH3gkMRcReAEnfAdYBmQj3HfuO8XdDRzk7Fvyru/8vAiI5d7FtAfNmdvDKayOTfs9Ubqe9vmaqNe31NVOtaa+vmlqh2BOw5baBmgZ8vbplFgMvlewfSI6Nk7RR0qCkweHh4TqVUR/b9x4lIsb3o+TcxbYDiCjvPVO5nfb6mqnWtNfXTLWmvb5qag3gzOgY2/cepZYaNqAaEZsioj8i+nt6ehpVRkUGVsylI5+jTdDRprK2O9tzfPp9q+hsr+z99d5Oe33NVGva62umWtNeX7W1tudz4900tVKvbpmDwNKS/SXJsUzo6+1my20DVfXHrVows+H9gs1aXzPVmvb6mqnWtNdXba217nNXafdCzX6plAdeAG6kGOpPAB+LiOcmen1/f38MDg7WvA4zsyyTtCMi+ic6V5eWe0SMSvr3wF8DbcB9Fwp2MzOrvXp1yxAR3we+X6/fb2ZmF+ZvqJqZZZDD3cwsgxzuZmYZ5HA3M8ugukyFLLsIaRjYV+Hb5wGv1LCcZtGK192K1wyted2teM1Q/nX3RsSE3wJNRbhXQ9LgheZ5ZlkrXncrXjO05nW34jVDba/b3TJmZhnkcDczy6AshPumRhfQIK143a14zdCa192K1ww1vO6m73M3M7O3ykLL3czMzuNwNzPLoKYO96w/hBtA0lJJ2yTtlPScpE8lx+dIelTSnuRnbReDTglJbZKelPS9ZH+5pMeTe/6ApI5G11hLkmZLekjSbkm7JL27Fe61pN9P/v1+VtL9kjqzeK8l3SfpiKRnS45NeH9V9LXk+p+WdG05n9W04V7yEO6bgTXARyWtaWxVdTEKfCYi1gADwB3Jdd4FbI2IlcDWZD+LPgXsKtn/EvCViLgcOAZsaEhV9fNV4AcRsRp4B8Vrz/S9lrQY+F2gPyKuprhM+K1k815/E1h73rEL3d+bgZXJn43A3eV8UNOGOyUP4Y6IEeDcQ7gzJSIORcTPku3XKP7HvpjitW5OXrYZuKUxFdaPpCXAB4B7kn0BNwAPJS/J1HVLmgVcB9wLEBEjEXGcFrjXFJcfn5Y86Gc6cIgM3uuIeAx49bzDF7q/64BvRdF2YLakhZP9rGYO94s+hDtrJC0DrgEeB+ZHxKHk1GFgfoPKqqc/Bj4LjCX7c4HjETGa7Gftni8HhoE/S7qi7pE0g4zf64g4CPwRsJ9iqJ8AdpDte13qQve3qoxr5nBvKZIuAf4c+L2IOFl6LorzWTM1p1XSB4EjEbGj0bVMoTxwLXB3RFwDvMF5XTAZvdfdFFupy4FFwAze2nXREmp5f5s53DP9EO5SktopBvuWiPhucviX5/6Klvw80qj66uQ9wIck/QPFLrcbKPZHz07+6g7Zu+cHgAMR8Xiy/xDFsM/6vb4J+EVEDEfEGeC7FO9/lu91qQvd36oyrpnD/QlgZTKi3kFxAOaRBtdUc0k/873Aroj4csmpR4D1yfZ64OGprq2eIuJzEbEkIpZRvLc/ioiPA9uADycvy9R1R8Rh4CVJq5JDNwI7yfi9ptgdMyBpevLv+7nrzuy9Ps+F7u8jwCeTWTMDwImS7puLi4im/QO8H3gBeBH4QqPrqdM1vpfiX9OeBn6e/Hk/xf7nrcAe4P8Acxpdax3/GVwPfC/ZXgH8PTAE/E+g0Oj6anyt/xQYTO73/wa6W+FeA/8V2A08C/wPoJDFew3cT3Fc4QzFv6ltuND9BURxRuCLwDMUZxNN+rO8/ICZWQY1c7eMmZldgMPdzCyDHO5mZhnkcDczyyCHu5lZBjnczcwyyOFuZpZB/w+ZYGzfuamgewAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ], "source": [ "plt.plot(np.linalg.svd(x0)[1], '.-')" ] }, { "cell_type": "markdown", "id": "baking-filling", "metadata": { "id": "baking-filling" }, "source": [ "Random selection of $p$ indices." ] }, { "cell_type": "code", "source": [ "def indices(n,p): return np.random.permutation(n**2)[0:p]" ], "metadata": { "id": "ib3a5jIycHdq" }, "id": "ib3a5jIycHdq", "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "id": "blind-sterling", "metadata": { "id": "blind-sterling" }, "source": [ "Associated forward (measurement) operator $\\Phi$ and its adjoint $\\Phi^\\top$." ] }, { "cell_type": "code", "source": [ "def phi(x,I): return x.flatten()[I]\n", "def phi_adj(y,I):\n", " x = np.zeros(n**2)\n", " x[I] = y\n", " return x.reshape((n,n))\n", "def dotp(x,y): return np.sum( x.flatten()*y.flatten() )" ], "metadata": { "id": "495hi7G-gnFC" }, "id": "495hi7G-gnFC", "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "Check adjointness of $\\Phi$ and $\\Phi^\\top$." ], "metadata": { "id": "xFYkXpdAhedA" }, "id": "xFYkXpdAhedA" }, { "cell_type": "code", "source": [ "p = 12\n", "x = np.random.randn(n,n)\n", "y = np.random.randn(p)\n", "I = indices(n,p)\n", "# must be 0\n", "print( dotp(y,phi(x,I)) - dotp(phi_adj(y,I),x) )" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "DRjpQVILet4Q", "outputId": "b659d88e-44fe-4ac4-bbed-8d3951ec8a5d" }, "id": "DRjpQVILet4Q", "execution_count": null, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "-8.881784197001252e-16\n" ] } ] }, { "cell_type": "markdown", "id": "intermediate-surveillance", "metadata": { "id": "intermediate-surveillance" }, "source": [ "Define the proximal operator of $F(x) = \\iota_{\\Phi \\cdot = y}(x)$ (ie the ortho-projector).\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "systematic-today", "metadata": { "id": "systematic-today" }, "outputs": [], "source": [ "def prox_F(x,y,I): return x+phi_adj(y-phi(x,I),I)" ] }, { "cell_type": "markdown", "source": [ "Check that prox$_F \\circ$ prox$_F$ = prox$_F$ (since it is a projector)." ], "metadata": { "id": "eYQHpFP2LnIm" }, "id": "eYQHpFP2LnIm" }, { "cell_type": "code", "source": [ "x1 = prox_F(x,y,I)\n", "print( np.linalg.norm(x1-prox_F(x1,y,I)) )" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "fc0C9tOxh7AN", "outputId": "0fd78b7c-c636-482e-b23f-68ab8d8d2d56" }, "id": "fc0C9tOxh7AN", "execution_count": null, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "2.1677797545656835e-16\n" ] } ] }, { "cell_type": "markdown", "source": [ "Compute the proximal operator of $G(x):=\\|x\\|_*$ the nuclear norm, which is the soft-thesholding of the singular values." ], "metadata": { "id": "mwOxb8ygLzRx" }, "id": "mwOxb8ygLzRx" }, { "cell_type": "code", "execution_count": null, "id": "upper-harvest", "metadata": { "id": "upper-harvest" }, "outputs": [], "source": [ "def soft_thresholding(x, u): return np.minimum(x+u, np.maximum(0,x-u))" ] }, { "cell_type": "markdown", "source": [ "Displays the soft thresholding." ], "metadata": { "id": "WWpANhUzL5Th" }, "id": "WWpANhUzL5Th" }, { "cell_type": "code", "source": [ "t = np.linspace(-3,3,1000)\n", "plt.plot(t, soft_thresholding(t,1.1));" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 184 }, "id": "qFBCditCiO1B", "outputId": "4a1c36c5-2c78-40c4-9b5d-a788fe106174" }, "id": "qFBCditCiO1B", "execution_count": 1, "outputs": [ { "output_type": "error", "ename": "NameError", "evalue": "ignored", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msoft_thresholding\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1.1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m;\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'np' is not defined" ] } ] }, { "cell_type": "code", "execution_count": null, "id": "flying-monitoring", "metadata": { "id": "flying-monitoring" }, "outputs": [], "source": [ "def prox_G(x, gamma): \n", " u,s,vh = np.linalg.svd(x)\n", " soft_s=soft_thresholding(s,gamma)\n", " return (u*soft_s)@vh" ] }, { "cell_type": "code", "source": [ "_,a,_ = np.linalg.svd(x)\n", "_,b,_ = np.linalg.svd(prox_G(x,6))\n", "plt.plot(a)\n", "plt.plot(b);" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 265 }, "id": "5jnjaU8Tiirx", "outputId": "3a62a387-cc2e-4bb4-94d2-475ce80a9586" }, "id": "5jnjaU8Tiirx", "execution_count": null, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3wVVf7/8dcnFUICCSUkAUJAQpNOQLogiBRXsIOIoGjsYlldXfe7rq6url0sKCJgASwIiogiRYr0AJEaCB0CKRBIg4SU8/tjrr/NsgmEtLnl83w87iP3zp175z2P0Q9zz5w5R4wxKKWUcl9edgdQSilVtbTQK6WUm9NCr5RSbk4LvVJKuTkt9Eop5eZ87A5Qkvr165uoqCi7YyillMvYtGnTCWNMg5Lec8pCHxUVRVxcnN0xlFLKZYjIodLeu2jTjYg0EZFfRWSniOwQkYmO5XVFZLGIJDr+hpTy+XGOdRJFZFz5d0MppVR5lKWNvgB4whjTFugBPCgibYGngaXGmGhgqeP1fxGRusBzwBVAd+C50v5BUEopVTUuWuiNMceNMZsdz7OAXUAjYATwqWO1T4GRJXz8GmCxMSbdGHMKWAwMqYzgSimlyuaSet2ISBTQGVgPNDTGHHe8lQw0LOEjjYAjxV4fdSwr6btjRSROROLS0tIuJZZSSqkLKHOhF5FA4FvgUWNMZvH3jDVgToUGzTHGTDHGxBhjYho0KPHCsVJKqXIoU6EXEV+sIj/TGDPXsThFRMId74cDqSV8NAloUux1Y8cypZRS1aQsvW4E+ATYZYx5s9hb84E/etGMA74v4eOLgMEiEuK4CDvYsUwppVQ1KcsZfW9gLHCViMQ7HsOAV4CrRSQRGOR4jYjEiMhUAGNMOvBPYKPj8YJjWaXLKyhkysp9bDxYJV+vlFIuS5xxPPqYmBhzqTdMnT1XyIDXl9Owtj/zHuiNl5dUUTqllHI+IrLJGBNT0ntuM9ZNTT9vnhjckt+PZrBg2/GLf0AppTyE2xR6gBu6NKZ1WBCv/pxAXkGh3XGUUsopuFWh9/YS/jqsDUdPneXztaUO+6CUUh7FrQo9QL+WDegbXZ93l+0l40y+3XGUUsp2blfoAf46rA2Zufk8/nU8p3LO2R1HKaVs5ZaFvk14bf42vC0r9qQx+O2VLEtIsTuSUkrZxi0LPcCEPs34/qHe1Kvlx10z4nhxwU6csSupUkpVNbct9ACXR9Th+4d6c0fPpkz97QD/WrhLi71SyuM45QxTlcnfx5vnr7scAT5edYAavt48MbiV3bGUUqrauH2hBxARnvvT5ZwrLOLdZXup5e/DfVdeZncspZSqFm7ddFOcl5fw0sj2DG8fzuuLdrM7OcvuSEopVS08ptCDVez/ObIdQTV8eHbeNoqKtL1eKeX+PKrQA9St5cczw9oQd+gU32w6cvEPKKWUi/O4Qg9wc9fGdG9Wl5d/SuBkdp7dcZRSqkp5ZKEXEf51fTty8gr4+/wd2oSjlHJrHlnoAVqEBvHooJb8uPU4T8/dSqEWe6WUm7po90oRmQZcC6QaY9o5ln0F/NEZPRg4bYzpVMJnDwJZQCFQUNqg+HZ5oP9l5BUUMWlpIgWFhtdu7oi3TliilHIzZelHPwN4D/jsjwXGmFv/eC4ibwAZF/j8AGPMifIGrEoiwuNXt8TXS3hj8R68vITXb+5odyyllKpUF226McasBEqciNUxcfgtwOxKzlWtHh4Yzb39mjNn01ESU7R/vVLKvVS0jb4vkGKMSSzlfQP8IiKbRCS2gtuqUrH9muPn7cUX63TCEqWUe6looR/Nhc/m+xhjugBDgQdFpF9pK4pIrIjEiUhcWlpaBWNdunqB/gzvEM63m5PIySuo9u0rpVRVKXehFxEf4Abgq9LWMcYkOf6mAvOA7hdYd4oxJsYYE9OgQYPyxqqQ23s0JTuvgO/jj9myfaWUqgoVOaMfBCQYY46W9KaI1BKRoD+eA4OB7RXYXpXrEhlM2/DafLb2oA5nrJRyGxct9CIyG1gLtBKRoyIywfHWKM5rthGRCBFZ6HjZEPhNRH4HNgA/GmN+rrzolU9EGNuzKQnJWWw+fMruOEopVSku2r3SGDO6lOXjS1h2DBjmeL4fcLm+iiM6RfCvH3fx+dpDdG1a1+44SilVYR57Z2xpAvx8uCmmMT9sPc5P247bHUcppSpMC30Jnhjcik5Ngnl49hYW7Ui2O45SSlWIFvoSBPr7MOPObrRrVIeHZm1m4bbj5OYX2h1LKaXKxSOmEiyPoBq+fDahO2OnrueBmZsBCA7wJTo0kHdHdyGsTg2bEyqlVNloob+A2jV8mXlPD37enkxyxlmSM3OZuzmJB2dtZvY9PfDz0R9ESinnp4X+IgL9fbipa+P///qKZvV4ePYWXv5pF8/96XIbkymlVNnoKekl+lPHCO7sHcX01QeZ/7veQauUcn5a6Mvhr8PaENM0hKe/3crqvU45ArNSSv1/WujLwdfbi/fHdKFRcE3umLaBT347oEMmKKWclhb6cmpYuwbzHuzNwNah/HPBTp745neyddRLpZQT0kJfAYH+Pnx4e1ceHRTN3M1J9H9tOV9uOKzzzyqlnIoW+gry8hIeHdSSeQ/0omm9AJ6eu43hk1axbv9Ju6MppRSghb7SdI4MYc59PXn/ti5k5RYwaso6Hv86nhPZeXZHU0p5OC30lUhEGN4hnCWPX8mDAy7jh9+PcdXry4k7WOKUu0opVS200FeBmn7ePHlNa36a2JfaNX3523fbtd1eKWUbLfRVqEVoEM8MbUNCchbfbipxIi6llKpyWuir2LD2YXSJDOb1X3Zz5px2v1RKVb+yTCU4TURSRWR7sWX/EJEkEYl3PIaV8tkhIrJbRPaKyNOVGdxViAjPDm9LalYeU1butzuOUsoDleWMfgYwpITlbxljOjkeC89/U0S8gfeBoUBbYLSItK1IWFfVtWkIw9uH89GK/exNzdK7aJVS1aosc8auFJGocnx3d2CvY+5YRORLYASwsxzf5fKeGtKKxbtSGPTmSoIDfGkTVptR3ZswolMju6MppdxcRdroHxKRrY6mnZAS3m8EHCn2+qhjWYlEJFZE4kQkLi0trQKxnFPTerVY+Ehf/jnicoa2C+dEdh4Tv4zn8a/jydGhE5RSVai849FPBv4JGMffN4C7KhLEGDMFmAIQExPjlm0bLUIDaREaCEBBYRHvLtvLu8sS2XL4NO/f1oW2EbVtTqiUckflOqM3xqQYYwqNMUXAx1jNNOdLApoUe93YsUwBPt5ePHZ1S2bd04Mz5wq4beo69qRk2R1LKeWGylXoRSS82Mvrge0lrLYRiBaRZiLiB4wC5pdne+6sR/N6fHNvL/y8vbh96noOnzxjdySllJspS/fK2cBaoJWIHBWRCcCrIrJNRLYCA4DHHOtGiMhCAGNMAfAQsAjYBXxtjNlRRfvh0iLrBfD5hCvIKyji9k/Wk5qZa3ckpZQbEWfs6hcTE2Pi4uLsjlHt4o+cZszH66hT05d/39SBvtEN7I6klHIRIrLJGBNT0nt6Z6wT6dQkmNmxPajh583YTzbw7LxtOpmJUqrCtNA7mQ6Ng1n4SF/u7tOMWRsOM/jNFSzYekxvslJKlZsWeidUw9ebv13blm/u7UmdAD8emrWF0R+vY9fxTLujKaVckBZ6JxYTVZcFD/fhxZHtSEjOYug7q4j9LI7fj5y2O5pSyoXoxVgXkXEmn2mrDzBjzUEyzuYzuG1D3r2tM/4+3nZHU0o5Ab0Y6wbqBPjy2NUtWf30VTxxdUt+2ZnCM3O3adu9UuqiyjsEgrJJoL8PDw+MxgBvLt5Di9BAHujfwu5YSiknpoXeRT18VQv2pmbz6s+7aV4/kCHtwuyOpJRyUu7TdJObCXPugm1z7E5SLUSEV2/qQKcmwdw/cxM3f7iGqav2k3T6rN3RlFJOxn0KvV8tSNsDS56Hgjy701SLGr7ezLizGxMHRpOdV8iLP+6i/2u/sirR/YZ5VkqVn/sUei9vGPwCZByGDVPsTlNtggP8eHRQS36a2Jflf+7PZQ0Cuf+Lzew8pn3ulVIW9yn0AJddBZcNhJWvwZl0u9NUu6j6tZh+ZzcC/X24c8YGjmkzjlIKdyv0AFe/YLXXr3rD7iS2CK9Tkxl3deNMXiF3Tt/IiWzPaMZSSpXO/Qp9WDvoPMZqvjl10O40tmgdVpuPxnblUHoOI99fTaJOaKKUR3O/Qg8w4FkQb5g1Co7/bncaW/RqUZ+vYnuSm1/EDZPX8FviCbsjKaVs4p6FvnYEjJoJZ9Ph46tgxWtQ6HnD/XZsEsx3D/Yiok5Nxk3fwNPfbuXQyRy7YymlqllZZpiaJiKpIrK92LLXRCRBRLaKyDwRCS7lswcdM1HFi0j1Dl7TYiA8sA7ajoRfX4SZN0K+512cbBwSwJz7ezLmikjmbkliwOvLeeyreG27V8qDlOWMfgYw5Lxli4F2xpgOwB7gmQt8foAxplNpg+1UqYC6cNMncN17sH8FfDXWY/rYFxdUw5cXRrTjt6cGcHff5izcdpwbJ6/h4Ak9u1fKE1y00BtjVgLp5y37xTEnLMA6oHEVZKs8XcbCn96BvYvhmzuhMN/uRLYIrV2Dvw5rw5exPcg8m8+Nk9fokMdKeYDKaKO/C/iplPcM8IuIbBKR2At9iYjEikiciMSlpVXBnZ1dx8Gw12H3j/DDo5X//S6kc2QIc+7vRU0/b0ZNWcdTc37n87UHiT9ymqIiHQ1TKXdTpvHoRSQKWGCMaXfe8meBGOAGU8IXiUgjY0ySiIRiNfc87PiFcEFVOh790hesPva3z7Xa8T1YamYu//f9djYcSOfUGetXTp8W9Xnvts4EB/jZnE4pdSkuNB59uQu9iIwH7gUGGmPOlOE7/gFkG2Nev9i6VVroC/Jgci8oKoQH1oJvzarZjgsxxnAsI5dfdiTz8sIEwurUYModXWkdVtvuaEqpMqr0iUdEZAjwFHBdaUVeRGqJSNAfz4HBwPaS1q1WPv4w/E04dQBWvWl3GqcgIjQKrsmdvZvx5b09yM0v5IYP1vCXOVt585fdfL7ukA6noJQLu+gZvYjMBvoD9YEU4DmsXjb+wEnHauuMMfeJSAQw1RgzTESaA/Mc7/sAs4wxL5UlVLVMJTg3FrbPhfvXQIOWVbstF5OSmcszc7exLSmDk9l5FBkIDfLny9geNG8QaHc8pVQJKtx0U92qpdBnp8F7XaFeNNzxPfhrAStJQWERu45nMX76Bny9vfjq3h40rVfL7lhKqfPonLElCWxg9a8/tgW+uNEaCE39Dx9vL9o3rsPMe64gr6CQ0VPWcfjkRS/JKKWciOcWeoC218FN0yApDj6/Hs5qn/LStA6rzcy7e3Amv5Ah76zkrcV7yM7zvGEllHJFnl3oAS4fCbd8Zg1+NutWcMKmLGfRNqI28x/sw4DWobyzNJH+r/3K1FX7yTjrmTegKeUqPLeN/nwbp8KPT8D4hRDVu3q37YK2HD7FKz8lsP5AOjV9vRnRKYI7ezejVViQ3dGU8kjaRl8WHW+DGnUg7hO7k7iEzpEhfHVvTxY83IcRnSL4Lj6J4ZNW8d6yRAoKi+yOp5QqRgv9H/wCoNPtsHM+ZKXYncZltGtUh1du7MDapwcytH04r/+yh1v1gq1STkULfXExd0FRPmz5zO4kLieklh/vju7MO6M6sScli+ve/409OrOVUk5BC31x9VtA8/4QN8MjJyqpDCM6NWLBw33w8/Zi7CfrOZKuZ/ZK2U0L/fm63Q2ZRyFxkd1JXFbTerX4fMIVnD1XyNhP1pOW5XlzACjlTLTQn6/lUAiKsCYXd8IeSa6iVVgQ0+/sTkpmHuOnbyBH+9wrZRst9Ofz9oEr7oX9y2HFv+1O49K6Ng3hg9u7sOt4Jo9/Ha9j3StlEy30Jen1CHQaA8tftiYWV+U2oFUozw5vy6IdKby9ZI/dcZTySD52B3BKXl5w3bvWmPW/vgj5OY4mnTAICgcfnZTjUtzVO4qE45lMWraX6IZB/KljhN2RlPIoemfshRQVwrz7YNvX/1lWo4412mVEZ/tyuaC8gkJu+3g9mw6dIjjAl6h6tWgdFsTDA6NpFKyTvyhVUTpMcUUYA6k7IfOY9VjxqnXGH7sCAuranc6lZJzJ55tNR9h/IodDJ3PYfOg03l7C34a34dZuTRARuyMq5bK00Femo3EwbYjV3/62r62ir8rlSPoZnpqzlbX7T9I3uj7je0XRu0V9avh62x1NKZdT4bFuRGSaiKSKyPZiy+qKyGIRSXT8DSnls+Mc6ySKyLjy7YITaRwDQ1+BvYthpV6orYgmdQOYefcVvDDicn4/cpoJn8YR8+ISHp69hTX7TuCMJyFKuaKyTg7eD8gGPvtjgnAReRVIN8a8IiJPAyHGmL+c97m6QBwQAxhgE9DVGHPqQttz6jN6sJpz5t0LW7+Ga16CHg+ANjtUyLmCItbsO8GiHcn8vD2ZU2fyaR0WxF29m3F9l0b4eusvJ6UupFKabkQkClhQrNDvBvobY46LSDiw3BjT6rzPjHasc6/j9UeO9WZfaFtOX+gBzp2BebGw6wfoOBqufQt89aJiZcjNL2R+/DGmrT5AQnIWHZsE8+6ozkTWC7A7mlJOq6qGKW5ojDnueJ4MNCxhnUbAkWKvjzqWuT6/ALj5MxjwLPw+G6YPhdQEu1O5hRq+3tzSrQk/TezLu6M7sz8tm2GTVvF9fJLd0ZRySZXye9hYPwsq1KAqIrEiEicicWlpaZURq+p5ecGVT8GoWXByP0zuBQsesyYeVxUmIvypYwQ/TexLq7AgJn4Zz/jpG9h6VKd8VOpSVKTQpziabHD8TS1hnSSgSbHXjR3L/ocxZooxJsYYE9OgQYMKxLJB6+HwyBboNgE2fQqTOsPa93UEzErSOCSAr2J78MzQ1sQfOc11763m7k83svnwKb1gq1QZVKSN/jXgZLGLsXWNMU+d95m6WBdguzgWbca6GJt+oW25RBt9adL2wKJnYO8SaNjeartv0s3uVG4jKzefT9ccZMrK/WTmFtChcR3G94ri2g4R+PnoBVvluSp8MVZEZgP9gfpACvAc8B3wNRAJHAJuMcaki0gMcJ8x5m7HZ+8C/ur4qpeMMdMvtj2XLvRg9crZNR9+ehqyjsENU6HDzXancis5eQXM3XyUGWsOsi8th1YNg3jr1k60jahtdzSlbKE3TNklLws+GwEZSfDIZvCrZXcit2OMYfHOFJ79bjsZZ/L58zUtubtPc7y8tLur8iw6Obhd/IPgmpchOxnWvGd3GrckIgy+PIxFj/bjqtah/GthAiM/WM2vCanafq+Ugxb6qhZ5BbQdAavfgaxku9O4rbq1/Jh8exfevrUT6TnnuHPGRq7/YA2r956wO5pSttNCXx0G/QMKz8Gv/7I7iVsTEUZ2bsSyJ/rz8g3tScvKY8zU9dw1YyN7U3WicuW5tNBXh7rNoXssbPncGjbh1CGdprAK+fl4Mbp7JEufuJKnh7Zm44F0rnl7FS8v3EVBYZHd8ZSqdnoxtrqcSYfJva1eOAD+daDLWBj4d/DxtzebmzuRncfri3bz5cYjDGjVgHdv60Kgv865o9yL9rpxFufOQMoOSNkGB1fD9jkQ3hFumg71LrM7ndv7Yt0hnpu/g5YNg5g2PobwOjo2kXIfWuidVcKP8P2DUJgPIydD2+vsTuT2lu9O5aFZWygsMgxo3YBh7cPp3ypUz/CVy9NC78wykuCb8XBsszVmTstr7E7k9vamZjNjzQF+3p7Ciew8wOq1E16nBm3Ca/PM0NbUC9TmNOVatNA7u9xM+Ow6SN0Ft38LUX3sTuQRCosMGw6ks+lQOscycjl++ixr9p2kbi0/3h/ThS6RJc6lo5RT0kLvCnJOWkMdZx6Dcd9Do652J/JI25MyeGDmZo5nnOWJwa0Y0CqUpvUCdHpD5fS00LuKzGMw7RrISoGrnoWeD4GXFpjqlnE2nye+jmfJLmtAVi+BZvVr8eQ1rRjSLtzmdEqVTAu9K8lOtca0T1gAjbvDyA+gfrTdqTyOMYYdxzLZl5bN/rQcftmZwq7jmQxtF8bzIy4nNKiG3RGV+i9a6F2NMbBtDiz8M5zLgZ4PQr8nwT/Q7mQeK7+wiI9X7eftJYnU9PXmoQEtGNuzqTbpKKehhd5VZaXA0uchfiYEhVs3V7W/Bby1K6Bd9qVl8/wPO1m5J42w2jV4ZGA0N8c01snLle200Lu6o3Gw8EmrC2ZwU+jzKHS8DXy1+cAu6/af5NWfE9h8+DRN6tbk4QHRXN+lkRZ8ZRst9O6gqAj2/AyrXoekTdb4OXctgsBQu5N5LGMMy3en8daSPWw9mkFk3QAeGRjNyE4R+GjBV9VMC707MQYSF8PXd0B4B7hjvp7Z28wYw7KEVN5asoftSZk0b1CLRwe15Nr24ToBiqo2VTLxiIi0EpH4Yo9MEXn0vHX6i0hGsXX+Xt7tKQcRaDkYrp8MR9bDDxN1JEybiQgD2zTkh4f68OHtXfHxEh6ZvYVbPlrL3tRsu+MpRbmv6hljdgOdAETEG0gC5pWw6ipjzLXl3Y4qxeXXw4lE+PUlCImCvk+Aj5/dqTyaiDCkXRiD2zbk281HefHHXQx7ZxUTB0VzR8+mBNXwtTui8lCV1X1jILDPGHOokr5PlUW/J+HEHljxCqybDK2GQPubIfpqu5N5NC8v4eaYJvRvFcpz87fz2qLdvLZoN0H+PoQH12B090ju7N3M7pjKg1RKG72ITAM2G2PeO295f+Bb4ChwDPizMWZHKd8RC8QCREZGdj10SP/NKJOiQqvNftcPsPtHOHsKhr0O3e+xO5lyWLP3BNuSMjiekcu2pAw2HTrFk9e04sEBLeyOptxIlV6MFRE/rCJ+uTEm5bz3agNFxphsERkGvGOMuehtnnoxtpwK862LtHt+hlu/gNbD7U6kzlNYZHji63i+iz+mxV5Vqiq5GFvMUKyz+ZTz3zDGZBpjsh3PFwK+IlK/ErapSuLtCzd+AhGdYc4Eq/+9cireXsIbt3RiZKcIXlu0m5cX7iI3v9DuWMrNVUahHw3MLukNEQkTEXE87+7Y3slK2KYqjV8AjP4KghrCrFvgWLzdidR5/ij2o7s34aOV+xny9kpW7kmzO5ZyYxUq9CJSC7gamFts2X0icp/j5U3AdhH5HZgEjDLO2HHf3QQ2gNvngm+ANfRxwo92J1Ln8fYSXr6hA19MuAIR4Y5pGxgzdR0fr9zPruOZ6P8mqjLpDVPuLCsFZo+CY1vgmpegxwNWP3zlVHLzC5m6aj/ztiSxLy0HgJYNA3nr1k5cHlHH5nTKVeidsZ7s3BmYF2v1ygnvaHXJbDUcvPQWfWd07PRZVu6xhlU4lZPPU0NacVfvZnqHrbooLfSerqjIGgHztzchfT80aAPXfwgRnexOpkpxKuccf/l2K7/sTOGKZnV5cEAL+kbXR/QXmSqFFnplKSyAHfNgyT8g9zTc+jlcdpXdqVQpjDF8ufEIby7eQ1pWHtGhgYzvHcV1HSP0Llv1P7TQq/+WeRxm3gRpCTDiA+h4q92J1AXkFRSy4PfjfPLbAXYez6SGrxdD24Uzunsk3ZvVtTuechJa6NX/ys2AL8fAwVUQ1h7CO1n979vdADVD7E6nSmCMIf7Iab7ZdJQf4o+RlVfAA/0v48+DW2kbvtJCr0pRkAerJ8HhNVZ/+7Pp1sQmo2ZBWDu706kLOHuukBcW7GD2hiNc3bYhb93aiUB/nXnMk2mhVxdnDBxeB3PutM72R06Gy0fanUpdgDGGGWsO8s8FO2nZMIgPb+9KVP1adsdSNtFCr8ouKxm+GgtHN1hNOkHhENgQOt0GTXvZnU6VYOWeNB6evYXCIsOrN3VgWPtwuyMpG1T1WDfKnQSFwfgF0O8pq8hnJcOu+fDZSNi3zO50qgT9Wjbgx0f60CI0kAdmbuYf83eQk1dgdyzlRPSMXl3cmXT49E9wci/c9jU0v9LuRKoE5wqKeOWnBKatPkDdWn7E9mvOHT2bEuCnbfeeQJtuVMXlnLCKffoBq/+9Tm7itLYcPsVbSxJZuSeNerX8GNW9CaO6RdKkboDd0VQV0kKvKkd2mlXs03ZByyEw6B8Q2sbuVKoUmw6d4oNf9/Lr7lQM0KdFfW7rHsmgtg3x9dZWW3ejhV5VnvyzsP5DWPUmnMu2pi7sHguNuuqAaU7q2OmzfB13hK82HuF4Ri4Ngvy5NaYJE/o0I6SWzjPsLrTQq8qXcxJWvQGbP7UKfnhH6PUItLtRC76TKigsYvnuNGZtOMzy3ak0qRvAJ+O60SI00O5oqhJooVdVJy8Ltn4FG6ZaTTpt/gTXvgO16tmdTF3A5sOniP0sjryCIj4Y04W+0Q3sjqQqSLtXqqrjHwTd7ob7V8PVL8Dun2FyT9i71O5k6gK6RIbw3YO9iahTk/HTN/LlhsN2R1JVqMKFXkQOisg2EYkXkf85DRfLJBHZKyJbRaRLRbepnJCXN/SeCLG/Qs26MPNm2LPI7lTqAhqHBPDtA73o3aI+T8/dxscr99sdSVWRyjqjH2CM6VTKz4ahQLTjEQtMrqRtKmcU1h7uXgzhHeDrO+DQGrsTqQsI9Pdh6h0xDG8fzksLd/HmL7t1GkM3VB1NNyOAz4xlHRAsInqPtjvzD4Ixc6BOE5g1CpK32Z1IXYCfjxeTRnfmlpjGTFq2l1s/WseUlftITMnSou8mKqPQG+AXEdkkIrElvN8IOFLs9VHHsv8iIrEiEicicWlpaZUQS9mqVn244zur6H9+A2QctTuRugBvL+GVGzrwzNDWZObm86+FCVz91kqGvrOK77YkUVBYZHdEVQEV7nUjIo2MMUkiEgosBh42xqws9v4C4BVjzG+O10uBvxhjSu1Wo71u3EhqAnxyNYREwV0/g5+OrugKjp0+y9KEVD5bc5DE1GwaBdfkqSGtGNHpf87RlJOo0l43xpgkx99UYB7Q/bxVkoAmxV43dixTniC0Ndz4idV8890D1nDIyulFBNdkbI+mLHq0H1PviKF+oB8Tv4xn0tJEbc5xQRUq9CJSS0SC/ngODAa2n7fafOAOR++bHkCGMeZ4RUy/j0IAABBRSURBVLarXEzLwVbXy53fWfPVnj5sTViunJ6XlzCobUPm3N+LG7o04s3Fe/jbd9spLNJi70oqOqxdQ2CeY2Z6H2CWMeZnEbkPwBjzIbAQGAbsBc4Ad1Zwm8oV9XoYUnfC6reth7e/dbZ/9QvQvL/d6dRF+Hp78cbNHQmrXYMPlu8jLSuPSaM7U8PX2+5oqgz0zlhVfYoK4fBaOJEI6fsh4UdI3wcxE6yC76+34ruCGasP8PyCnXSNDGHquBiCA3S8HGegQyAo55R/Fpa9CGvfhzqNofNYawiF0DY6Xo6T+3HrcR77Kp7IegF8eld3GgXXtDuSx9NCr5zb4XVW2/3hdYCBei2g63ir8NcMtjmcKs3afSeJ/SyOAH9vpt7RjfaN69gdyaNpoVeuISvZas7ZNgcOrwHfWtZctVc+BYGhdqdTJUhIzmTCjDhO5uTx5i2ddL5aG2mhV67n+O+w7kPYPgdqBMP1H0KLgXanUiVIy8rj3s/j2Hz4NI9c1YL7+7egpp9epK1uWuiV60rZCXPusoZA7vUIXPU38PG3O5U6T25+IX+du425W5KoV8uPu/s2Z2zPpgT663y11UULvXJt587AL89C3DSo3Rj6/Rk6jQEf7e3hbDYeTOe9ZXtZsSeNIH8fru0Yzo1dGtO1aQiiF9irlBZ65R72/Wr10kmKg+BI6wy/42jtlumEfj9ymk/XHOSn7cmczS+kWf1a3Nk7ipu7NtFmnSqihV65D2MgcTGseAWSNoF/beh8O/R8COroOCzOJievgJ+3J/PF+kNsOXyakABfxvaM4oH+l+nNVpVMC71yP8bA0TjY8BHsmGf10Bn+ujVZuTYROKW4g+l8tHI/i3em0DosiPfHdOGyBvprrLJooVfuLX0/zLsPjqyHy2+A4W9AQF27U6lSLN+dyuNf/05ufiH/ur49IzvrL7HKoHPGKvdWtzmMX2j1yNk1HyZ1ghWvQW6m3clUCfq3CuXHR/pweURtHv0qnnHTNrAnJcvuWG5Nz+iVe0neDr++BLsXQs0Q6PMYdL8XfGvYnUydp6CwiOmrDzJpWSI5eQXc2q0JXSJDqF3Tl5AAPzpHBuPrreeiZaVNN8rzJG2CX/8Fe5dAcFO4+nloO1Lb751Qes45Ji1N5It1hygoNvxxt6gQ3h/ThdAg/Ue6LLTQK8+171f45W+Qst0aDvmm6dp+76Ry8gpIzzlHxtl8tiVl8PwPOwiu6cfk27vQOTLE7nhOTwu98mxFhbBpBvz8tNX/fsw3Vru+cmo7j2US+3kcqZl53Nf/Mu7sFUVILb1JrjR6MVZ5Ni9v6DYB7vgezqTDxwNh/3Kd1tDJtY2ozQ8P9WFQ21AmLU2k97+X8eKCnaRm5dodzeWU+4xeRJoAn2HNMmWAKcaYd85bpz/wPXDAsWiuMeaFi323ntGrKnNyH8y82ZrwJKQZtB0B0YOtM/2gcPDWsVmc0Z6ULCYv38f8349Rw8eLB69qwV29m+lNV8VUSdONiIQD4caYzY55YzcBI40xO4ut0x/4szHm2kv5bi30qkrlZlo3We38DvavAFNoLRcvaNAGRr4PEZ3tzahKdOBEDi/9uIslu1JoUrcmsf0uY3DbhjSsrRdsq6WNXkS+B94zxiwutqw/WuiVMzuTbvXQyTgKmUkQPxtyUmHYa9BlnPbScVK/JZ7gxR93kpBs9b/v2LgOfaLr06FxMB0bBxNWx/MKf5UXehGJAlYC7YwxmcWW9we+BY4Cx7CK/o6LfZ8WemWbnJMw927Yt8y6yzZ6MNSPth41dAYlZ2KMITE1m8U7U1iyK4WtRzModHTPbB0WxCMDoxlyeRheXp7xj3WVFnoRCQRWAC8ZY+ae915toMgYky0iw4B3jDHRpXxPLBALEBkZ2fXQoUMVyqVUuRUVwop/w6o3oSjfWublA70nQr+n9OYrJ5WbX8iOY5nEHznN7A2H2ZuaTeuwIO7p25y+0fUJdfPmnSor9CLiCywAFhlj3izD+geBGGPMiQutp2f0yikUnIPTh+BEIuz8HrZ+CfWiYcR7ENnD7nTqAgqLDAu2HuOdJYnsP5EDQPMGtejfMpTxvaKIrBdgc8LKV1UXYwX4FEg3xjxayjphQIoxxohId2AO0NRcZKNa6JVT2rsEfngMMg5Du5tgwF+h3mV2p1IXUFhk2Hksk7X7T7B230l+23uCwiLD8A4R3NuvOe0auU9zXFUV+j7AKmAbUORY/FcgEsAY86GIPATcDxQAZ4HHjTFrLvbdWuiV08rLhlVvwPoPofAcdB4L/Z7UsfBdRHJGLtNXH2Dm+sNk5xUwqE0oEwe2pH1j1y/4emesUpUtKwVWvQ5x062eOV3GQd/HoXaE3clUGWTm5vPp6oNM/e0AGWfzGdSmIU8PbUWL0CC7o5WbFnqlqsrpw7DydYifCeINMXdCn8chqKHdyVQZ/FHwp6zcz5n8Qsb2aMrEgdEuOdSCFnqlqtqpQ7DyNYifBd5+0P0e6wy/pg7G5QpOZufx1pI9zFp/mEB/H+698jLG9Yoi0N917pTWQq9UdTm5z+qaufVrCAyFYa9D2+vsTqXKaHdyFq/+nMDShFRCAny598rLmNCnmUuMi6+FXqnqdvx3+P4hSN4Kba6zzvBrN7Ie2g/f6cUfOc3bS/awfHcaA1uH8v6YLk4/ro4WeqXsUJgPaybB8n9DYd5/lne8Da57VwdQcwGfrzvE37/fTkzTEKaO60admr52RyrVhQq9/pemVFXx9oW+T0Cn2yFtF2QkwbHNsHEqFOTCDR9rsXdyY3s0JSTAl8e+iufWj9byyfhuNAquaXesS6b/lSlV1YIa/qcXTucx1pDIi/9uvdZi7/Su7RBBnZq+3P/FZoa9s4p/39iBIe3C7I51SZz/CoNS7qb3RLj6BdgxF2bdYvXYUU6tb3QDFjzch8i6Adz3xSb+77vtnMzOu/gHnYS20Stll7hpsOhvYIqg/1+gx4Pg43r9tz3JuYIiXluUwMerDiACXSNDuLptQ3q3qE+b8Np42zhSpl6MVcpZZRyFn/4CCQugVgNrAvPmA6zhkQMb2J1OlSIhOZOftiWzeGcKO49bI7MH+vvQOTKYW2KacG2HcKSa5zLQQq+Us0tcDFu/suayzUkDvyAY/gZ0vNXuZOoijp0+y8aD6Ww8mM5viSc4ePIMXZuG8H/XtqVTk+Bqy6GFXilXUVRk9b3/+Wk4vBba32wVfJ30xCUUFhnmbDrCa4v2cCI7j3v6NuOZoW2qZfKTCxV6vRirlDPx8oKITjBuAQx4FrbPhfe6WyNmnkm3O526CG8v4dZukSx/sj9jrojk41UHeOrbrRQUFl38w1VIC71SzsjbB658Cib8AqGtYekL8GZb+GGiNcetE/4SV/8R6O/DiyPb8diglszZdJQHZ20mr6DQtjzadKOUK0jZCes+gG1zoOAsNGgDXcdbQyt4Ofet+Z5u2m8HeGHBTpo3qMVt3SO5oUtj6lbB6JjaRq+Uu8jNgB3zYMsXcHSj1Tvnxqnahu/kftmRzIcr9rH58Gn8vL0Y2TmCP1/TitCgyhv3SAu9Uu5o4yfw01NQ9zK47Uuo29zuROoididnMWv9IWZtOIy/jzePDopmXK+oShkds8ouxorIEBHZLSJ7ReTpEt73F5GvHO+vF5GoimxPKVVMtwkwdh7kpMLHV1ln+sqptQoL4vkR7Vj0aD9iokJ48cddXDvpN3Yey6zS7Za70IuIN/A+MBRoC4wWkbbnrTYBOGWMaQG8Bfy7vNtTSpWgWT+4ZxmERME342HOXdo7xwU0bxDI9PHdmDK2K+lnzjHi/d/4cMU+CouqpoWlIpOD9wT+YYy5xvH6GQBjzMvF1lnkWGetiPgAyUADc5GNatONUpeosAB+e8ua9MQvAILC7U6kyqigyJCamUt2XgEF/iFEPbmCAL9LH+iuqoYpbgQcKfb6KHBFaesYYwpEJAOoB5woIWQsEAsQGRlZgVhKeSBvH7jySWg1BNa8Z/XMUS7BBwhvCEmnznIs14+aVTDBidOMj2qMmQJMAeuM3uY4SrmmsPZww0d2p1CXSIDGjkdVqMjF2CSgSbHXjR3LSlzH0XRTBzhZgW0qpZS6RBUp9BuBaBFpJiJ+wChg/nnrzAfGOZ7fBCy7WPu8UkqpylXuphtHm/tDwCLAG5hmjNkhIi8AccaY+cAnwOcishdIx/rHQCmlVDWqUBu9MWYhsPC8ZX8v9jwXuLki21BKKVUxOqiZUkq5OS30Sinl5rTQK6WUm9NCr5RSbs4pR68UkTTgUDk/Xp8S7rx1c564z+CZ++2J+wyeud+Xus9NjTElzijvlIW+IkQkrrTxHtyVJ+4zeOZ+e+I+g2fud2XuszbdKKWUm9NCr5RSbs4dC/0UuwPYwBP3GTxzvz1xn8Ez97vS9tnt2uiVUkr9N3c8o1dKKVWMFnqllHJzblPoLzZRubsQkSYi8quI7BSRHSIy0bG8rogsFpFEx98Qu7NWNhHxFpEtIrLA8bqZY9L5vY5J6P3szljZRCRYROaISIKI7BKRnu5+rEXkMcd/29tFZLaI1HDHYy0i00QkVUS2F1tW4rEVyyTH/m8VkS6Xsi23KPRlnKjcXRQATxhj2gI9gAcd+/o0sNQYEw0sdbx2NxOBXcVe/xt4yzH5/CmsyejdzTvAz8aY1kBHrP1322MtIo2AR4AYY0w7rCHQR+Gex3oGMOS8ZaUd26FAtOMRC0y+lA25RaEHugN7jTH7jTHngC+BETZnqhLGmOPGmM2O51lY/+M3wtrfTx2rfQqMtCdh1RCRxsBwYKrjtQBXAXMcq7jjPtcB+mHN64Ax5pwx5jRufqyxhk+v6ZiVLgA4jhsea2PMSqx5Ooor7diOAD4zlnVAsIiUeQZ4dyn0JU1U3simLNVGRKKAzsB6oKEx5rjjrWSgoU2xqsrbwFNAkeN1PeC0MabA8dodj3kzIA2Y7miymioitXDjY22MSQJeBw5jFfgMYBPuf6z/UNqxrVCNc5dC73FEJBD4FnjUGJNZ/D3HdI1u029WRK4FUo0xm+zOUs18gC7AZGNMZyCH85pp3PBYh2CdvTYDIoBa/G/zhkeozGPrLoW+LBOVuw0R8cUq8jONMXMdi1P++Cnn+JtqV74q0Bu4TkQOYjXLXYXVdh3s+HkP7nnMjwJHjTHrHa/nYBV+dz7Wg4ADxpg0Y0w+MBfr+Lv7sf5Dace2QjXOXQp9WSYqdwuOtulPgF3GmDeLvVV8IvZxwPfVna2qGGOeMcY0NsZEYR3bZcaYMcCvWJPOg5vtM4AxJhk4IiKtHIsGAjtx42ON1WTTQ0QCHP+t/7HPbn2siynt2M4H7nD0vukBZBRr4rk4Y4xbPIBhwB5gH/Cs3XmqcD/7YP2c2wrEOx7DsNqslwKJwBKgrt1Zq2j/+wMLHM+bAxuAvcA3gL/d+apgfzsBcY7j/R0Q4u7HGngeSAC2A58D/u54rIHZWNch8rF+vU0o7dgCgtWzcB+wDatXUpm3pUMgKKWUm3OXphullFKl0EKvlFJuTgu9Ukq5OS30Sinl5rTQK6WUm9NCr5RSbk4LvVJKubn/B7b8Fh1uvuZKAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "markdown", "source": [ "Compute the symmetrized proximal operator (rprox)." ], "metadata": { "id": "9d5IuO_sL9Rw" }, "id": "9d5IuO_sL9Rw" }, { "cell_type": "code", "execution_count": null, "id": "designing-first", "metadata": { "id": "designing-first" }, "outputs": [], "source": [ "def rprox_F(x,gamma,y,I):\n", " return 2*prox_F(x,y,I)-x\n", "def rprox_G(x,gamma):\n", " return 2*prox_G(x,gamma)-x" ] }, { "cell_type": "code", "source": [ "def n_norm(x): return np.linalg.norm(x,'nuc')" ], "metadata": { "id": "CGN3cn10jAqE" }, "id": "CGN3cn10jAqE", "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "id": "latest-request", "metadata": { "id": "latest-request" }, "source": [ "Douglas-Rachford's algorithm. " ] }, { "cell_type": "code", "execution_count": null, "id": "appropriate-adventure", "metadata": { "id": "appropriate-adventure" }, "outputs": [], "source": [ "def DR(y, mu, gamma,n_iter,I): \n", " x_t = np.zeros((n,n))\n", " x_list = []\n", " x_t_list = []\n", " nucl_norm=[]\n", " for i in range(n_iter):\n", " x_t_list.append(x_t.copy())\n", " x_t=(1-mu/2)*x_t+mu/2*rprox_G(rprox_F(x_t,gamma,y,I),gamma)\n", " x=prox_F(x_t,y,I)\n", " x_list.append(x.copy())\n", " nucl_norm.append(n_norm(x))\n", " return x, nucl_norm " ] }, { "cell_type": "markdown", "id": "magnetic-property", "metadata": { "id": "magnetic-property" }, "source": [ "Test the DR algorithm." ] }, { "cell_type": "code", "execution_count": null, "id": "interpreted-search", "metadata": { "id": "interpreted-search" }, "outputs": [], "source": [ "mu = 1\n", "gamma = 1\n", "p = int( .5*n*n )\n", "I = indices(n,p)\n", "x,nucl_norm = DR(phi(x0,I), mu, gamma, 500, I )" ] }, { "cell_type": "code", "execution_count": null, "id": "missing-scenario", "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 282 }, "id": "missing-scenario", "outputId": "e54afb42-5dcc-4bc0-ad36-7c7467daf4b0" }, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "[]" ] }, "metadata": {}, "execution_count": 20 }, { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxdVbn/8c+TeU7azG3SZmjSNh2hoQMtpUCBAgIyCiiCUisqityrXhGvV/2JKHrxAiIWEQdQwYkZWqZCQUonOw9pk85t0qRjkraZ1++Pc0hPaIqW5JwkO9/365VXzll7Z++1yuHJyrPWXsucc4iIiDeF9XQFREQkeBTkRUQ8TEFeRMTDFORFRDxMQV5ExMMieroCgdLS0lxeXl5PV0NEpE9Zvnz5PudcemfHelWQz8vLY9myZT1dDRGRPsXMtp/smNI1IiIepiAvIuJhCvIiIh6mIC8i4mEK8iIiHqYgLyLiYQryIiIe5okgf/hoM/e8vIH566qormvo6eqIiPQavephqI9qc3Udj72zlbmtWwDIGRDLaUMGMD43hZLsJEZkJTIgPqqHaykiEnrWmzYNKS0tdR/1ideG5lbW7allxY6DrNhxiBU7DrLn8PFefVZSDMOzEhmRncjIrCSGZyVSkB5PdER4d1VfRKRHmNly51xpZ8c80ZMHiIkMZ8LQAUwYOqC9rLqugY2VdWysqvV/r2NRxX6aWtsACDMYMjCOYRkJFKYnUOj/PiwjgeTYyJ5qiohIt/FMkO9MRmIMGYkxTC8+vm5Pc2sbW/cdYUNlLRXV9ZTX1FNRfYSFm/a1B3+AtIRohmXE+4J/+y+AeAYlxxIWZj3RHBGRU+bpIN+ZyPAwijMTKc5M7FDe0trGroPHqKipp7y6noqaeipqjvDC6koOH2tuPy82MpyC9MDg73udnxZPTKRSPyLSu/S7IH8yEeFh5KXFk5cWz3kjM9vLnXPsP9JERbUv6PuCfz0rdh7k+dV7eH9Iw8w34Nse/NN9Pf/CjARS46MwU+9fREJPQf5fMDPSEqJJS4hmUkFqh2MNza1s3ecP/NXHfwEs3nKAY82t7eclx0b6An5A3r8wPZ7cgXFEhntiFquI9FIK8l0QExnOyOwkRmYndShva3NU1jb4e//17b8E3tpUw1+W72o/LzLcGJoaT2F6PMMyEtqvlZcaT7jy/iLSDRTkgyAszBicEsvglNgOg74AtQ3NbKk50uEXQHl1Pa9vqKalzZf7iY0MpzgrkZLsxPbAPyIrkcQYzfgRkVOjIB9iSTGRjM9NYXxuSofyxpZWNu+tZ0NlLRsq69hQWcvLa6v405Kd7efkDIhlzOBkxuemMC43hTGDk4mP1n9CETk5RYheIjoinNGDkxk9OLm9zDlHVW1De+Bfv6eWNbsP8/LaKsA3z784M5FxOSmMH5LCuJwUhmclKtUjIu0U5HsxMyM7OZbs5FjOHXF8xs/++kZW7zrMip2HWLXzEPPXV/HUMl+PPzE6gtK8AUwqSGVyQSqjByURocFdkX5LQb4PSk2I5pwRGZwzIgPw9fi37z/Kyp2HWLLtAIu37GdBWQ0A8VHhTMgbyOSCgUwvSqckO0kPc4n0I55Zu0Y6qq5rYMnWAyzecoDFW/ezaW894HuSd8bwdGYMT+esYekkx2kwV6Sv+7C1axTk+4maukYWbqrhzU01LNxUw+FjzYSHGROGDmDWqCxmjc5iUEpsT1dTRD4CBXnpoKW1jVW7DvFmWQ2vrNtL2d46AMblpjBrVBYXjc4iLy2+h2spIv8uBXn5UFtq6pm3rop5a6tYveswACOyEpk1OosLSrIYmZ2oZRlEejEFefm37T50jHlrq5i/toql2w/gHAxOieX8kkwuGJXJxLyBmq0j0ssoyMtHUlPXyBsb9/LKur28Xb6PppY2kmMjOXdEBueXZHJ2cboexhLpBRTkpcuONrWwcNM+XllfxRsbqzl0tJmoiDCmFqZywagsZo7MJD0xuqerKdIvKchLt2ppbWPptoO8un4vr26oYueBY5jBGUMHcnVpDpeOHURslNbWFwmVHg3yZjYLuB8IBx51zv3oZOcqyPc9zjk2VNYxf10VL6zeQ0XNEZJiIrhqQg5XnZ7D8KxELacsEmQ9FuTNLBzYBJwP7AKWAtc759Z3dr6CfN/mnGPJ1gP8YfEOXl5bSXOrIyo8jDPyB3DTlDzOL8nULB2RIOjJjbwnAuXOuS3+ijwJXA50GuSlbzMzJhWkMqkglX31JfyjfB/r9tTy4upK5jy+nCkFqcyZXsBZRWmaoSMSIsHuyV8NzHLOzfa/vxGY5Jy7LeCcOcAcgCFDhkzYvn170OojPaOltY0/LdnB/a9vZl99EwPiIpk1OosvzhhG7sC4nq6eSJ/Xkz35f8k59wjwCPjSNT1cHQmCiPAwbpySxyfOGMKCsmpeWlPJMyv28PSK3dwyLZ8xg1NIiI5gaGocEeHGoaPNDMtIUC5fpBsEO8jvBnID3uf4y6QfiooI48JRWVw4KovKw8f4/vPreWhBRafnpidG87mz8pk9rUCrZop0QbCD/FKgyMzy8QX364AbgnxP6QOyk2N5+FMTOHikid2HjlHf2EJFTT1tbY7EmEj+9s9d/PCljazZXct9145r79U3trTyp8U7eLdiP/HREVxTmsOZhWk93BqR3iuoQd4512JmtwHz8U2hfMw5ty6Y95S+ZUB8FAPiowCYXJDaXn75+EH88q0t/HjeRuobmvnWxSMZGB/F7N8vY8WOQxSkxXPoWDNPr9jNN2YN54szhvVUE0R6taDn5J1zLwEvBfs+4i1mxhdmFJIUG8F3nl3Hgp8txAyiwsN46IbTuWRsNg3NrXzjr6u5d14ZsZHhfGZqfk9XW6TX6fGBV5EP88lJQzl/ZCYvr62ioqaea0tz2/fBjYkM575rx9HY0sr3nl/P2t21ZCRFU3usmYgw45qAc0X6Ky1rIH1eY0srd7+4gb8u39W+iNrRplYaW1qZM72Q/zi/mKiIMA4eaWJBWTVHmlqZNiyNfK2ZLx6htWukX2hrc5j5Uj2HjzVzz0sbeHLpTkZkJVKQHs9rG6ppamkDIDzM+M8LivnC2YV6Clf6vF49T16kuwROtUyOjeRHV43lnBEZ3PfKJpZvP8gnSnO5tjSXlLhIfjxvI/fOKyMlNoobJg3pwVqLBJeCvHja+/PyP+j+606jtqGF7z63jjGDkxmelcgfF29n0Zb9DBkYx+yzCshMiumBGot0L6VrpN86cKSJjz3wNg3+PP7WfUcYmhrHnkPHSIqJ5PFbJlEyKAnnHAs372Pp1gMUZyVy8egsrb0jvYpy8iInsW7PYb7+l9W0Occ3Zg3n3BGZlFfXceOvl+AcPP2lM3nwjXL+uHhH+8+cWZjKozeVEhelP4Sld1CQFzlF6/Yc5uqHF3GsuRWAW88u5I7zi3hmxW6++fc1XHHaYO67dnwP11LERwOvIqdo1KBknpg9kd++u52ZIzO4fPxgAD5xxhB2HzzGA2+Uc8Vpg5lamMZTy3Yyb20VCdER3DhlaIcnd0V6mnryIqeoobmVi+5/m6NNLRSkJbBoy34K0+M5fKyFffWNfOdjJXx2mu/p250HjvLimkoiwozLxg0iQ4O5EgTqyYt0o5jIcO69eiyf+c1Slu84yHcvLeGmM/NobGnj9idX8P0X1pMYE0FiTCRf+8sq6htbALj/9c386tOl6ulLSKknL/IR1TU00+Z8c/Lf19jSyuzfLePtzfsAGJeTzIPXn05zWxuff3w51bUNvPzV6QxOie2paosHaeBVJISONLbwwBubiY+KYM70AmIiwwHYvv8IF93/tn92zhkcbWrh/tc2s2jLfoozE/nyucMYmqqlFuTUKciL9BJz36rgnpc3ctOUoby5qYbt+48yMX8g63YfJjzMeHLOlPa5+a+s38uqnYcYNSiZi0ZnafMUOSkFeZFeoqmljUsffIeyvXUMGRjHvVePZXJBKjsPHOXauYsAeGL2JH7wwnoWlNW0/9x5IzJ4+FMTiIrQQ1hyIgV5kV6krqGZpdsOMLkgtcMDVev31HLNL9/lSFMr4WHGf18ykhsmDeUPi7fzvefXM3taPt/+WEn7+a1tjpa2NqIjwnuiGdKLaHaNSC+SGBPJuSMyTygvGZTE3BtLeWrZTq6fmNu+reFnpuazubqex/6xlWvPyKUoI4Ffvb2FB14v51hzK7NGZfHDK8aQHBd5wjVFFORFepFpRWlMKzpxz9qvXzCc51fu4afzyyjKTOChBRWcOyKDvNR4nnhvO7sOHuXPt04hOiKcbfuO8LtF22hobuPqCYOZMHRg6BsivYaCvEgfMCA+is9NL+C+Vzfxyvq9XD8xl7s/PoawMGNi/kBufWI5D7y+melF6Xz2t0tpaXNEhYfx1NId/PSacVx5ek5PN0F6iIK8SB8x+6x8auoaGZoax2en5rfPtpk1OosrTxvMQwsqeGhBBcMyEnj8lokkxUTyud8v45t/X8NpQwa074S16+BR3i3fT0F6PBOGDtCmKR6ngVcRD6htaOY/nlpFZLjx/ctHk54YDUB1XQMzfvIm543M5MHrT2Pe2kq+8qeVNLX6dsi6tjSHH105VtMz+zgNvIp4XFJMJI/edOL/4xmJMXxq8lB+/c5WrjsjlzueWsWowUncc+UYnl25h4ffrKA4M5HZZxVwtKmFH7+8kQVlNRSmx3PXJSMZlpHYA62R7qRJtyIed+PkoRjwyUcXExMZxi8/NYERWUl848LhnF2czgOvb+bQ0Sa+/tfV/P697RRnJrJy5yGunfseuw8d6+nqSxcpyIt4XO7AOL59yUhGZCXyi09OaN/W0Mz45kUjqGtsYfz3X+XF1ZV87YLhPHpTKX/7wpk0NLfyP8+ua7/Oih0Hmf27pdz02BIWbqo52e2kl1GQF+kHbp6az7yvTmdKYccVMEdmJ/H56YUAXDI2my+c7XtdkJ7ArWcX8tqGvWzaW8fa3Yf5xCPvsWrXYcqr67n5N0uYv66qw7X21jawff8RetM4n2jgVaTfc86xubqeYekJHQZgDxxpYtIPX+PTU/JYtu0AVbUNvPSVs4iNCucT/lTOm1+fQWJ0BD+at5FHFm7BOZgxPJ2f33A6CdEa8guVDxt4VU9epJ8zM4ozE0+YYTMwPooZwzP49TtbWbXrMF+dWUxqQjRxURH84OOjOXCkiSeX7OCZlbuZ+9YWrjo9h69fOJyFm2q48+9rOlyrtc1x8EhTKJslfvpVKyIn9fnpBby6fi9pCdFcFfBA1bjcFCbmDeTx97YTbkZJdhL3XuWbitnU0sb9r2/mM1PzOH3IAJZvP8Btf1xB5eEGJgwdwIPXn8YgracfMurJi8hJleYN5O9fPJMXvjzthBUwrynNYeeBY2zbf5SvnFfU/pfA588uIDE6gj8u3sHBI03c8rtlREWE8dWZRWyqquOW3y2jqaWt/TpLth7goQXlLN12IKRt6y/UkxeRD3X6kAGdll86bhCLtuynsbmNC0qOL7gWFxXBrNFZvLy2itT4KA4fa+apOVMYnpXIiKwkbn1iOX9etpNPTR7K44u28d8BM3j+a9YIvjCjMNhN6lfUkxeRjyQmMpz7rh3PQ588/YR8/gWjsqhvbGHuwi2cMzyD4Vm+h6ouHJXJuJzk9kXVfvDiBs4Zns7yb8/kkrHZ3Dt/Iyt2HGy/zqGjTcxbW8mGytqQts1LFORFpNtNLji+8uWl47LbX5sZl44bxMaqOu6dV0Zzaxs/uGIMqQnR/PiqsQyIi+IXb1YAUF5dx/k/W8itT/yTi+5/m/99pSzk7fACBXkR6XaJMZF8/cLh3HbOMC4dO6jDsQtHZQHw3Ko9TB2W1r6peUJ0BNedkcvrG/ZSU9fI1/6ymrY2xxO3TOLqCTk8+EY5r63f236dN8uq+fRjS/jyn1ZQUVMfusb1MQryIhIUXzpnGF+7cDgR4R3DTO7AuPbXV0/ouATyrNFZtDl4aEE5K3ce4j8uKGZaURr3XDmGwvR4fvpKGc45FlXs5zO/XUpFdT1vllVz1cPvsmP/0ZC0q69RkBeRkPt/l49ibE4yF4/J7lA+alAyA+Ii+e2724iOCOPScb6/AiLDw/jM1Hw2VtWxdnctP3hxPUMGxvHKHdN5/rZptLY6vv/C+vbrHD7azNy3KnhoQTk1dY0hbVtvE7Qgb2bfNbPdZrbS/3VxsO4lIn3LjVPyeO62aUR+oJcfHmZMHebbGevMwlSSYo5vafixsdmYwc8XbGbdnlpmT8snPjqCvLR45kwv4LUNe9lSU09tQzNXPvwP7nl5Iz+ZX8YlD7zNroP9t5cf7J78z5xz4/1fLwX5XiLiAWMGJwNQnNVxmeOUuCiKMxKZv86Xlw/8K+DqUl/a56U1lcx9q4It+47wh9mTePEr0zjW1Mp//W11+7n/3HGQ8+97izHfnc/dL66npbUNL9M8eRHpVa6bOIRt+4+0L5wWqDRvAGV764iLCic1Ibq9PDs5lrE5ySwoq6Gipp4LSjLb/yK4fWYRP3hxA//ccZChA+O4+bElJMVGMr0onV+9vZXWNvjOpSUha1+oBbsnf5uZrTazx8ys0ycqzGyOmS0zs2U1NVq+VKS/S46N5J4rxzIwPuqEY5+dlg/A5eMHn3BsXE4Ky7cf5NDRZi4bd/z49ROHEB0RxvOr9vDrd7ZS39jCb24+g4c+eTo3Th7Kb97dSnm1b3ZOS2sbv/nHVr7x11UdZvL0ZV0K8mb2mpmt7eTrcuBhoBAYD1QC/9vZNZxzjzjnSp1zpenp6V2pjoh4XGF6AkvuOo//6aTnPXpwUvvrwHn68dERTClMZcHGal5cU8n04nSKMn2poNtnFhEZFsbji7YB8K2n1/C959fz0poqZv9+GY+9s/WE+wQuydAXdCld45yb+e+cZ2a/Al7oyr1ERMC3pWFnJuanEma+XH1gKgfgrKJ03izzZQpu8f81AJCWEM2M4em8UVbNx3cc5M/LdnHr2YX85wXFfOGJf/Kjlzdy3sgMhqbGc+hoE199aiVvltVQkp3E/103nuLM3r89YjBn1wTOjboCWBuse4mI5KfFs/77s/j5DaefcGxsTnL76xnFGR2OTS5IZeeBY8x9awtREWF8+dxhRIaHcfcVowH4zT+2AfDtZ9byj/J93HxmHtV1jdz82BIOH2sOXoO6STBz8vea2RozWw2cA9wRxHuJiBATGd5peUn28VTOkNS4Dscm5vtSO/PWVTG1MJV4/2YnmUkxnF+SyXOr9lBeXc8Lqyu59exCvnvZKB69qZQ9hxv41cIt7depPHyMR9/ewqvr9/aq3bGCNrvGOXdjsK4tInIq4qMjuGRsNmMHJ59wbGR2EokxEdQ1tHDeyMwOxy4YlcmLayq5d95GzOCGSUMAGJ+bwvklmTy5dCd3nF/MzgNHuezn71Db0AL4zvvhFWOC37B/g554FZF+4aEbTufzZ584LTM8zJiY5+vNnzuiYyqn1F/+yvq9jBmcTHby8c1OLhs3iH31jSzbdqD9aduXbz+Lz52Vzx8X7+DVgNk5W2rquf3JFdzx1Eq2hHidHQV5Een3Pje9gDtmFp+wY9XglFiyknwDvaflpnQ4dvZw32zAp1fsZkFZNTdPzWdkdhL/NWsEealxzH3Lt5rmoaNNXP+r93ht/V5eWVfFtXMXUV3bEIJW+SjIi0i/N7kgldtnFnV67JwRvmD+wSdwk2IiKUyP58mlO3EOrjjNNzc/IjyMGyYNYdn2g+w8cJTfL9rO3tpGnpwzhWe+NJXaYy38NITLJivIi4h8iG9eNJI50wu4bNygE46Ny/H17lPiIskLGNCdMdyX9lm0ZT/PrNjN5IKBjMlJpigzkWtKc3hu1R6ONPry9845nl25mzfLqoNSfwV5EZEPkRwbybcuHkliwGJp77vSv7l5Xmo8Zsd3xyrKSCA1PooXVleyZd8RZgYM6H78tME0NLexwB/Uf/5GObc/uZKnlu4MSv21do2IyEc0rSiNX35qAsWZCR3KzYySQUks3OR7AKtk0PEpnKflphAVHsaaXYc5szCNB98o5+IxWTx4/Ynz+7uDevIiIl0wa3QWBekJJ5QXBpSNyDoe5CPCwyjKTGBDVR0LNlbT1NrGrWcXEv6BfXK7i3ryIiJBkJ8WD/imaH5wsbURWUks3FzDwLhI0hOjGT3oxPn73UU9eRGRIMgd6JuO+cHdrwAGp8RQU9fIzoPHGJaeQFiQevGgnryISFCcVZTOty4ewQ2Thp5w7P1B2qrDDe3LKgSLgryISBBEhocxp5ONTwDen4izt7aBzKTOV9XsLkrXiIiEWJg/yre0ObKSov/F2V28V1CvLiIiJwjMwGclqycvIuIpAc9NkZagnryIiKcEPh0bER7cMKwgLyISYoE9+eBNnvRRkBcRCTELCO1hFtwwryAvIhJigc8+BTnGK8iLiIRasAN7IAV5EZEQU7pGRMTDTOkaERHvCpxCqZ68iIjHaOBVRMTD7CSvg0FBXkQkxALTNaZ0jYiItyhdIyLiZRp4FRHxLuXkRUQ8LEw9eRER79LDUCIiHhbCpWsU5EVEQq1DuiZM6RoREW/RpiEiIt4VGNh79cCrmV1jZuvMrM3MSj9w7E4zKzezMjO7sGvVFBHxjrAOT7wG914RXfz5tcCVwNzAQjMrAa4DRgGDgNfMrNg519rF+4mI9Hl9ZnaNc26Dc66sk0OXA0865xqdc1uBcmBiV+4lIuIVHTfy7sXpmg8xGNgZ8H6Xv+wEZjbHzJaZ2bKampogVUdEpPfoVekaM3sNyOrk0F3OuWe7WgHn3CPAIwClpaWuq9cTEelLgj3w+i+DvHNu5ke47m4gN+B9jr9MRKTf67DUcJDvFax0zXPAdWYWbWb5QBGwJEj3EhHpUwKff+rtUyivMLNdwBTgRTObD+CcWwf8GVgPzAO+pJk1IiI+FsKnobo0hdI59zTw9EmO3Q3c3ZXri4h4UZ+ZQikiIqeuz6RrRETko+j7A68iInISpp68iIh3hfJhKAV5EZEQ67DHq4K8iIi3hAVE3r66do2IiJxEYGBXT15ExGs08Coi4l1hHli7RkRETkIDryIiHtZxWQOla0REPCXYefgO9wrZnUREBAh+Hj6QgryISKiFMMoryIuIhJjSNSIiHqZ0jYiIhwV7Rk0gBXkRkRALU05eRMS7QtiRV5AXEQk9pWtERDxL6RoREQ/TwKuIiIdpCqWIiIfpYSgREQ/T7BoREekWCvIiIiGmdI2IiIcpXSMi4mHqyYuIeJh68iIiHqZ58iIiHqYnXkVEPEzpGhERD+sz6Rozu8bM1plZm5mVBpTnmdkxM1vp//pl16sqIuINoZxdE9HFn18LXAnM7eRYhXNufBevLyLiOaFM13QpyDvnNkBoBxFERPo688imIflmtsLM3jKzs4J4HxGRPqVX9eTN7DUgq5NDdznnnj3Jj1UCQ5xz+81sAvCMmY1yztV2cv05wByAIUOG/Ps1FxHpo3pVkHfOzTzVizrnGoFG/+vlZlYBFAPLOjn3EeARgNLSUneq9xIR6Wv6/Dx5M0s3s3D/6wKgCNgSjHuJiPQ1fWaPVzO7wsx2AVOAF81svv/QdGC1ma0E/grc6pw70LWqioh4QygHXrs6u+Zp4OlOyv8G/K0r1xYR8So98Soi4mEK8iIiHuaVefIiItIJ9eRFRDxMO0OJiHhYn1mFUkRETp3SNSIiHtbnn3gVEZGTU09eRMTDNPAqIuJhGngVEfEwpWtERDxM6RoREekWCvIiIiGmdI2IiIcpXSMi4mGaXSMi4mF64lVExMP6zB6vIiJy6tSTFxGRbqEgLyLiYQryIiIepiAvIuJhCvIiIh6mIC8i4mEK8iIiHqYgLyLiYQryIiIepiAvIuJhCvIiIh6mIC8i4mEK8iIiHqYgLyLiYQryIiIepiAvIuJhXQryZvYTM9toZqvN7GkzSwk4dqeZlZtZmZld2PWqiojIqepqT/5VYLRzbiywCbgTwMxKgOuAUcAs4BdmFt7Fe4mIyCnqUpB3zr3inGvxv30PyPG/vhx40jnX6JzbCpQDE7tyLxEROXUR3XitzwJP+V8Pxhf037fLX3YCM5sDzAEYMmRIN1ZHRKT3+tGVYyjKTAj6ff5lkDez14CsTg7d5Zx71n/OXUAL8IdTrYBz7hHgEYDS0lJ3qj8vItIXXTcxNJ3afxnknXMzP+y4md0MfAw4zzn3fpDeDeQGnJbjLxMRkRDq6uyaWcA3gMucc0cDDj0HXGdm0WaWDxQBS7pyLxEROXVdzcn/HIgGXjUzgPecc7c659aZ2Z+B9fjSOF9yzrV28V4iInKKuhTknXPDPuTY3cDdXbm+iIh0jZ54FRHxMAV5EREPU5AXEfEwBXkREQ+z41Pbe56Z1QDbu3CJNGBfN1Wnr1Cb+we1uX/4qG0e6pxL7+xArwryXWVmy5xzpT1dj1BSm/sHtbl/CEabla4REfEwBXkREQ/zWpB/pKcr0APU5v5Bbe4fur3NnsrJi4hIR17ryYuISAAFeRERD/NEkDezWf4Nw8vN7Js9XZ/uYmaPmVm1ma0NKBtoZq+a2Wb/9wH+cjOzB/z/BqvN7PSeq/lHZ2a5ZrbAzNab2Tozu91f7tl2m1mMmS0xs1X+Nn/PX55vZov9bXvKzKL85dH+9+X+43k9Wf+uMLNwM1thZi/433u6zWa2zczWmNlKM1vmLwvqZ7vPB3n/BuEPARcBJcD1/o3EveC3+DZCD/RN4HXnXBHwuv89+Npf5P+aAzwcojp2txbgP51zJcBk4Ev+/55ebncjcK5zbhwwHphlZpOBHwM/86/2ehC4xX/+LcBBf/nP/Of1VbcDGwLe94c2n+OcGx8wHz64n23nXJ/+AqYA8wPe3wnc2dP16sb25QFrA96XAdn+19lAmf/1XOD6zs7ry1/As8D5/aXdQBzwT2ASvicfI/zl7Z9zYD4wxf86wn+e9XTdP0Jbc/xB7VzgBcD6QZu3AWkfKAvqZ7vP9+TxbRC+M+D9STcN94hM51yl/3UVkOl/7bl/B/+f5KcBi/F4u/1pi5VANfAqUAEccs61+E8JbFd7m/3HDwOpoa1xt/g/fJd6slwAAAHmSURBVDvLtfnfp+L9NjvgFTNbbmZz/GVB/Wx3dWco6UHOOWdmnpwDa2YJwN+Arzrnav07jwHebLfz7Zw23sxSgKeBET1cpaAys48B1c655WY2o6frE0LTnHO7zSwD3456GwMPBuOz7YWefH/bNHyvmWUD+L9X+8s98+9gZpH4AvwfnHN/9xd7vt0AzrlDwAJ8qYoUM3u/IxbYrvY2+48nA/tDXNWumgpcZmbbgCfxpWzux9ttxjm32/+9Gt8v84kE+bPthSC/FCjyj8pHAdfh20jcq54DbvK/vglfzvr98k/7R+QnA4cD/gTsM8zXZf81sME5d1/AIc+228zS/T14zCwW3xjEBnzB/mr/aR9s8/v/FlcDbzh/0ravcM7d6ZzLcc7l4ft/9g3n3CfxcJvNLN7MEt9/DVwArCXYn+2eHojopsGMi4FN+PKYd/V0fbqxXX8CKoFmfPm4W/DlIV8HNgOvAQP95xq+WUYVwBqgtKfr/xHbPA1f3nI1sNL/dbGX2w2MBVb427wW+I6/vABYApQDfwGi/eUx/vfl/uMFPd2GLrZ/BvCC19vsb9sq/9e692NVsD/bWtZARMTDvJCuERGRk1CQFxHxMAV5EREPU5AXEfEwBXkREQ9TkBcR8TAFeRERD/v/j89/rfpC7+gAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ], "source": [ "plt.plot(np.log10( nucl_norm-np.min(nucl_norm)+1e-20 ))" ] }, { "cell_type": "markdown", "source": [ "Test with varying number $p$ of measurements." ], "metadata": { "id": "0Hp7W0IabuFC" }, "id": "0Hp7W0IabuFC" }, { "cell_type": "code", "source": [ "r = 4\n", "x0 = np.random.randn(n,r) @ np.random.randn(r,n)\n", "# we varie p\n", "n_test = 10 # number of tests\n", "p_list = np.round( np.linspace(.02,.5,n_test)*n**2 )\n", "err = []\n", "n_iter = 500\n", "import progressbar\n", "for i in progressbar.progressbar(range(n_test)):\n", " p = int( p_list[i] )\n", " I = indices(n,p)\n", " y = phi(x0,I)\n", " x,nucl_norm = DR(y, mu, gamma, n_iter, I )\n", " err.append( np.linalg.norm(x-x0) )" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "2IsPW9-YzjJ_", "outputId": "7dfd8c35-6dca-4d42-fa11-95933071d043" }, "id": "2IsPW9-YzjJ_", "execution_count": null, "outputs": [ { "output_type": "stream", "name": "stderr", "text": [ "100% (10 of 10) |########################| Elapsed Time: 0:00:27 Time: 0:00:27\n" ] } ] }, { "cell_type": "code", "source": [ "plt.plot(p_list/(n**2), err/np.linalg.norm(x0));" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 282 }, "id": "J7kQcVFPIQiV", "outputId": "e099cce5-fe34-48dc-fa70-12e222a5f52f" }, "id": "J7kQcVFPIQiV", "execution_count": null, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "[]" ] }, "metadata": {}, "execution_count": 25 }, { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAdDUlEQVR4nO3deXxU9b3/8dcnk4QQSIJA2JJAgAQhCFaJgBetFrTiArTVWqlasSjaCrbX/nprf7a2pe2jrd7bWxdcqLhXrdW2F/cq7lSW4MK+hH0RCFsIS8j2vX9k6I0YyEBm5jvL+/l45PHIzDnMeX8fg2+/nDPzPeacQ0RE4l+K7wAiIhIeKnQRkQShQhcRSRAqdBGRBKFCFxFJEKm+Dty5c2dXWFjo6/AiInFpwYIFO5xzuc1t81bohYWFlJWV+Tq8iEhcMrP1R9umUy4iIglChS4ikiBU6CIiCaLFQjezh81su5ktPsp2M7O7zazczBaa2enhjykiIi0JZYb+KDD6GNsvBIqDP5OA+1sfS0REjleLhe6cexfYdYxdxgGPu0ZzgA5m1j1cAUVEJDThOIeeB2xs8nhT8LnPMbNJZlZmZmUVFRVhOLSIiBwW1YuizrnpzrlS51xpbm6zn4tv0bod+7nj1eXU1TeEOZ2ISHwLR6FvBgqaPM4PPhcRry3Zyn1vr+aaR+axe39NpA4jIhJ3wlHoM4FvBT/tMhyodM59GobXbdYN5/TljksHM3/tbsZNm82KrVWROpSISFwJ5WOLTwMfACeb2SYzm2hmN5rZjcFdXgbWAOXAH4HvRixt0OVnFPD0pOEcrK3nq/fN5tXFEfv/h4hI3DBft6ArLS11rV3LZWtlNTc8uYBPNu7he6OK+d6oYlJSLEwJRURij5ktcM6VNrctrr8p2i0ngz9PGs6lp+dz16xV3PjkAvYdqvMdS0TEi7gudICMtAD/+fXB/PSSEmYt387X7pvN+p37fccSEYm6uC90ADNj4lm9eezaoWzbe4ix987mvVX6nLuIJJeEKPTDziruzMzJI+iWncE1D8/joffW4OsagYhItCVUoQP06tSOv3733zi/pCu/emkZP3j2E6pr633HEhGJuIQrdIB2bVK5/8oh/Pt5/fjrR5v5xoMfsLWy2ncsEZGISshCB0hJMb53XjEPXj2E8u37GHPv+yxYf6w1xkRE4lvCFvphFwzsxt9uGkFmeoArps/hmXkbfEcSEYmIhC90gH5ds/ifm0YwvE8nbv3rIm7/n8XUanEvEUkwSVHoAB0y03lkwhlcf3ZvHv9gPVfPmMvOfYd8xxIRCZukKXSA1EAKt11cwu8vP5UPN+xh7L2zWbKl0ncsEZGwSKpCP+xrp+fzlxvOpL7Bcdn9H/Diwi2+I4mItFpSFjrAqQUdmDllBCU9spn81Efc+dpyGhr0JSQRiV9JW+gAXbIyeOr6YXyjtIBpb63m+sfL2Ftd6zuWiMgJSepCB2iTGuC3lw5i6riBvLOygq9Mm82ain2+Y4mIHLekL3RoXNzrW2cW8sTEYew5UMu4abN5a8V237FERI6LCr2JM/t2YubkEeSflMm3H53P/W+v1uJeIhI3VOhHyD8pk+e/cyYXDerO715dzs3PfMzBGi3uJSKxT4XejMz0VO4dfxo/vOBkXly4hcse+Ceb9xz0HUtE5JhU6EdhZtz0pSJmXFPKhp0HGHvP+8xds9N3LBGRo1Kht2Bk/6787aYR5LRN48qH5vLknPW+I4mINEuFHoKiLu35200jOKu4Mz/5+2L+/98WUVOnxb1EJLao0EOU0zaNGdecwXfO7ctTczdw5UNzqKjS4l4iEjtU6MchkGL8aHR/7h5/Gos2V3L1jLnUa7kAEYkRKvQTMPbUHtxx2aks31rFK4s/9R1HRARQoZ+wiwd1p6hLe+6ZVa5FvUQkJqjQT1AgxZgysogV26r4x9KtvuOIiKjQW+OSwT3o07kdd80q1xIBIuKdCr0VAimNXz5a9ule3limxbxExC8VeiuN+0IPenbM5O5ZqzRLFxGvVOitlBpIYfKXili0uZK3V1T4jiMiSSykQjez0Wa2wszKzezWZrb3NLO3zOwjM1toZheFP2rs+urpeeR1aMtdmqWLiEctFrqZBYBpwIVACTDezEqO2O0nwLPOudOAK4D7wh00lqUFUrjpS0V8vHEP763a4TuOiCSpUGboQ4Fy59wa51wN8Aww7oh9HJAd/D0H2BK+iPHh0iF59MjJ0CxdRLwJpdDzgI1NHm8KPtfUz4GrzGwT8DIwpbkXMrNJZlZmZmUVFYl1vrlNaoDvnNuXBet388FqLbMrItEXroui44FHnXP5wEXAE2b2udd2zk13zpU650pzc3PDdOjY8fXSArpmt+GuWat8RxGRJBRKoW8GCpo8zg8+19RE4FkA59wHQAbQORwB40lGWoAbz+nL3LW7mKObYYhIlIVS6POBYjPrbWbpNF70nHnEPhuAUQBmNoDGQk+scyohGj+0J53bt+GeNzVLF5HoarHQnXN1wGTgNWAZjZ9mWWJmU81sbHC3HwDXm9knwNPABJekVwYbZ+l9mF2+k7J1u3zHEZEkYr56t7S01JWVlXk5dqQdqKnj7N+9xcC8HB7/9lDfcUQkgZjZAudcaXPb9E3RCMhMT+X6L/bh3ZUVfLRht+84IpIkVOgRcvXwXpyUmcY9b5b7jiIiSUKFHiHt2qRy3dl9eHP5dhZtqvQdR0SSgAo9gr51Zi+yM1K5W594EZEoUKFHUFZGGhPP6sPrS7exZItm6SISWSr0CJswopCsNqncq3PpIhJhKvQIy2mbxrUjCnll8VZWbK3yHUdEEpgKPQq+fVZv2qUH9O1REYkoFXoUdMhM55p/K+SlRZ9Svl2zdBGJDBV6lFx3dh/apgV0Ll1EIkaFHiUd26Vz9fBezPxkC2sq9vmOIyIJSIUeRded3Yf01BSmvbXadxQRSUAq9CjKzWrDlcN68fePN7N+537fcUQkwajQo+yGL/YhkGLcp1m6iISZCj3KumRn8M2hPXn+w01s3HXAdxwRSSAqdA9uOKcPKWbc/45m6SISPip0D7rntOXyM/L5S9lGtuw56DuOiCQIFbon3zm3CIAHNEsXkTBRoXuS16Etlw3J55l5G9laWe07jogkABW6R989t4h653jwXc3SRaT1VOgeFXTM5Gun5fHU3A1sr9IsXURaR4Xu2U1fKqK2voE/vrvGdxQRiXMqdM8KO7fjK1/I48k5G9ix75DvOCISx1ToMeCmkUVU19Xz0HtrfUcRkTimQo8BfXPbM2ZwDx7/YB279tf4jiMicUqFHiMmjyziYG09D7+vWbqInBgVeozo1zWLi07pzqP/XEflgVrfcUQkDqnQY8jkkUXsO1THw7M1SxeR46dCjyEDumdzwcCuPDx7LXurNUsXkeOjQo8xU0YWU1Vdx2Oz1/mOIiJxRoUeY07Jy+G8AV146P217DtU5zuOiMSRkArdzEab2QozKzezW4+yz+VmttTMlpjZU+GNmVymjCym8mAtj3+wzncUEYkjLRa6mQWAacCFQAkw3sxKjtinGPgxMMI5NxD4fgSyJo1TCzpw7sm5PPTeWvZrli4iIQplhj4UKHfOrXHO1QDPAOOO2Od6YJpzbjeAc257eGMmnykji9m1v4Y/zV3vO4qIxIlQCj0P2Njk8abgc031A/qZ2Wwzm2Nmo5t7ITObZGZlZlZWUVFxYomTxJBeJ3F2cWemv7uGgzX1vuOISBwI10XRVKAYOBcYD/zRzDocuZNzbrpzrtQ5V5qbmxumQyeum0cVs2NfDU/N2+A7iojEgVAKfTNQ0ORxfvC5pjYBM51ztc65tcBKGgteWuGMwo6c2acTD7yzmupazdJF5NhCKfT5QLGZ9TazdOAKYOYR+/ydxtk5ZtaZxlMwWuA7DG4eVUxF1SH+PH9jyzuLSFJrsdCdc3XAZOA1YBnwrHNuiZlNNbOxwd1eA3aa2VLgLeCHzrmdkQqdTIb36cjQwo7c//ZqDtVpli4iR2fOOS8HLi0tdWVlZV6OHW/eX7WDq2bM5VdfOYWrhvfyHUdEPDKzBc650ua26ZuicWBEUSdO79mB+99eTU1dg+84IhKjVOhxwMy4eVQxm/cc5K8fbvIdR0RilAo9TpzTL5dT83OY9nY5tfWapYvI56nQ48ThWfrGXQf5+0dHfmpURESFHldG9u/CwB7ZTHurnDrN0kXkCCr0OHJ4lr5u5wFeWLjFdxwRiTEq9Dhz/oCu9O+WxT1vllPf4OcjpyISm1TocSYlpXGWvqZiPy8t+tR3HBGJISr0ODR6YDeKu7TnnlmraNAsXUSCVOhxKCXFmDKqmFXb9/Hqkq2+44hIjFChx6mLB3WnT2477tYsXUSCVOhxKpBiTBlZxPKtVby+bJvvOCISA1TocWzM4B4Udsrkv/6xQp9LFxEVejxLDaRw64X9WbltH0/rrkYiSU+FHucuGNiNM/t04r9eX8meAzW+44iIRyr0OGdm3D6mhL0Ha/nDG6t8xxERj1ToCWBA92y+OawnT8xZz6ptVb7jiIgnKvQEccv5J9MuPcDUF5fi6y5UIuKXCj1BdGyXzvfP68d7q3Ywa9l233FExAMVegK5+sxe9M1tx69eWqobSoskIRV6AkkLpHD7mIGs23mAx/65znccEYkyFXqCOadfLiP7d+HuWeVUVB3yHUdEokiFnoB+cvEAqmvr+c/XVviOIiJRpEJPQH1y23PtiEKeXbCRxZsrfccRkShRoSeoKaOK6ZiZzi9eWKKPMYokCRV6gsrOSOP/XXAy89ft5sWFurORSDJQoSewy0sLKOmezW9eXsbBGn2MUSTRqdATWCDF+NmYErZUVvPgu6t9xxGRCFOhJ7hhfTpx8eDuPPDOarbsOeg7johEkAo9Cfz4wv44B799ZbnvKCISQSr0JJB/UiY3fLEPMz/Zwvx1u3zHEZEICanQzWy0ma0ws3Izu/UY+11qZs7MSsMXUcLhxnP70i07g6kvLNVNpUUSVIuFbmYBYBpwIVACjDezkmb2ywK+B8wNd0hpvcz0VH58UX8Wba7kuQ83+Y4jIhEQygx9KFDunFvjnKsBngHGNbPfL4HfAdVhzCdhNPbUHpzeswN3vLqCqupa33FEJMxCKfQ8YGOTx5uCz/2LmZ0OFDjnXjrWC5nZJDMrM7OyioqK4w4rrWNm/GzMQHbsO8S9b5X7jiMiYdbqi6JmlgL8HvhBS/s656Y750qdc6W5ubmtPbScgFMLOnDZkHwefn8ta3fs9x1HRMIolELfDBQ0eZwffO6wLOAU4G0zWwcMB2bqwmjs+o8LTiY9kMKvX1rmO4qIhFEohT4fKDaz3maWDlwBzDy80TlX6Zzr7JwrdM4VAnOAsc65sogkllbrkp3BTSOLeGPZNt5bpVNfIomixUJ3ztUBk4HXgGXAs865JWY21czGRjqgRMa3R/SmZ8dMpr6wlLr6Bt9xRCQMQjqH7px72TnXzznX1zn36+BztzvnZjaz77mance+jLQAt108gFXb9/GnuRt8xxGRMNA3RZPYl0u6MqKoE79/fSW799f4jiMiraRCT2Jmxk8vKaGqupY/vLHSdxwRaSUVepLr3y2bK4f14sm5G1ixtcp3HBFpBRW6cMv5/WjfJpWpL+p2dSLxTIUunNQunX8/r5jZ5Tt5fek233FE5ASp0AWAK4f3oqhLe3798jIO1el2dSLxSIUuAKQFUvjpJSWs33mAR2av8x1HRE6ACl3+5Zx+uZw3oAv3zFrF9iotmikSb1To8hm3XVxCTX0Dd766wncUETlOKnT5jN6d23HtiN489+EmFm7a4zuOiBwHFbp8zuSRRXRql84vXliqjzGKxBEVunxOdkYaP7zgZBas383MT7b4jiMiIVKhS7MuG1LAKXnZ/PaV5RyoqfMdR0RCoEKXZgVSjNsvGcinldU8+M4a33FEJAQqdDmqob07csng7jzwzmo27znoO46ItECFLsf044sGAPCbl3W7OpFYp0KXY8rr0JYbz+nLiws/Zd7aXb7jiMgxqNClRTee05fuORn84oUl1DfoY4wisUqFLi1qmx7g1gv7s2TLXp5bsNF3HBE5ChW6hGTsqT0o7XUSd762gr3Vtb7jiEgzVOgSEjPjZ2MGsnN/Dfe+We47jog0Q4UuIRuUn8Nlp+fzyOy1rN2x33ccETmCCl2Oyw9Hn0x6IIVfv7TUdxQROYIKXY5Ll6wMpowq5o1l23lnZYXvOCLShApdjtu1Iwrp1SmTX764lNr6Bt9xRCRIhS7HrU1qgNsuGkD59n38ac5633FEJEiFLifk/JKunFXUmd+/vpJd+2t8xxERVOhygsyMn15Swv6aev779ZW+44gIKnRphZO7ZXHVsJ78ae56lm/d6zuOSNJToUurfP+8fmRlpDFVt6sT8U6FLq1yUrt0bjm/H/9cvVO3qxPxLKRCN7PRZrbCzMrN7NZmtt9iZkvNbKGZzTKzXuGPKrHqymE9Oa1nB3743EL+uXqH7zgiSavFQjezADANuBAoAcabWckRu30ElDrnBgPPAXeEO6jErtRACo9MOIPCTplc/1gZH2/c4zuSSFIKZYY+FCh3zq1xztUAzwDjmu7gnHvLOXcg+HAOkB/emBLrOmSm88TEYXRq34YJj8xj5bYq35FEkk4ohZ4HNF0Ee1PwuaOZCLzS3AYzm2RmZWZWVlGhr40nmq7ZGTw5cRjpgRSunjGXjbsOtPyHRCRswnpR1MyuAkqBO5vb7pyb7pwrdc6V5ubmhvPQEiN6dsrkiYnDqK5t4MqH5rJ9b7XvSCJJI5RC3wwUNHmcH3zuM8zsPOA2YKxz7lB44kk8OrlbFo9eewY79h3i6hnz2HNA3yQViYZQCn0+UGxmvc0sHbgCmNl0BzM7DXiQxjLfHv6YEm9O63kSf/xWKWt37GfCI/PZf6jOdySRhNdioTvn6oDJwGvAMuBZ59wSM5tqZmODu90JtAf+YmYfm9nMo7ycJJERRZ2555unsWhzJZOeKKO6tt53JJGEZr6+3VdaWurKysq8HFui6/kFm/jBXz7hyyVdue/K00kN6PtsIifKzBY450qb26b/siTiLh2Sz8/GlPCPpdv40fOLaGjQEgEikZDqO4Akh2tH9KbyYC1/eGMV2W1Tuf2SEszMdyyRhKJCl6j53qhiKg/W8sjsdeS0TeP75/XzHUkkoajQJWrMjJ9eXEJVdV3jTD0jjW+f1dt3LJGEoUKXqEpJMX77tUFUVdcy9cWlZLdN47IhWilCJBx0UVSiLjWQwl1XnMaIok786PmFvLZkq+9IIglBhS5eZKQFmH51KYPycpjy1EfMLteyuyKtpUIXb9q1SeXRa8+gd+d2XP94GR9t2O07kkhcU6GLV43L7g6lc/s2THhkPiu2atldkROlQhfvumRn8KfrhtEmtXHZ3Q07teyuyIlQoUtMKOiYyZPXDaOmvoErZ8xhm5bdFTluKnSJGf26ZvHotUPZta+Gq2fMZfd+LbsrcjxU6BJTvlDQgT9+q5R1Ow8w4ZF57NOyuyIhU6FLzPm3os7cO/40Fm/Zy/WPadldkVCp0CUmfXlgN+64dDAfrNnJlKc/oq6+wXckkZinQpeYdemQfH4+poTXl27jP55bqGV3RVqgtVwkpk0Y0ZvKg3X89xsryW6bxs/GaNldkaNRoUvMu3lUEZUHa3l49lqy26Zxy/ladlekOSp0iXlmxk8uHkBVdS13z1pFTts0JmrZXZHPUaFLXEhJMX7ztUFUVdfxyxeXkpWRyuWlBb5jicQUXRSVuJEaSOGu8V/g7OLO3Pr8Ql5d/KnvSCIxRYUucaVNaoAHrhrCqQUduPnpj3l/lZbdFTlMhS5xp12bVB6dMJQ+ue2Y9EQZH2rZXRFAhS5xKiczjccnDiU3qw0THp7H8q17fUcS8U6FLnGrS1YGT04cRtv0AFfPmMe6Hft9RxLxSoUuca2gYyZPThxGXX0DV82Yy9ZKLbsryUuFLnGvOLjs7u79NVw1Yy7PL9jEym1VWv9Fko4552d9jNLSUldWVubl2JKYPli9kxufXEDlwVoAMtJSKOmezaC8HE4J/hR3aU9qQPMYiV9mtsA5V9rsNhW6JJL6Bsfqin0s3lzJos2VLN5cyZItezlQ07gEb5vUFAZ0z+aUvMaiH9gjh35ds0hPVclLfFChS1Krb3Cs3bH/cyV/+OYZ6YEU+nfPYmCPHAblNf7069aeNqkBz8lFPk+FLnKEhgbHup37Wbxlb2PRb6pk8ZZKqqobSz4tYPTrmvWZ0zX9u2WRkaaSF79aXehmNhq4CwgADznnfnvE9jbA48AQYCfwDefcumO9pgpdYo1zjg27DgRn8Xv/NaM/fE4+NcUo7prFKT2yGZTfWPIDumXTNl0lL9HTqkI3swCwEjgf2ATMB8Y755Y22ee7wGDn3I1mdgXwVefcN471uip0iQfOOTbtPvh/p2uCM/pdwRtYB1KMotz2wVl843n5kh7ZZKZr3TuJjGMVeih/64YC5c65NcEXewYYByxtss844OfB358D7jUzc77O54iEiZlR0DGTgo6ZXDioO9BY8lsqq1kcPB+/aHMl76zczvMfbgr+GSg4KTNmLrTqdiCx5+ZRxYw5tUfYXzeUQs8DNjZ5vAkYdrR9nHN1ZlYJdAI+s3KSmU0CJgH07NnzBCOL+GVm5HVoS16HtlwwsBvQWPLb9h7610XX1RX7iIXpjCMGQsjn5LRNi8jrRvXfhc656cB0aDzlEs1ji0SSmdEtJ4NuORmcX9LVdxxJUqH8m3Az0PROAvnB55rdx8xSgRwaL46KiEiUhFLo84FiM+ttZunAFcDMI/aZCVwT/P0y4E2dPxcRia4WT7kEz4lPBl6j8WOLDzvnlpjZVKDMOTcTmAE8YWblwC4aS19ERKIopHPozrmXgZePeO72Jr9XA18PbzQRETkesfG5KhERaTUVuohIglChi4gkCBW6iEiC8LbaoplVAOuDDztzxLdKk4jGnrySefzJPHZo3fh7Oedym9vgrdA/E8Ks7GiLzSQ6jT05xw7JPf5kHjtEbvw65SIikiBU6CIiCSJWCn267wAeaezJK5nHn8xjhwiNPybOoYuISOvFygxdRERaSYUuIpIgolboZjbazFaYWbmZ3drM9jZm9ufg9rlmVhitbNEQwvi/aGYfmlmdmV3mI2OkhDD2W8xsqZktNLNZZtbLR85ICWH8N5rZIjP72MzeN7MSHzkjoaWxN9nvUjNzZpYwH2UM4X2fYGYVwff9YzO7rtUHdc5F/IfGZXdXA32AdOAToOSIfb4LPBD8/Qrgz9HIFkPjLwQGA48Dl/nOHOWxfwnIDP7+nSR877Ob/D4WeNV37miNPbhfFvAuMAco9Z07iu/7BODecB43WjP0f91o2jlXAxy+0XRT44DHgr8/B4wys0S5v22L43fOrXPOLQQafASMoFDG/pZz7kDw4Rwa74qVKEIZ/94mD9tBwtwINJT/7gF+CfwOqI5muAgLdexhFa1Cb+5G03lH28c5VwccvtF0Ighl/InqeMc+EXgloomiK6Txm9lNZrYauAO4OUrZIq3FsZvZ6UCBc+6laAaLglD/3l8aPNX4nJkVNLP9uOiiqMQMM7sKKAXu9J0l2pxz05xzfYEfAT/xnScazCwF+D3wA99ZPHkBKHTODQZe5//OUJywaBV6st9oOpTxJ6qQxm5m5wG3AWOdc4eilC0ajve9fwb4SkQTRU9LY88CTgHeNrN1wHBgZoJcGG3xfXfO7Wzyd/0hYEhrDxqtQk/2G02HMv5E1eLYzew04EEay3y7h4yRFMr4i5s8vBhYFcV8kXTMsTvnKp1znZ1zhc65Qhqvn4x1zpX5iRtWobzv3Zs8HAssa/VRo3jV9yJgJY1Xfm8LPjeVxjcQIAP4C1AOzAP6+L5SHeXxn0Hjebb9NP7LZInvzFEc+xvANuDj4M9M35mjPP67gCXBsb8FDPSdOVpjP2Lft0mQT7mE+L7/Jvi+fxJ83/u39pj66r+ISILQRVERkQShQhcRSRAqdBGRBKFCFxFJECp0EZEEoUIXEUkQKnQRkQTxv3aAhT9iYavEAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] } ], "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.9.1" }, "colab": { "name": "Trace Norme Douglas-Rachford", "provenance": [] } }, "nbformat": 4, "nbformat_minor": 5 }