Example_Optimize_GenBlockCoheherence.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Example_Optimize_GenBlockCoheherence.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "w45mRgCsCdgS"
      },
      "source": [
        "import cvxpy as cp\n",
        "import numpy as np\n",
        "from scipy.linalg import toeplitz, circulant\n",
        "import matplotlib.pyplot as plt"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "amuaG12qOq7A"
      },
      "source": [
        "def R(D, n, d):\n",
        "    n_y, n_x  = D.shape\n",
        "    R = np.zeros((n,n))\n",
        "    for k in range(0,n):\n",
        "      for l in range(k,n):\n",
        "          I = np.zeros((n_x,n_x))\n",
        "          I_s = np.zeros((n,n))\n",
        "          I[l*d:(l+1)*d, k*d:(k+1)*d] = np.ones(d)\n",
        "          I_s[l,k] = 1 \n",
        "          if k==l:\n",
        "            R = R+cp.norm(cp.multiply(I,D), 2)*I_s\n",
        "          else:\n",
        "            R = R+cp.norm(cp.multiply(I,D), 2)*I_s+cp.norm(cp.multiply(I.T,D), 2)*I_s.T\n",
        "        \n",
        "    return R"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fWRsuCMXZljB"
      },
      "source": [
        "# Test for a random Matrix D"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UxEJ7t3-JXx2"
      },
      "source": [
        "d = 2 # number of block length\n",
        "m = 5 \n",
        "n = 10 # number of blocks\n",
        "\n",
        "# creating a random matrix with orthonormal blocks\n",
        "D = np.zeros((m*d, n*d))\n",
        "for l in range(0,n):\n",
        "    D[:,l*d:(l+1)*d] = np.random.random(size=(m*d, d))\n",
        "    D[:,l*d:(l+1)*d], _ = np.linalg.qr(D[:,l*d:(l+1)*d])\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2R36jZTBWSsz",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "e69ac23b-bcd7-41e6-8f69-b2ca98e159d9"
      },
      "source": [
        "\n",
        "# Construct the problem.\n",
        "B_cvx = cp.Variable((m*d,n*d))\n",
        "\n",
        "I = np.kron(np.eye(n), np.ones((d, d)))\n",
        "k = np.tile(np.eye(d), (1, n)).T\n",
        "b = cp.multiply(B_cvx.T@D, I)\n",
        "b = cp.matmul(b, k) #extracting the diagonal blocks of D^TB\n",
        "constraints = [b == k]\n",
        "\n",
        "objective = cp.Minimize(1/d*cp.max(R(B_cvx.T@D-np.eye(n*d), n, d)))\n",
        "\n",
        "#objective = cp.Minimize(cp.max(cp.abs(A.T@x-np.eye(n))))\n",
        "\n",
        "prob = cp.Problem(objective, constraints)\n",
        "\n",
        "# The optimal objective value is returned by `prob.solve()`.\n",
        "result = prob.solve()\n",
        "# The optimal value for x is stored in `x.value`.\n",
        "#print(x.value)\n",
        "# The optimal Lagrange multiplier for a constraint is stored in\n",
        "# `constraint.dual_value`.\n",
        "#print(constraints[0].dual_value)\n",
        "print(result)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "0.46081424546188576\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "J4SFQXJoaNge"
      },
      "source": [
        "first block of $B^TD$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HCU2jSzZD2c0",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "941d6367-a887-43ce-c223-cc0ab6e223da"
      },
      "source": [
        "# Construct the problem.\n",
        "B_up = cp.Variable((m*d,n*d))\n",
        "\n",
        "I = np.kron(np.eye(n), np.ones((d, d)))\n",
        "k = np.tile(np.eye(d), (1, n)).T\n",
        "b = cp.multiply(B_up.T@D, I)\n",
        "b = cp.matmul(b, k) #extracting the diagonal blocks of D^TB\n",
        "constraints = [b == k]\n",
        "\n",
        "objective = cp.Minimize(1/d*(cp.norm(B_up.T@D, 'fro')))\n",
        "\n",
        "#objective = cp.Minimize(cp.max(cp.abs(A.T@x-np.eye(n))))\n",
        "\n",
        "prob = cp.Problem(objective, constraints)\n",
        "\n",
        "# The optimal objective value is returned by `prob.solve()`.\n",
        "result = prob.solve()\n",
        "# The optimal value for x is stored in `x.value`.\n",
        "#print(x.value)\n",
        "# The optimal Lagrange multiplier for a constraint is stored in\n",
        "# `constraint.dual_value`.\n",
        "#print(constraints[0].dual_value)\n",
        "print(result)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "3.6164415773707432\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TN8tWw_aMkPm"
      },
      "source": [
        "# Comparison of \"direct\" solution and frobenius estimation"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G5cA7HyqtBis"
      },
      "source": [
        "generalized coherence:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DETiKU70tHje",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "cf1c8476-9298-4619-81e6-f641eefd0971"
      },
      "source": [
        "1/d*np.max(R(B_up.T@D-np.eye(n*d), n, d).value)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0.5662568328312224"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 6
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "oVDu7wWlMvR5",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "cf293d74-fc18-4193-ab60-aff015245d5c"
      },
      "source": [
        "1/d*np.max(R(B_cvx.T@D-np.eye(n*d), n, d).value)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0.4608047108265468"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 7
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zk_wQ-B0MjHj"
      },
      "source": [
        "images"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "uAsQDh-6D-A3",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 266
        },
        "outputId": "d61ff9b5-9934-47f2-f0b0-5014dbf08c57"
      },
      "source": [
        "plt.imshow(B_up.value)\n",
        "plt.colorbar()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.colorbar.Colorbar at 0x7f7fe8f0f810>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 8
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWIAAADnCAYAAAAkVlylAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAZX0lEQVR4nO3dfbRddX3n8ffHm5vnJyBAQkCgGC1oC+iVBx8o8qCQukA7SmOrBZVJ6cCIrna6oMxIm8502elUp7JYaBRGdBzEoVDTmhJ5EK12iIRnQkBCBJMQEpJgEkhI7sN3/jg79HA559792+ecu89OPq+19sp52N/9+2Wfe773d3/7t38/RQRmZlaeN5RdATOz/Z0TsZlZyZyIzcxK5kRsZlYyJ2Izs5KNK7sCZmat+MD7psSWrYO59r3/kd3LIuKcDlcpmROxmVXalq2D/GzZG3Pt2zPnqVkdrk4hTsRmVmkBDDFUdjVa4kRsZpUWBP2Rr2uiWzkRm1nluUVsZlaiIBis+FQNTsRmVnlDOBGbmZUmgEEnYjOzcrlFbGZWogD63UdsZlaeINw1YWZWqoDBaudhJ2Izq7banXXV5tnXzKzixGDObdQjSTdI2iTpsSbvS9KXJa2W9Iikt7fjf+BEbGaVVrtYp1xbDt8ARpqd7VxgXrYtBK5rtf7gRGxmFVcbR9yeFnFE/BjYOsIu5wPfjJp7gZmS5rT6f3AfsZlV3lC+1i7ALEkr6p4vjojFCUXNBdbWPV+XvbYh4Riv40RsZpW2t0Wc0+aI6OtgdQpxIjazSgvE4Nj1sq4Hjqh7fnj2WkvcR2xmlTcUyrW1wRLgD7LRE6cA2yKipW4JcIvYzCouEHuipy3HknQTcDq1vuR1wNVAL0BEfAVYCswHVgM7gU+2o1wnYjOrtNoNHe354z4iPjbK+wFc2pbC6jgRm1nlJVys60pOxGZWaRFiMKp9ucuJ2Mwqb8gtYjOz8tQu1lU7lVW79ma232vnxbqyOBGbWeUNtmeMcGmciM2s0sb4zrqOcCI2s8ob8qgJM7Py1Cb9cSI2MytNIPrbdItzWZyIzazSIvANHWZm5ZJv6DAzK1PgFrGZWel8sc7MrERB2yZ9L40TsZlVWgD9nmvCzKxM8nzEZmZlCnxnnZlZ6dwiNjMrUYTcIm5k4syJMXXO1KSYl/rHd6IqDfX2DBWK2/PK2PzeUn+x3+5vmDKQHBMvFfs/HXjgjuSYLb+alhzzhv7kEAB6dkdyzMS5uwqVFQWu2O94eVJ6QW9I/z8B9PSm/7zHjmK3DA9NSqvjwOYXGdzxckvN2drFOt/i/DpT50zlt288Lylm+bojO1GVhubM3F4obs3q2W2uSWOT1hb7WCaevCU5Zs9PDypU1oLfvzs55ltL3pccM/m5Yt/RmWvSM/ibFz1WqKw9Q+mf149+dlxyTEws1oCYeWj6L82BnxxYqKyXjt2TtP/zf3FNoXJey2vWmZmVqnaxzn3EZmalqvqddblqL+kcSU9KWi3pik5Xyswsr7131uXZutWoLWJJPcC1wNnAOuA+SUsi4vFOV87MLI/9YfHQk4DVEbEGQNJ3gPMBJ2IzK10E9A/t+4l4LrC27vk64OThO0laCCwEmDJ7SlsqZ2Y2mlrXRLUTcdtqHxGLI6IvIvomzpzYrsOamY1qMJtvYrStW+VJxOuBI+qeH569ZmZWur3D19p1sW60wQmSLpL0gqSHsu3iVv8Pebom7gPmSTqaWgJeAPxeqwWbmbVH+7omEgYn3BwRl7WlUHIk4ogYkHQZsAzoAW6IiJXtqoCZWavauGZdKYMTct3QERFLgaWdrIiZWRG1URO555qYJWlF3fPFEbG47nmuwQnAv5N0GvBz4HMRsbbBPrl15M66XQO9rNySNi/Drk2TC5XVszP9T5K1e9InnwEYf8zL6UE/Tx9BsvugYnMK/PrMrckxDx49o1BZN33njOSY8QUm8OndWWyim61v6U2OeXjz3EJlHT8r/ZLJpOfTJ6nZfVCxVt/AU+nzRgxOKFQUvJL4fSz2o/4aiUslbY6IvhaL/EfgpojYLekPgRuB9C9EnWqP+TAzo9Y1kWfLYdTBCRGxJSJ2Z0+/Dryj1fo7EZtZpbV51MSrgxMkjac2OGFJ/Q6S5tQ9PQ9Y1er/wZP+mFnltWvURLPBCZIWASsiYgnwGUnnAQPAVuCiVst1IjazSosQA228s67R4ISI+Hzd4yuBK9tWIE7EZrYP6OaZ1fJwIjazSvPE8GZmXcCJ2MysRInjiLuSE7GZVV4bb3EuhROxmVVaBAzsBxPDm5l1NXdNmJmVyH3EZmZdIJyIXy8Q/QNps0tNXlesKr/14QeSY57fVWz2tZUb5oy+0zBDx+xKjln+3muTYwDO+pv/lB507EChsvqnp8+K1rst/csy4xPrkmMAntl4UHLM7AmvFCrr0b89Pjlm9zvTpx075L7kEAA+8p9vT4657pHTCpV1yW/8NGn/a6/ZUaic4XyxzsysRBHuIzYzK5kYrPioiVFrL+kIST+U9LiklZIuH4uKmZnlFaFcW7fK0yIeAP44Ih6QNA24X9IdDRbTMzMbc/vFXBMRsQHYkD3eIWkVtXWdnIjNrHxR6yeusqQ+YklHAScCyztRGTOzIvabUROSpgJ/D3w2IrY3eH8hsBCg9+DpbaugmdlIYn+4WAcgqZdaEv52RNzaaJ+IWBwRfRHRN25G+srFZmZFReTbutWoLWJJAq4HVkXEFztfJTOzNN08IiKPPC3idwOfAM6Q9FC2ze9wvczMcqm1dvfx4WsR8ROoeE+4me3T9vnha2Zm3a6b+3/z6EgiHvccHLoobdKfI697sFBZT20/ODlm7ZaZhco6ZtHu5Jh1/y39au47lxW7efFHf/Lfk2PeOG5qobLm3XNRcsz0A9IneNm0o1j9Ji+fnBzzyprDCpW17W3pn/HBb92YHHPngpuSYwDeevt/SI45bFna93evr24+I2n/F3Y8WqiceoEYqvioCbeIzazyKt4gdiI2s4qL6o+acCI2s+qreJPYidjMKs8tYjOzEgUwNOREbGZWngAq3iKu9pgPMzPaO9eEpHMkPSlptaQrGrw/QdLN2fvLs1kpW+JEbGbVFzm3UUjqAa4FzgWOAz4m6bhhu30aeDEi3gR8CfjrVqvvRGxmFZdvnomcF/ROAlZHxJqI2AN8Bzh/2D7nAzdmj28BzswmRyvMidjMqi9/i3iWpBV128JhR5oLrK17vi57reE+ETEAbAMOaqX6vlhnZtUWEPlHTWyOiL5OVqcIt4jNbB+gnNuo1gNH1D0/PHut4T6SxgEzgC3F6+5EbGb7gjZdrAPuA+ZJOlrSeGABsGTYPkuAC7PHHwHujmht/reOdE0Mzg1+9ZdpM5XN7J9UqKw5k7clx6zbWmz2tWcXpZ+ugUfT1++b3bcpOQbggiv/JDnmxWOLXWOISek/d89vH58cc/l7f5AcA/B3G89Ojhk6/eVCZY2/K/3naev9hyTHzL/m0uQYgJOvejo5ZuAtxdpoGx48plBcy9p0i3NEDEi6DFgG9AA3RMRKSYuAFRGxhNqKRd+StBrYSi1Zt8R9xGZWbW2+oSMilgJLh732+brHrwAfbVuBOBGb2T6g6hPD5/77Q1KPpAcl/VMnK2RmlmxI+bYuldIivhxYBaR3epqZdZD2hxaxpMOB3wa+3tnqmJklyjtioouTdd4W8f8E/hSY1myH7A6VhQDjD3Gj2czGivb92dckfRDYFBH3j7RfRCyOiL6I6Bs3I33hRjOzwvaDFvG7gfMkzQcmAtMl/e+I+Hhnq2ZmltNQ2RVozagt4oi4MiIOj4ijqA1cvttJ2My6xt5xxHm2LuVxxGZWeVUfNZGUiCPiHuCejtTEzKyoiidiT/pjZlayzkz6s3McWx85OClGxxf7lbbx6VnJMdN/3lOorPd8YsSBIw3dc987kmOeX3tgcgxA7wd3JcdodbERLpOfS/8dfvC5zyXHfGXle5NjgEL9gb09g4WKemVGesy5596XHPO9g09MLwiYvGdicswvfnRUobIOWZ121eyFYvMsvc5+1TVhZtZ1gq6+fTkPJ2Izqz63iM3MyuWuCTOzsjkRm5mVzInYzKw8CndNmJmVz6MmzMzK5RaxmVnZnIjNzErkPmIzsy7gRGxmVi5VfGL4jiTi8etf5ugr/l9SzNqr3lWorFnr038VTv699MlnAJ669M3JMbP+a3pZe1bMSY4B6P1Fb3LM9f/+mkJlffaJ302Oefbhw5JjTjrlyeQYgClH7kmOufup9M8XYN6ZzybHHDPxhfSCCs6VuPGmI5Njeue/WKisXS8ekLT/kJuCgFvEZrYvcNeEmVmJfLHOzKwLVDwR5+p1kjRT0i2SnpC0StKpna6YmVlukXPrUnlbxH8H3B4RH5E0Hii2rIOZWZuJsRk1IelA4GbgKOAZ4IKIeN1VTUmDwKPZ019GxHmjHXvUFrGkGcBpwPUAEbEnIn6Vt/JmZh0V/zbxz2hbi64A7oqIecBd2fNGdkXECdk2ahKGfF0TRwMvAP9L0oOSvi5pyvCdJC2UtELSin525ynbzKw9xqZr4nzgxuzxjcCHWj5iJk8iHge8HbguIk4EXqbBb4KIWBwRfRHR18uEdtXPzGx0+RPxrL0NxmxbmFDKoRGxIXv8PHBok/0mZse+V1KuZJ2nj3gdsC4ilmfPb6F5k9zMbMwldDtsjoi+pseR7gRmN3jrqvonERFS01KPjIj1kn4NuFvSoxHx9EiVGjURR8TzktZKektEPAmcCTw+WpyZ2Zhp04iIiDir2XuSNkqaExEbJM0BNjU5xvrs3zWS7gFOBEZMxHlvmvyPwLclPQKcAPxVzjgzs86K2qiJPFuLlgAXZo8vBL43fAdJB0iakD2eBbybHA3XXMPXIuIhoGlz3sysVGMzRvgLwHclfRp4FrgAQFIfcElEXAwcC3xV0hC1hu4XIqI9idjMrJuNxS3OEbGFWtfs8NdXABdnj/8V+I3UY3ckEce0yQyc/I6kmCPO+GWhslY/cnhyzPihYtNYvfiuackxf3H0bckxn1u3IDkGYP4HHkqOuXfXMYXKmjYhfYjitH9Ij3ngiPTPF2DShP7kmA8d+3Chss6esTI55rLvfTI5ZsIRLyfHAEzY/rrRpqMaKFQSbD82LXJoYrs6d9tzmLK4RWxm1dblty/n4URsZpUmPPuamVnpnIjNzMrmRGxmVjInYjOzEnmFDjOzLuBEbGZWrrGYGL6TnIjNrPLcNWFmVibf0GFm1gWciM3MyuM765o44I3b+dCX70yK+eK9Zxcq69ivbU2OWf35qYXKuuhTdyXHPPnKnOSYeV9Ln7AGYOnlxyXHvOnqnYXKeuZTByfHfO0bX02Oufi2lJVs/s3Lc3clx9y6vNhMr3c/fkpyzITp6eWceOra9CDgkQXpP4P3vfObhcq6ZO0ZSft/f/IrhcoZTkPVzsRuEZtZtbmP2MysfO6aMDMrW8UTca4Z0iV9TtJKSY9JuknSxE5XzMwsL0W+rVuNmoglzQU+A/RFxNuAHqDYEhJmZp0QObculbdrYhwwSVI/MBl4rnNVMjNLENW/xXnUFnFErAf+B/BLYAOwLSJ+MHw/SQslrZC04qWtxYZfmZml2juOeF/vmjgAOB84GjgMmCLp48P3i4jFEdEXEX1TD+xtf03NzJqJyLd1qTwX684CfhERL0REP3Ar8K7OVsvMLL99vkVMrUviFEmTJQk4E1jV2WqZmeWU90JdFyfiUS/WRcRySbcADwADwIPA4k5XzMwsr6pfrMs1aiIirgau7nBdzMwK2S8SsZlZ1wq6+kJcHh1JxFt2T+Fbz5yUFKOeYidyzZ9PSI6ZvmxKobK+wcnJMf07xyfHTD49PQZg9gEbkmOe/vjsQmWNf9P25JiL770wOWZwarGmzpSfpX/G/b9ebNjluHM2J8dM600v66cPvTk5BqBnZ64baF/jtFs/U6isHeftSNr/pf5iP+vDdfOFuDzSPyEzs24zBhfrJH00m+phSFLTOVMlnSPpSUmrJV2R59hOxGZWaWN4Q8djwO8AP25aF6kHuBY4FzgO+JikUScKdx+xmVVbxJhMDB8RqwBqo3ibOglYHRFrsn2/Q+2GuMdHCnKL2Myqr3vGEc8F6pdSWZe9NiK3iM2s8hK6HWZJWlH3fHFEvHpfhKQ7gUZXsK+KiO8Vr+HInIjNrNoCyN81sTkiml5oi4izWqzNeuCIuueHZ6+NyF0TZlZ93dM1cR8wT9LRksZTm7t9yWhBTsRmVnljMWpC0oclrQNOBb4vaVn2+mGSlgJExABwGbCM2pw8342IlaMd210TZlZ5YzRq4jbgtgavPwfMr3u+FFiacmwnYjOrti6fWS0PJ2Izq7TaDR3VzsROxGZWfRWffU3Rgd8kkl4Anm3w1iwgfYaU9uuGenRDHaA76tENdYDuqEc31AHGrh5HRsTBrRxg+vTD4519l+ba9+4f/tn9Iw1fK0tHWsTNTqykFd1wErqhHt1Qh26pRzfUoVvq0Q116KZ65OI+YjOzso3NXBOd5ERsZtXni3VJumWtu26oRzfUAbqjHt1QB+iOenRDHaB76jG6qP5SSR25WGdmNlamT50bJx//R7n2vfNf/8v+c7HOzGxMVbw96URsZpWnoWr3TXRk0p/R1mySNEHSzdn7yyUd1YE6HCHph5Iez9aZurzBPqdL2ibpoWz7fAfq8YykR7Pjr2jwviR9OTsXj0h6ewfq8Ja6/+NDkrZL+uywfdp+LiTdIGmTpMfqXjtQ0h2Snsr+PaBJ7IXZPk9JSl91dPR6/I2kJ7JzfpukmU1iR/z8WqzDn0taX3fO5zeJTV4DLbEeN9fV4RlJDzWJbcu5aLugdkNHnq1LtT0R51yz6dPAixHxJuBLwF+3ux7AAPDHEXEccApwaZO1o/4lIk7ItkUdqAfA+7LjN+qbOheYl20LgevaXXhEPLn3/wi8A9hJg8lLaP+5+AZwzrDXrgDuioh5wF3Z89eQdCBwNXAytaVnrm6WsFuoxx3A2yLiN4GfA1eOED/S59dKHQC+VHfOXzdRTNE10FLqERG/W/fz8ffArSPEt+NctJUIFPm2btWJFvGrazZFxB5g75pN9c4Hbswe3wKcqVEWgkoVERsi4oHs8Q5qU9KNumRJCc4Hvhk19wIzJc3pYHlnAk9HRKM7H9sqIn4MbB32cv1nfyPwoQahHwDuiIitEfEitaTZKIkVrkdE/CCbshDgXmoTeHdMk3ORR57vU1vqkX0HLwBuKnr80kTk27pUJxJxnjWbXt0n+zJsAw7qQF0AyLo+TgSWN3j7VEkPS/pnSW/tQPEB/EDS/ZIWNni/0BpXLVhA8y9ap88FwKERsSF7/DxwaIN9xvqcfAr45ybvjfb5teqyrHvkhiat/rE8F+8FNkbEU03e7/S5KM6JuLtJmkrtz63PRsT2YW8/QO1e9+OBa4B/6EAV3hMRb6f2p+Wlkk7rQBm5qLZiwHnA/23w9lici9eI2tjJUr8dkq6i1o317Sa7dPLzuw44BjgB2AD8bRuPXcTHGLk13DU/y6/hPuKG8qzZ9Oo+ksYBM4At7a6IpF5qSfjbEfG6fq+I2B4RL2WPlwK9kma1sw4RsT77dxO1ftmThu1SaI2rgs4FHoiIjQ3q2fFzkdm4t+sl+3dTg33G5JxIugj4IPD70WRAfY7Pr7CI2BgRgxExBHytybHH6lyMA34HuLnZPp08F63S0FCurVt1IhHnWbNpCbD3SvhHgLubfRGKyvq7rgdWRcQXm+wze2/ftKSTqJ2Ptv1CkDRF0rS9j4H3A48N220J8AfZ6IlTgG11f7q3W9MWT6fPRZ36z/5CoNHKuMuA90s6IPtz/f3Za20j6RzgT4HzImJnk33yfH6t1KH+WsCHmxy70BpoBZwFPBER6xq92elz0Zqc3RJd3DXR9nHEETEgae+aTT3ADRGxUtIiYEVELKGWIL8laTW1CwcL2l0P4N3AJ4BH64bj/BnwxqyeX6H2S+CPJA0Au4AFbf6FcChwW5bfxgH/JyJul3RJXR2WUltmZTW10QyfbGP5r8q+PGcDf1j3Wn092n4uJN0EnE5tCfN11EZCfAH4rqRPU5sq9YJs3z7gkoi4OCK2SvpLakkIYFFEFLnQNVI9rgQmAHdkn8+9EXGJpMOAr0fEfJp8fm2sw+mSTqD2x/UzZJ9NfR2afZ+K1KFZPSLiehpcO+jUuWi7oKuTbB6+xdnMKm3GpDlx6q99Kte+yx7/K9/ibGbWCd08RjgPJ2Izqz4nYjOzEkXAYPeOiMjDidjMqs8tYjOzkjkRm5mVKACvWWdmVqaAcB+xmVl5Al+sMzMrXcX7iPf52dfMbD8wBnNNSPqoaqv9DGW35DfbL3klE7eIzazixmxCn8eozVD31Rz7vi8iNuc9sBOxmVVbAGMwxWVErAJo82JCgLsmzGxf0F3TYCavZOIWsZlVXNItzrOG9dsujojFe59IuhOY3SDuqohoNHd2I++JiPWSDqE2zeoT2VqBTTkRm1m1BUT+ccSbR5oGMyLOark6dSuZSNq7ksmIidhdE2ZWfUORb+uwoiuZOBGbWfWNzfC1D2ermpwKfF/Ssuz1wyQtzXY7FPiJpIeBnwHfz7OSibsmzKzaIsZq1MRt1BZNHf76c9SWOyMi1gDHpx7bidjMqq/id9Y5EZtZxQUxOFh2JVriRGxm1eZpMM3MuoCnwTQzK08A4RaxmVmJwhPDm5mVruoX6xQVH/ZhZvs3SbcDs3LuvjkizulkfYpwIjYzK5lvcTYzK5kTsZlZyZyIzcxK5kRsZlYyJ2Izs5L9fw6MK8TELXCUAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 2 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ySNcyi9JMy8r",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 272
        },
        "outputId": "42a3d52c-9559-42d6-92d3-18c6a78999c5"
      },
      "source": [
        "plt.imshow(B_cvx.value)\n",
        "plt.colorbar()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.colorbar.Colorbar at 0x7f7fe8ea0fd0>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 9
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWIAAADtCAYAAACFzX/DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAZkklEQVR4nO3dfZxdVX3v8c83kwfyAEkgMeSBR81VsRWDKaBoi/Jg4LZQWqHh1hYt3qhXrFbvpXB9KX3Rh5e9tSJeLTZFruClCAWpuRqFgCJUX1AGCCQQKOExEwIhDySEhCQz87t/nD28DsOZmb32nDP77Mz3/XrtV86Zs397rdnnnF/WrL3W2ooIzMysPGPKroCZ2WjnRGxmVjInYjOzkjkRm5mVzInYzKxkTsRmZiUbW3YFzMyG40MfmBybt/Tk2ve+h3bfEhGLWlylZE7EZlZpm7b0cM8t83LtO272EzNaXJ1CnIjNrOKCnugtuxLD4kRsZpUWQC/VniHsRGxmlRYEeyNfH3G7ciI2s8pzi9jMrEQB9DgRm5mVyy1iM7MSBdBT8eV8nYjNrNKCYK9bxGZmJQroqXYe9loTZlZttXHE+bahSLpK0kZJqwd4XZK+IWmtpIckHdOM38GJ2MwqTvTk3HL4LjDYWhSnAfOzbQlwxbCrjxOxmVVcAL2RbxvyWBF3AlsG2eVM4JqouRuYJmn2cH8H9xGbWeXlbO02w1xgXd3zruxnG4ZzUCdiM6u0APZG7j/uZ0jqrHu+NCKWNr9WaZyIzazSajPrcreIN0XEwmEUtx44pO75vOxnw+I+YjOrtED0MCbX1gTLgD/ORk8cD2yLiGF1S4BbxGa2D+iN5vQRS7oOOJFaF0YXcAkwDiAivg0sB04H1gI7gY81o1wnYjOrtMSuicGPFXHuEK8H8OmmFFbHidjMKk705L9Y15aciM2s0gLYS0fZ1RgWJ2Izq7QIt4jNzErXO3ITOlrCidjMKq12sc4tYjOzErlrwsysVLUpzr5YZ2ZWmr6ZdVXmRGxmldfrrgkzs/L4Yp2ZWckC0dOktSbK4kRsZpXX6xaxmVl5IuRRE2ZmZQrwOGIzs7L5Yp2ZWYkCNW1h+LI4EZtZ5blFbGZWosATOszMShV41ISZWemadc+6sjgRm1mlRchdE43sN22/mDJ7SlLMju0TC5WliT3JMeM60mMA9uwYnxwT4yI5puOVYv+79x7QmxwzZkexD/CkA3clxxR5j9WdHALAuF3p5333jPQYgFmTXk6O2bh1anJMFPy2qsjHPf2jBEBqD0H3li30vFLwA1/H44gbmDJ7Cr9z9e8kxfzy1l8vVFbHO7Ynx8yZlh4DsO6X85Jjds/dmxxz0N3jkmMAXjl1R3LMhLv2L1TWb/zhg8kxd/306OSYCVuTQwCY+eCryTFPnV8sEX/+mNuSY755Q9r3A2D3zGINiLHb05PU2B3FcuOeA9My+Pqvfb1QOfUC3yrJzKxUgdjb64t1Zmalqvo44ly1l7RI0mOS1kq6qNWVMjPLq29mXZ6tXQ3ZIpbUAXwLOAXoAu6VtCwiHml15czM8hgNy2AeC6yNiCcBJH0fOBNwIjaz0kUwKhaGnwusq3veBRzXfydJS4AlAJMPntyUypmZ5dHO3Q55NO1iXUQsBZYCzHh7wQGZZmaJRssU5/XAIXXP52U/MzMrXW3Rn32/RXwvMF/SEdQS8GLgv7S0VmZmuVV/ivOQtY+IbuAC4BZgDXBDRDzc6oqZmeXVi3JteQw1XFfSRyW9KGlltn18uPXP1UccEcuB5cMtzMys2Zo5aiJhuO71EXFBUwqlRTPrduyZwF1PvTkpZuKWYifylScPSI5Ztys9BmD3rPQVaCY/nr5Q0MuHFbvWedy8Z5Jj7lowv1BZD/zTO5NjJk5IL2fi5mKrz6w7Kb0wbShUFLdveltyTM/49Pd48jPFLkhN6Uo/h5vS314Axm9J6yIotCBRP4Hobt4U51KG61a7Y8XMjKSuiRmSOuu2Jf0O1Wi47twGRf6+pIck3SjpkAavJ/FaE2ZWaYmjJjZFxMJhFvn/gOsiYrekTwBXAx8czgHdIjazyuuNMbm2HIYcrhsRmyNid/b0SuDdw62/E7GZVVvOBX9ytppfG64raTy14brL6neQNLvu6RnURpMNi7smzKzSmrkwfER0S+obrtsBXBURD0u6FOiMiGXAn0o6A+gGtgAfHW65TsRmVmkBdPc274/7RsN1I+LLdY8vBi5uWoE4EZvZPmA0THE2M2tbfQvDV5kTsZlVnm8eamZWpnDXhJlZqUbLMphmZm2rttZEtadEOBGbWeWFW8RvFHvH0P3ifkkxrx5UbMWx979/dXLM1HG7CpW17JfpMxlfOTR9xbZVZ3wjOQbg1P/+ueSYcUcX+wDvmZoeN2ZvejkHfjJ9RTmAzb88PDlm0vyXCpW16e+PSI6ZemB6OZOfL3ACgWv+8bLkmN9akf5ZAjjlpLTv40037SxUTn++WGdmVqLwxTozs/Lt810T2Vqb1wCzqF2gXBoRl7e6YmZm+YieUXCxrhv4QkTcL2l/4D5JKxrcOsTMbMSNiuFrEbEB2JA9flnSGmor1jsRm1n5otZPXGVJfcSSDgcWAPe0ojJmZkWMmlETkqYANwGfi4jtDV5fAiwB6Jg+rWkVNDMbTDAKLtYBSBpHLQlfGxE/aLRPRCwFlgJMOPSQiv+hYGbVMQpWX5Mk4DvAmoj4WuurZGaWpre32ok4z5iPE4A/Aj4oaWW2nd7iepmZ5RJR65rIs7WrPKMm/g0q3hNuZvu0fb5rwsys3Y2q4Wt57fdiN2+7YktSzOHXrCtU1l1dRybH7HzqgEJlveWmV5NjjrjsseSYhVf9WXIMwFf+8nvJMadOTHuf+nz7pbclx/zDTz6UHLOre1xyDMDkrvSYSfcX+1zsnZzeGjvmUyuTY/7i4BXJMQDv/dHnk2OO+EFPobIeP3xm0v67e5qTgtq52yEPt4jNrNKC9u7/zcOJ2MyqzauvmZm1AfcRm5mVy10TZmYl86gJM7MSjZq1JszM2lZAVHyKsxOxmVWfuybMzMrkccRmZuWreIu42nfcMzNr8uprkhZJekzSWkkXNXh9gqTrs9fvye5cNCxOxGZWfZFzG4KkDuBbwGnAUcC5ko7qt9v5wNaIeAtwGfC3w62+E7GZVV8o3za0Y4G1EfFkROwBvg+c2W+fM4Grs8c3AidlN9AorCV9xHvmimf+anxSzPlTi90UesueSckxO6ZvLVTWf8xJW1kKYOP3FiTH7FmYvsobwOUXLE6O+dLRae9THxVYnKv36PTf66Q3pa9eB3DDuHnJMc/9VqGimL46/Tu4+qtHJ8ec//ChyTEA/+3625Njbn1r/0ZgPk88Oidp/92vFvv8vUHz+ojnAvVLQXYBxw20T0R0S9oGHARsKlqoL9aZWbUFeVu7ADMkddY9X5rdb7NUTsRmVnkJU5w3RcTCQV5fDxxS93xe9rNG+3RJGgtMBTbnrkEDufuIJXVIekDSj4ZToJlZ0zXpYh1wLzBf0hGSxgOLgWX99lkGnJc9/jDws4jhrXaR0iL+LLAGKHYbAzOzVmnShI6sz/cC4BagA7gqIh6WdCnQGRHLqN3V/nuS1gJbqCXrYcmViCXNA/4z8NdA+n1XzMxaJUC9TTxcxHJgeb+ffbnu8avA2c0rMX+L+OvAhcD+zSzczGz4cg9Na1tD9hFL+m1gY0TcN8R+SyR1Surs3r6zaRU0MxtS8/qIS5HnYt0JwBmSnqY2uPmDkv5v/50iYmlELIyIhWMPSB/ba2ZW2L6eiCPi4oiYFxGHU+uU/llEfKTlNTMzy6viidjjiM2s2gI0mhaGj4g7gDtaUhMzs6LauLWbhxf9MTMrWUu6Jsa81MGkH6XN+/jmAR8oVNa6lWmLjAActKpQUfz5xTcnx1yx4qzkmIlTdifHAOz4zN7kmL33H1SorIkb0/8UPP3tDyfHXPmLE5NjAMYctys55vT/VGzhqV90vTs55oxP/yI55rud702OAdj9/DuSY164LX3RJIAj70lb2GnL1uY0ZVXxFrH7iM2s+io+jtiJ2Myqrc1HROThRGxmldfMKc5lcCI2s+pzi9jMrGROxGZm5VF41ISZWfk8asLMrFy+WGdmVjZ3TZiZlch9xGZmbcCJ2MysZE7EDQ66Yy8z73o+KWbdQXMLlTVjQ3ov/c6ztxUq6/qPL0qO2fyRnuSYN/1rsRtlj9+Rfi6+c9llhcr67ON/kBxz+4/TF8fpmFjsG/aOdz6XHLNkZvpCPAA/PvLo5JgHXjokOebIwzYmxwBs/+f079Yrv9FdqKxnp01I2n/PY80Z7eCuCTOzsjkRm5mVyBfrzMzaQMUTca47dEiaJulGSY9KWiPpPa2umJlZbqPk5qGXAz+NiA9LGg9MamGdzMxyE6Oga0LSVOA3gY8CRMQeYE9rq2VmlqDiiThP18QRwIvA/5H0gKQrJU3uv5OkJZI6JXXu6Um/X5iZWSFRW2siz9au8iTiscAxwBURsQB4Bbio/04RsTQiFkbEwvEdE5tcTTOzQVS8jzhPIu4CuiLinuz5jdQSs5lZW+hbk3iorV0NmYgj4nlgnaS3Zj86CSh233Ezs1aoeIs476iJzwDXZiMmngQ+1roqmZklaPMkm0euRBwRK4GFLa6LmVkhI3EhTtKBwPXA4cDTwDkRsbXBfj3AquzpsxFxxlDHzjWhw8ysnY1QH/FFwO0RMR+4nQaDFjK7IuJd2TZkEoYWTXHumTSObQvelBSz4y3FVnt6dWZHckzvy/sVKmvbkelxJx/zUHLM3U+kr+YFMP1jXckxv9r55kJljR2T3gQ58spnkmPWfbPYSnSrnp2THHPtAccXKuu/vjd91bZrrzspOWb/9xdbfW3Si+krAGpvsVXRuqekfS4i/es7wIGadJzBnQmcmD2+GrgD+PNmHNgtYjOrtrwX6mrJekbffIdsW5JQ0qyI2JA9fh6YNcB++2XHvlvS7+Y5sBf9MbNKU7bltCkiBrzeJek24OAGL32x/klEhDRgZ8dhEbFe0pHAzyStiognBquUE7GZVV+TuiYi4uSBXpP0gqTZEbFB0mygYV9RRKzP/n1S0h3AAmDQROyuCTOrvBGa4rwMOC97fB7wwzfUQ5ouaUL2eAZwAjnmXTgRm1n1jcyEjq8Ap0h6HDg5e46khZKuzPZ5O9Ap6UHg58BXImLIROyuCTOrthGavhwRm6nNLO7/807g49njXwG/nnpsJ2Izq77RMLPOzKydtfOCPnk4EZtZ9TkRm5mVKNp70fc8nIjNrPrcIjYzK8+ouHloEQfP3cz/+Jtrk2JW75pXqKxfnXJocswz/zCzUFnnXvTz5JhxSl9wZdWG5NEvAKx9Pv33uvmDxc7FU1+dmxxz+R3fTY75/H3nJMcA7PdI+u26lk85qlBZvf8+LTlm17z0z8XkKLYQz7Nnpf/d/tSiK4feqYEF9y5O2n/j+PTz0JATsZlZuRTVzsROxGZWbb5YZ2bWBqrdIM631oSkP5P0sKTVkq6TVGxldTOzFtjn7+IsaS7wp8DCiPg1oANI65E3M2ulUXIX57HAREl7gUnAc62rkplZgjZv7eYxZIs4W+T4q8CzwAZgW0Tc2n8/SUv6bj+ybUux+8+ZmRVS8RZxnq6J6dRumncEMAeYLOkj/feLiKURsTAiFk490NcAzWxkCFBv5NraVZ6LdScDT0XEixGxF/gB8N7WVsvMLL99/mIdtS6J4yVNkiRqCyOvaW21zMxySruLc1sasg8hIu6RdCNwP9ANPAAsbXXFzMzyGhUTOiLiEuCSFtfFzKyYNm7t5uGramZWee3c/5tHSxJx187pXNj5+0kx3RvTV8sCGPv1XckxM/6lWFmXn3JycszMg7clx2x9Z7FVti5c8IZRhUP6X393ZqGypr9tc3LMF+4/Ozlm74vF3qs98/ckx+jZAwqVddAJm5Jj5kxM/9w+8eic5BiAiRs6kmPe95lPFCpr5zlpv1dvb7HP+usEbT0iIg+3iM2s+qqdh52IzazavDC8mVnZImpbhTkRm1nluUVsZlayUTGO2MysbQXgURNmZiWrdh52Ijaz6nMfsZlZ2TxqwsysXFVvEee6eaiZWbtSjMzC8JLOzm6i3Ctp4SD7LZL0mKS1ki7Kc2wnYjOrvt6c2/CsBn4PuHOgHSR1AN8CTgOOAs6VdNRQB25J18TuJ5/b9MTiLz3T4KUZQPoKKU32RNF6XJde1rMDv9T0c/GpC4tE3VmoHk8VKWpgbfG5oGA9Gn3QR7oOLVCsHjcmRxyWHNGARqCPOCLWANTujzGgY4G1EfFktu/3qd1q7pHBglqSiCNiZqOfS+qMiAGb9COlHerRDnVol3q0Qx3apR7tUId2qkcuaXffmCGps+750oho5o0u5gLr6p53AccNFeSLdWZWcUlrTWwa7D8YSbcBBzd46YsR8cMitcvDidjMKq9Z6xFHRPqi46+3Hjik7vm87GeDGumLde1yr7t2qEc71AHaox7tUAdoj3q0Qx2gfeoxtKitNZFnGwH3AvMlHSFpPLAYWDZUkKLiA6HNbHQ7YMrcOO7oT+Xa97Zffem+on3fks4C/jcwE3gJWBkRH5I0B7gyIk7P9jsd+DrQAVwVEX891LHdNWFm1TcC7cmIuBm4ucHPnwNOr3u+HFiecmwnYjOrvJEYvtZKLekjHmpmiaQJkq7PXr9H0uEtqMMhkn4u6ZFsNsxnG+xzoqRtklZm25dbUI+nJa3Kjt/Z4HVJ+kZ2Lh6SdEwL6vDWut9xpaTtkj7Xb5+mnwtJV0naKGl13c8OlLRC0uPZv9MHiD0v2+dxSee1oB5/J+nR7JzfLGnaALGDvn/DrMNfSFpfd85PHyA2eaZWYj2ur6vD05JWDhDblHPREn136Rhqa1NNT8Q5Z5acD2yNiLcAlwF/2+x6AN3AFyLiKOB44NMDzHC5KyLelW2XtqAeAB/Ijt+ob+o0YH62LQGuaHbhEfFY3+8IvBvYSYM/sWj+ufgusKjfzy4Cbo+I+cDt2fPXkXQgcAm18ZfHApcMlLCHUY8VwK9FxDuB/wAuHiR+sPdvOHUAuKzunL/hz9miM7VS6hERf1D3+bgJ+MEg8c04F02lCNSTb2tXrWgRvzazJCL2AH0zS+qdCVydPb4ROElDTFdJFREbIuL+7PHLwBpqg63bzZnANVFzNzBN0uwWlncS8ERENHlC2BtFxJ3Aln4/rn/vrwZ+t0Hoh4AVEbElIrZSS5qNkljhekTErRHRnT29m9owo5YZ4Fzkkef71JR6ZN/Bcyg0h7RkbhG/QaOZJf0T4Gv7ZF+GbcBBLagLAFnXxwLgngYvv0fSg5J+IukdLSg+gFsl3SdpSYPX85yvZlrMwF+0Vp8LgFkRsSF7/Dwwq8E+I31O/gT4yQCvDfX+DdcFWffIVQO0+kfyXLwfeCEiHh/g9Vafi+KciNubpCnU/tz6XERs7/fy/cBhEXE0tWEp/9qCKrwvIo6h9qflpyX9ZgvKyCUb13gG8C8NXh6Jc/E6URs7Weq3Q9IXqXVjXTvALq18/64A3gy8C9gA/H0Tj13EuQzeGm6bz/LrBCO16E/LtCIR55lZ8to+ksYCU4HNza6IpHHUkvC1EfGGfq+I2B4RO7LHy4FxkmY0sw4RsT77dyO1ftlj++1SaCZOQacB90fECw3q2fJzkXmhr+sl+3djg31G5JxI+ijw28AfxgAD6nO8f4VFxAsR0RMRvcA/DXDskToXY6mtLHb9QPu08lwMlyJybe2qFYk4z8ySZUDflfAPAz8b6ItQVNbf9R1gTUR8bYB9Du7rm5Z0LLXz0bT/ECRNlrR/32PgVGpL6dVbBvxxNnrieGBb3Z/uzTZgi6fV56JO/Xt/HtBo/v4twKmSpmd/rp+a/axpJC0CLgTOiIidA+yT5/0bTh3qrwWcNcCxC83UKuBk4NGI6Gr0YqvPxbBVvGui6eOII6Jb0gXUvjh9M0selnQp0BkRy6glyO9JWkvtwsHiZtcDOAH4I2BV3XCc/wkcmtXz29T+E/iUpG5gF7C4yf8hzAJuzvLbWOCfI+Knkj5ZV4fl1AaDr6U2muFjTSz/NdmX5xTgE3U/q69H08+FpOuAE6mteNVFbSTEV4AbJJ1PbQXJc7J9FwKfjIiPR8QWSX9JLQkBXBoRRS50DVaPi4EJwIrs/bk7Ij6p18+Savj+NbEOJ0p6F7U/rp8me2/q6zDQ96lIHQaqR0R8hwbXDlp1LpouAnrbuN8hB09xNrNKmzpxdrznyD/Jte8tj/xN4SnOreSZdWZWee3c/5uHE7GZVZ8TsZlZiQJo0nrEZXEiNrOKq/7FOidiM6s+d02YmZXIXRNmZmULCHdNmJmVy10TZmYlcteEmVkb8KgJM7MytfeCPnk4EZtZtQVuEZuZlc4tYjOzkjkRm5mVKTxqwsysVAHR01N2LYbFidjMqs9dE2ZmJdoHbpXkRGxm1ecWsZlZucItYjOzEkVAT7UT8ZiyK2BmNmzRm28bBklnS3pYUq+kAe8ELelpSaskrZTUmefYbhGbWaUFECMzjng18HvAP+bY9wMRsSnvgZ2IzazaYmQWho+INQCSmn5sd02YWeVFb+TaRqo6wK2S7pO0JE+AW8RmVmkvs/WW23pvmJFz9/369dsujYilfU8k3QYc3CDuixHxw5xlvC8i1kt6E7BC0qMRcedgAU7EZlZpEbGoicc6uQnHWJ/9u1HSzcCxwKCJ2F0TZmZNImmypP37HgOnUrvINygnYjOzHCSdJakLeA/wY0m3ZD+fI2l5ttss4N8kPQj8O/DjiPjpkMeOik8NNDOrOreIzcxK5kRsZlYyJ2Izs5I5EZuZlcyJ2MysZE7EZmYlcyI2MyuZE7GZWcn+P14K27yElBxXAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 2 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Rmgoai5TM3m-"
      },
      "source": [
        "condition number, with respect to the $\\|\\cdot\\|_2$ operator norm."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cKOy3mz7M1Yv",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "fc707ad8-1c17-44a2-d447-7181de5fc59b"
      },
      "source": [
        "np.linalg.cond(B_up.value)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "4.95279241640092"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 10
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mFgpPVm1bvZW",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "7e0a9e8b-9be2-416a-97a4-f0fb1e9d7d76"
      },
      "source": [
        "np.linalg.cond(B_cvx.value)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "5.417888775286857"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 11
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "m43VVod_NDiY"
      },
      "source": [
        "We observe $ \\mu_b(B_{cvx}) < \\mu_b(B_{up}) $ but on the other hand $ \\kappa(B_{cvx}) > \\kappa(B_{up}) $\n"
      ]
    }
  ]
}