notebooks/portfolio_optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1553f843",
   "metadata": {},
   "source": [
    "# Portfolio optimization under dominance constraints\n",
    "\n",
    "The problem solved here comes from Example 1 in [Xue et al. (2020)](https://www.aimsciences.org/article/doi/10.3934/jimo.2019071):\n",
    "\n",
    "$$\n",
    "\\begin{equation*}\n",
    "    \\min_{z} \\ \\ -\\mathbb{E}_P[G_2(\\xi)] \\ \\ \\\n",
    "    \\mathrm{s.t.} \\ \\ G_2(\\xi) \\succeq_2 G_1(\\xi), \\ 1\\leq z \\leq 2\n",
    "\\end{equation*}\n",
    "$$\n",
    "\n",
    "where $\\xi \\in \\mathbb{R} \\sim P$ is drawn uniformly from $[0,1]$, the benchmark portfolio is defined as:\n",
    "\n",
    "$$\n",
    "\\begin{equation*}\n",
    "    G_1(\\xi) = \n",
    "    \\begin{cases}\n",
    "        \\frac{i}{20} &\\xi \\in [0.05 \\times i, 0.05 \\times (i+1)) \\ \\ \\ \\ i=0, \\ldots, 19\\\\\n",
    "        1 &\\xi=1\n",
    "    \\end{cases}\n",
    "\\end{equation*}\n",
    "$$\n",
    "\n",
    "and the optimization is over the parametrized portfolio $G_2(\\xi) := G(\\xi;z) = z \\xi$.\n",
    "\n",
    "$Y_2 \\succeq_2 Y_1$ denotes 2nd order stochastic dominance, which holds if and only if $\\mathbb{E}[(\\eta-Y_2)_{+}] \\leq \\mathbb{E}[(\\eta - Y_1)_{+}]$ for any $\\eta \\in \\mathbb{R}$, or equivalently, if $\\mathbb{E}[\\tilde{u}(Y_2)] \\geq \\mathbb{E}[\\tilde{u}(Y_1)]$ for all concave non-decreasing $\\tilde{u} : \\mathbb{R} \\to \\mathbb{R}$ ([Dentcheva and Ruszczynski, 2004](https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.90.9157&rep=rep1&type=pdf)).\n",
    "\n",
    "We solve a relaxed version of this constrained optimization and replace the class of **non-decreasing concave** functions with that of **decreasing convex** functions.\n",
    "\n",
    "This is equivalent to using convex decreasing $u$, i.e. flipping the inequality in the dominance requirement as follows:\n",
    "\n",
    "$$\\mathbb{E}[u(Y_2)] \\leq \\mathbb{E}[u(Y_1)]$$\n",
    "\n",
    "In other words $Y_2 \\preceq_{dcvx} Y_1$, where $\\preceq_{dcvx}$ is the order generated by the class of decreasing convex functions in 1D.\n",
    "\n",
    "The relaxed problem is therefore:\n",
    "$$\\min_{z}-\\mathbb{E}_P[G(z\\xi)] + \\lambda\\sup_{u \\in dcvx}\\mathbb{E}_P[u(G(z\\xi))] - \\mathbb{E}_P[u(G_1(\\xi))] $$\n",
    "\n",
    "This can be optimized as the following bilevel optimization problem\n",
    "\n",
    "$$\\min_{z}-\\mathbb{E}_P[G(z\\xi)] + \\lambda[\\mathbb{E}_P[u(G(z\\xi))] - \\mathbb{E}_P[u(G_1(\\xi))]] $$\n",
    "\n",
    "and\n",
    " \n",
    "$$\\min_{u \\in dcvx} \\mathbb{E}_P[u(G_1(\\xi))] - \\mathbb{E}_P[u(G(z\\xi))] $$\n",
    "\n",
    "Finally, $dcvx$ is approximated with 3-layer, fully-connected, input convex max-out networks with non-positive input to hidden weights and non-negative weights for all subsequent layers."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d111f15",
   "metadata": {},
   "source": [
    "## 0. Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1b7fedf4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import os\n",
    "import sys\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import pytorch_lightning as pl\n",
    "import torch\n",
    "from torch import nn\n",
    "from torch.optim import SGD\n",
    "from tqdm.notebook import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "01e21149",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Global seed set to 0\n"
     ]
    }
   ],
   "source": [
    "seed = 0\n",
    "pl.seed_everything(seed=seed);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de6d2804",
   "metadata": {},
   "source": [
    "## 1. Modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "adf1345e",
   "metadata": {},
   "outputs": [],
   "source": [
    "class ConvexNN(nn.Module):\n",
    "    def __init__(self, input_dim=1, hidden_dim=128, output_dim=1, max_out=4, bias=False):\n",
    "        super().__init__()\n",
    "        max_out_adusted_hidden_dim = hidden_dim // max_out\n",
    "        self.input_to_hidden = nn.Sequential(\n",
    "            nn.Linear(input_dim, hidden_dim, bias=bias),\n",
    "            nn.MaxPool1d(kernel_size=max_out)\n",
    "        )\n",
    "        self.main = nn.Sequential(\n",
    "            nn.Linear(max_out_adusted_hidden_dim, hidden_dim, bias=bias),\n",
    "            nn.MaxPool1d(kernel_size=max_out),\n",
    "            nn.Linear(max_out_adusted_hidden_dim, output_dim, bias=bias)\n",
    "        )\n",
    "        self.project_params_for_decreasing_cvxty()\n",
    "    \n",
    "    def forward(self, x):\n",
    "        hidden = self.input_to_hidden(x)\n",
    "        return self.main(hidden)\n",
    "    \n",
    "    def project_params_for_decreasing_cvxty(self):\n",
    "        for module in self.input_to_hidden.modules():\n",
    "            if module and hasattr(module, 'weight'):\n",
    "                with torch.no_grad():\n",
    "                    module.weight.data.clip_(-torch.inf, 0)\n",
    "        for module in self.main.modules():\n",
    "            if module and hasattr(module, 'weight'):\n",
    "                with torch.no_grad():\n",
    "                    module.weight.data.relu_()\n",
    "\n",
    "    def get_weight_norms(self):\n",
    "        param_norm = 0\n",
    "        numel = 0\n",
    "        for name, module in self.named_modules():\n",
    "            param_norm = 0\n",
    "            numel = 0\n",
    "            if name == '':\n",
    "                continue\n",
    "            for param in module.parameters():\n",
    "                if param.requires_grad:\n",
    "                    with torch.no_grad():\n",
    "                        param_norm += param.norm(2).detach().item() ** 2\n",
    "                        numel += param.numel()\n",
    "        return math.sqrt(param_norm) / numel\n",
    "                    \n",
    "    def test_convexity(self, batch_size=64, eps=1e-6):\n",
    "        x = torch.randn((batch_size, 1))\n",
    "        with torch.no_grad(): \n",
    "            # Sample lambda uniformly\n",
    "            lmbds = torch.rand((batch_size, 1))\n",
    "\n",
    "            # Shuffle indices to randomly select images\n",
    "            indices1 = torch.randperm(x.shape[0])\n",
    "            indices2 = torch.randperm(x.shape[0])\n",
    "\n",
    "            # Check that `f(lambda x + (1-lambda)y) <= lambda f(x) + (1-lambda) f(y)`\n",
    "            convex_combo_of_inputs = self(lmbds * x[indices1] + (1 - lmbds) * x[indices2])\n",
    "            convex_combo_of_outputs = lmbds * self(x[indices1]) + (1 - lmbds) * self(x[indices2])\n",
    "            return torch.sum(convex_combo_of_inputs <= convex_combo_of_outputs + eps).item() / batch_size"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd4fbcfb",
   "metadata": {},
   "source": [
    "## 2. Optimization loop\n",
    "In the notation below, we use $Y$ to denote the benchmark."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dc6f3686",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_ξ_Y(batch_size):\n",
    "    ξ = torch.rand((bsz, 1))\n",
    "    Y = torch.floor(ξ*20) / 20\n",
    "    return ξ, Y\n",
    "\n",
    "def optimization_loop(z, u, z_opt, u_opt, choquet_weight, iterations, batch_size):\n",
    "    pbar = tqdm(range(iterations), desc='Iter')\n",
    "    z_steps, u_steps = 0, 0\n",
    "    steps_log, z_log = [], []\n",
    "    for i in pbar:\n",
    "        # Alternate training z and u\n",
    "        if i % 2:\n",
    "            u.eval()\n",
    "        else:\n",
    "            u.train()\n",
    "\n",
    "        # Sample ξ and calculate Y(ξ)\n",
    "        ξ, Y = sample_ξ_Y(batch_size)\n",
    "        \n",
    "        # Calculate objectives\n",
    "        portfolio = z*ξ\n",
    "        portfolio_objective = -portfolio.mean()\n",
    "        u_Y = u(Y).mean()\n",
    "        u_G = u(portfolio).mean()\n",
    "        choquet_objective = u_G - u_Y\n",
    "        \n",
    "        # Take gradient steps: alternate training z and u\n",
    "        if i % 2:  # optimize z\n",
    "            optimization_objective = portfolio_objective + choq_weight*choquet_objective\n",
    "            z_opt.zero_grad()\n",
    "            optimization_objective.backward()\n",
    "            z_opt.step()\n",
    "            z.data.clip_(1, 2)\n",
    "            z_steps += 1\n",
    "            steps_log.append(i)\n",
    "            z_log.append(z.data.item())\n",
    "        else:  # optimize u\n",
    "            choquet_objective_for_u = -choquet_objective\n",
    "            u_opt.zero_grad()\n",
    "            choquet_objective_for_u.backward()\n",
    "            u_opt.step()\n",
    "            u.project_params_for_decreasing_cvxty()       \n",
    "            u_steps += 1\n",
    "        \n",
    "        # Test u convexity\n",
    "        cvxty = u.test_convexity(batch_size=batch_size)\n",
    "        \n",
    "        # Log values to pbar\n",
    "        with torch.no_grad():\n",
    "            postfix_str = f'G(ξ,z): {portfolio_objective.item():0.2f}, ' + \\\n",
    "                          f'z: {z.data.item():0.2f}, ' + \\\n",
    "                          f'u(Y): {u_Y.item():0.2f}, ' + \\\n",
    "                          f'u(G(ξ,z)): {u_G.item():0.2f}, ' + \\\n",
    "                          f'cvx: {cvxty:0.2f}'\n",
    "            pbar.set_postfix_str(postfix_str)\n",
    "        pbar.set_description(f'u: {u_steps:,d}, z: {z_steps:,d}')\n",
    "\n",
    "    return {'portfolio': portfolio_objective.item(), 'choquet': choquet_objective.item(),\n",
    "            'steps_log': steps_log, 'z_log': z_log}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7dac348d",
   "metadata": {},
   "source": [
    "## 3. Run optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2114aab9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hyperparams\n",
    "u_hidden_dim = 32\n",
    "u_max_out = 4\n",
    "u_bias = True\n",
    "choq_weight = 1\n",
    "bsz = 512\n",
    "steps = 5000\n",
    "u_lr = 1e-3\n",
    "z_lr = 1e-3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "89d433ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Models and optimizer\n",
    "u = ConvexNN(hidden_dim=u_hidden_dim, max_out=u_max_out, bias=u_bias)\n",
    "u_opt = SGD(params=u.parameters(), lr=u_lr, maximize=False)\n",
    "\n",
    "z = nn.Parameter(torch.ones(1))\n",
    "z_opt = SGD(params=[z], lr=z_lr, maximize=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "eeb048dd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4433b5ce284d498ab7cb09142831151d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Iter:   0%|          | 0/5000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Run loop\n",
    "opt_output = optimization_loop(z=z, u=u, z_opt=z_opt, u_opt=u_opt,\n",
    "                               choquet_weight=choq_weight, iterations=steps, batch_size=bsz);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3bfdcfc0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "z           : 2.0000\n",
      "E[G(ξ,z)]   : 1.0424\n",
      "E[Y(ξ)]     : 0.4964\n",
      "E[u(G(ξ,z))]: -0.1116\n",
      "E[u(Y)]     : -0.1058\n"
     ]
    }
   ],
   "source": [
    "# Display results\n",
    "u.eval()\n",
    "ξ_val, Y = sample_ξ_Y(bsz*4)\n",
    "G = (z*ξ_val)\n",
    "print(f'z           : {z.data.item():0.4f}')\n",
    "print(f'E[G(ξ,z)]   : {G.mean().item():0.4f}')\n",
    "print(f'E[Y(ξ)]     : {Y.mean().item():0.4f}')\n",
    "print(f'E[u(G(ξ,z))]: {u(G).mean().item():0.4f}')\n",
    "print(f'E[u(Y)]     : {u(Y).mean().item():0.4f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "257272c4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAhl0lEQVR4nO3deXhU1f3H8fcXCFtAtgAiEAOyKSBbwiIuaN2K/krdqojIHlCkWm3VVlvb2sXWLtpaRSqYsMiiIGKrtu5YKZAECPsS9rAlEHYIZDm/P+ZiI7KTm5uZ+byeJ09m7r1MvidM5jP3njPnmHMOERGJXhWCLkBERIKlIBARiXIKAhGRKKcgEBGJcgoCEZEoVynoAs5WXFycS0hICLoMEZGwkpGRsdM5V/9E+8IuCBISEkhPTw+6DBGRsGJmG0+2T5eGRESinIJARCTKKQhERKKcgkBEJMopCEREopxvQWBmTc3sUzNbbmbLzOzhExxjZvYXM8sys8Vm1tmvekRE5MT8HD5aCDzmnFtgZjWBDDP70Dm3vMQx3wZael/dgFe87yIiUkZ8CwLn3DZgm3d7v5mtABoDJYOgDzDehebCnmtmtc2skfdvRSQKfLYqhwUbdwddRlhITKjL1a1O+Jmw81ImHygzswSgEzDvuF2Ngc0l7md7274WBGaWDCQDxMfH+1aniJSd4mLHnz9azV8/yQLALOCCwsCIay4JzyAwsxrAdOAR59y+c3kM59wYYAxAYmKiVtIRCXOHjhby6NRMPli2nbsTm/Lsd9tRuZLGrgTF1yAwsxhCITDJOTfjBIdsAZqWuN/E2yYiEWrLnsMMS01n5fZ9/PTWyxjcMwHT6UCgfAsCC/3PjgVWOOf+dJLDZgEPmdkUQp3Ee9U/IBK5MjbuZviEDI4UFDFuYBK9WjcIuiTB3zOCnkB/YImZLfK2/QSIB3DOjQbeA3oDWcAhYJCP9YhIgGYsyObJ6Uu4qHZVpiR3o0WDmkGXJB4/Rw39Bzjl+Z43WmikXzWISPCKih3P/2sVoz9fyxWX1OPlfp2pXb1y0GVJCWE3DbWIhI8DRwp5ZMpCPlqRw33d43nm/9oSU1GdwuWNgkBEfLE57xBDU9PJyj3As33a0r9HQtAlyUkoCESk1M1bt4sHJi2gqNgxfnBXeraIC7okOQUFgYiUqqlpm3h65lKa1q3O2AFJNIuLDbokOQ0FgYiUisKiYn7z3krGfbmeq1vV5699O1GrWkzQZckZUBCIyHnbl1/AqDcW8vnqXAb1TOCp3pdSSZ3CYUNBICLnZf3OgwxNTWPjrkP89vb29O2q+cDCjYJARM7Zl1k7eXDSAioYTBzaje7N6wVdkpwDBYGInJMJ/93Az99dziX1Yxk7IImmdasHXZKcIwWBiJyVgqJifvnucibM3ci32jTghXs6UrOqOoXDmYJARM7YnkNHeXDSAuas3cXwa5rz+E1tqFhBM4eGOwWBiJyRrJwDDE1NY+uefP5wVwfu7NIk6JKklCgIROS0PluVw6jJC6lSqQKTk7vR5eK6QZckpUhBICIn5Zxj3Jcb+PU/l9P6wgt4bUAijWtXC7osKWUKAhE5oaOFxfzsnaVMSdvMTW0b8qfvdSS2il4yIpH+V0XkG/IOHmXExAzmr89j1HUt+MH1raigTuGIpSAQka9ZtX0/Q1LTyN1/hBfv6Uifjo2DLkl8piAQka98vGIH35+8kNgqlZg6vAcdm9YOuiQpAwoCEcE5x5jZ63jug5W0u6gWf78/kQtrVQ26LCkjCgKRKJdfUMRP3l7CjAVbuPXyRjx/ZweqVa4YdFlShhQEIlEsd/8Rhk9IZ8GmPTx6QytGXdcCM3UKRxsFgUiUWrZ1L8NS08k7dJSX+3Wmd/tGQZckAVEQiEShD5Zu4wdTM6ldPYa3RlxBu8a1gi5JAqQgEIkizjle+iSLP364mo5NazPm/i40qKlO4WinIBCJEvkFRfzorcW8m7mV2zo15re3t6dqjDqFRUEgEhV27Mtn2Ph0lmzZyxM3t2HENc3VKSxfURCIRLjF2XsYNj6d/fmFjOmfyA2XNQy6JClnFAQiEezdzK388M1M4mpUYfoDV3BpowuCLknKIQWBSAQqLna88NFq/vJJFkkJdRh9Xxfq1agSdFlSTikIRCLMoaOFPDYtk/eXbud7iU341XfbU7lShaDLknJMQSASQbbuOczQ1HRWbt/H07dcypArm6lTWE7Lt7cJZjbOzHLMbOlJ9tcys3fNLNPMlpnZIL9qEYkGGRt3852XvmRz3iHGDkxi6FUaGSRnxs/zxRTg5lPsHwksd851AHoBfzSzyj7WIxKxZizIpu+YucRWqciMB6/g2tYNgi5Jwohvl4acc7PNLOFUhwA1LfSWpQaQBxT6VY9IJCoqdjz/r1WM/nwtPZrX4+V+nakTq/dTcnaC7CN4CZgFbAVqAnc754pPdKCZJQPJAPHx8WVWoEh5duBIIY9MWchHK3Lo1y2en3+nLTEV1SksZy/IZ81NwCLgIqAj8JKZnXCQs3NujHMu0TmXWL9+/bKrUKSc2px3iDtensOnq3L5ZZ+2/Oq77RQCcs6CPCMYBDznnHNAlpmtB9oA8wOsSaTcm78+jxETMygsKiZ1UFeubBkXdEkS5oJ8C7EJ+BaAmTUEWgPrAqxHpNybmraJfq/NpXa1GGaO7KkQkFLh2xmBmU0mNBoozsyygWeAGADn3GjgWSDFzJYABjzhnNvpVz0i4aywqJjfvr+Ssf9Zz1Ut43ipb2dqVY8JuiyJEH6OGup7mv1bgRv9+vkikWJffgGj3ljI56tzGXhFAk/fcimV1B8gpUifLBYpxzbsPMiQ1DQ27jrEb25rz73dNGpOSp+CQKScmpO1kwcmLcAMJgzpRo9L6gVdkkQoBYFIOTR5/iZ+OnMpzeJiGTsgifh61YMuSSKYgkCkHCkqdjz3/gr+/sV6rmlVn5fu7UTNquoUFn8pCETKiUNHC3l4yiI+XL6DAT0u5qe3XqZOYSkTCgKRcmD73nyGpKaxYts+fv5/lzGwZ7OgS5IooiAQCdjSLXsZkprGgfxCxg5I4to2mjlUypaCQCRA/162nYenLKJO9Rje0prCEhAFgUgAnHO89sV6fvP+Ci5vUpu/39+FBjWrBl2WRCkFgUgZKygq5mfvLGXy/M3c0r4Rf/xeB6rGVAy6LIliCgKRMrT3cAEPTsrgy6xdjLz2Eh67oTUVKmg5SQmWgkCkjGzcdZDBKWlsyjvE83dezl2JTYMuSQRQEIiUifQNeSRPyKDYOSYM6Ub35pouQsoPBYGIz2Yu3MLjby2mcZ1qjBuYRLO42KBLEvkaBYGIT5xz/PmjNfzl4zV0b16X0fd1oXZ1LSwv5Y+CQMQH+QVFPP7WYmZlbuWuLk349W3tqVxJ00VI+aQgECllOw8cIXl8Ogs27eHxm1vzwDWXYKaRQVJ+KQhEStHqHfsZnJLGzgNHeKVfZ77dvlHQJYmcloJApJTMXp3LyEkLqFq5IlOTe9Chae2gSxI5IwoCkVIwce5Gnpm1jJYNajB2YBKNa1cLuiSRM6YgEDkPRcWOX/9zBeO+XM91bRrwl76dqFFFf1YSXvSMFTlHB44U8vDkhXy8ModBPRN4+pbLqKjpIiQMKQhEzsHWPYcZkprO6h37ebZPW/r3SAi6JJFzpiAQOUtLskMLyRw6WsTYAYn0aq2FZCS8KQhEzsIHS7fxyNRF1IutwvQHutH6wppBlyRy3hQEImfAOcers9fx3Psr6RRfmzH9E6lfs0rQZYmUCgWByGkcLSzm6ZlLmJaeza2XN+IPd2khGYksCgKRU9h7qIAREzP477pdfP+6FjxyfSstJCMRR0EgchIbdoYWksnefZg/fa8Dt3duEnRJIr5QEIicwLx1uxg+MQMDJg7tRtdmdYMuScQ3vs2La2bjzCzHzJae4pheZrbIzJaZ2ed+1SJyNqZnZHPf2HnUja3MzJE9FQIS8fw8I0gBXgLGn2inmdUGXgZuds5tMjMNxpZAFRc7/vThal76NIsrLqnHK/26UKt6TNBlifjOtyBwzs02s4RTHHIvMMM5t8k7PsevWkROJ7+giMfezOSfi7dxT1JTnv1uO2IqaiEZiQ5B9hG0AmLM7DOgJvCic+5kZw/JQDJAfHx8mRUo0SF3/xGGjU8nM3sPP+ndhmFXNddCMhJVggyCSkAX4FtANeC/ZjbXObf6+AOdc2OAMQCJiYmuTKuUiLZy+z6GpKSTd/Aoo+/rwk1tLwy6JJEyF2QQZAO7nHMHgYNmNhvoAHwjCET88NmqHB56YyGxVSry5ogetGtcK+iSRAIR5EXQd4ArzaySmVUHugErAqxHosj4/25gcEoa8XWrM3NkT4WARDXfzgjMbDLQC4gzs2zgGSAGwDk32jm3wsw+ABYDxcBrzrmTDjUVKQ2FRcX86p8rSJmzgesvbcCL93QiVgvJSJTzc9RQ3zM45nngeb9qECnpwJFCRr2xgE9X5TL0ymb8uPelWkhGBH2yWKLElj2HGZKSxpqcA/z6tnb063Zx0CWJlBsKAol4izbvYWhqOkcKi0gZlMRVLesHXZJIuaIgkIj23pJt/GDqIhpcUIXJw7rRsqEWkhE5noJAIpJzjpc/W8vz/1pFl4vrMKZ/F+rV0EIyIieiIJCIc7SwmB/PWML0Bdn06XgRv7vjci0kI3IKCgKJKLsPHmXExAzmrc/jketb8vC3Wmq6CJHTUBBIxFiXe4Ahqels2X2YF+/pSJ+OjYMuSSQsKAgkIvx37S5GTMygYgXjjWHdSEzQGgIiZ0pBIGFvWvpmnnp7CRfXi2XcgCTi61UPuiSRsKIgkLBVXOx4/t+reOWztVzVMo6X7u1MrWpaSEbkbCkIJCwdPlrEo9MW8f7S7dzbLZ5ffKetFpIROUcKAgk7OfvyGTY+ncVb9vL0LZcy5MpmGhkkch4UBBJWlm/dx9DUNPYcLmBM/0RuuKxh0CWJhL0zOpc2s4/NrPdx28b4U5LIiX2ycgd3jZ5DsYNpw3soBERKyZleVG0GPGFmz5TYluhDPSLf4Jzj9S/XMzQ1nWb1Y3nnIS0kI1KazjQI9hBaW7ihmb1rZvorlDJRWFTMz95Zxi/eXc71lzZk2vAeNLygatBliUSUM+0jMOdcIfCgmQ0E/gPU8a0qEWBffgGj3ljI56tzGX51c564uQ0VtJCMSKk70yAYfeyGcy7FzJYAI/0pSQQ25x1iSGoa63IP8tzt7bmna3zQJYlErDMKAufcq8fdzwAG+1KRRL0Fm3aTPD6do4XFpA7uSs8WcUGXJBLRNHxUypV3M7fy2JuZNKpVlSnJSbRoUCPokkQinoJAygXnHC99ksUfP1xNUkIdXu2fSN3YykGXJRIVFAQSuCOFRfx4+hJmLNzCbZ0a89wd7alSSQvJiJQVBYEEKu/gUYZPSCdtw24eu6EVD13XQtNFiJQxBYEEJivnAENS09i2N5+/9u3E/3W4KOiSRKKSgkACMSdrJyMmZlC5UgWmJHenc7w+liISFAWBlLlpaZv5ydtLaF4/lrEDkmhaVwvJiARJQSBlprjY8ccPV/G3T0MLyfytX2cuqKqFZESCpiCQMpFfUMSP3lrMu5lb6du1Kb/s004LyYiUEwoC8V3ewaMMG59OxsbdPPntNgy/urlGBomUIwoC8dW63AMMSklj+958Xu7Xmd7tGwVdkogcR0Egvpm3bhfJEzKoVMGYrJFBIuWWbxdpzWycmeWY2dLTHJdkZoVmdqdftUjZe3thNveNnUdcjcq8/WBPhYBIOeZnb10KcPOpDjCzisDvgH/7WIeUIeccL3y0mh9MzSTx4rrMeKAn8fU0PFSkPPPt0pBzbraZJZzmsFHAdCDJrzqk7BwpLOLJ6Ut4e+EW7uzShN/c1p7KlTQySKS8C6yPwMwaA7cB13KaIDCzZCAZID5eC5SUR3sOHSV5Qgbz1+fxwxtbMfJazRkkEi6C7Cx+AXjCOVd8uhcM59wYYAxAYmKi8780ORsbdx1k0OtpZO8+zIv3dKRPx8ZBlyQiZyHIIEgEpnghEAf0NrNC59zMAGuSs5SxMY9h4zNwzjFpWDeSEuoGXZKInKXAgsA51+zYbTNLAf6hEAgv/1i8lUenZdK4djXGDUyiWVxs0CWJyDnwLQjMbDLQC4gzs2zgGSAGwDk32q+fK/5zzvHq7HU89/5KkhLqMKZ/InW0mphI2PJz1FDfszh2oF91SOkqLCrmmVnLmDRvE7de3og/3NWBqjFaTUwknOmTxXLGDh4pZNTkhXyyMocR11zC4ze1pkIFjQwSCXcKAjkjOfvyGZyaxvKt+/j1be3o1+3ioEsSkVKiIJDTWr1jP4NeT2P3oaOMHZDEtW0aBF2SiJQiBYGc0pysnQyfmEG1mIpMG96Ddo1rBV2SiJQyBYGc1PSMbJ6csZhmcbG8PqgrjWtXC7okEfGBgkC+wTnHXz/J4k8fruaKS+rxyn1dqFVNS0qKRCoFgXxNQVExP525lClpm7m9U2Oeu+NyTRwnEuEUBPKVA0cKGTlpAZ+vzmXUdS149IZWmjhOJAooCAQIDQ8dlJLGyu37+e3t7enbVbO8ikQLBYGwZsd+BnrDQ18bkMi1rTU8VCSaKAii3Nx1u0gen07lShWZmtyD9k00PFQk2igIotiszK38cFomTetWI2VQV5rW1ZKSItFIQRCFSs4e2rVZXf7eP5Fa1TU8VCRaKQiiTFGx4+ezljFh7kbNHioigIIgqhw+WsSoyQv5aMUOhl/dnCdubqPZQ0VEQRAtdh44wpDUdBZn7+GXfdpyf4+EoEsSkXJCQRAF1uUeYODraeTsz+fV+7pwY9sLgy5JRMoRBUGEy9iYx9DUdMyMycO60ym+TtAliUg5oyCIYB8s3cbDUxbRqFZVUgZ1JUGLy4vICSgIItS4/6zn2X8up2PT2rx2fyL1alQJuiQRKacUBBGmuNjx6/dWMPY/67nxsoa8eE8nqlXW8FAROTkFQQTJLyji0WmLeG/JdgZekcBPb72MihoeKiKnoSCIELsPHmXY+HTSN+7m6VsuZciVzTSFtIicEQVBBNicd4gBr88ne/dh/nZvZ265vFHQJYlIGFEQhLnF2XsYnJJGQZFj0tBuJCXUDbokEQkzCoIw9vGKHTz0xkLq1ajMlOSutGhQI+iSRCQMKQjC1KR5G/npzKW0vagWYwcm0qBm1aBLEpEwpSAIM8XFjuf/vYpXPlvLta3r89K9nYmtov9GETl3egUJI0cLi3n8rUxmLtpK367xPNunLZUqVgi6LBEJcwqCMLH3cAEjJmTw33W7+NFNrXmw1yUaHioipcK3t5NmNs7Mcsxs6Un29zOzxWa2xMzmmFkHv2oJd1v2HOau0XNI35jHn+/uwMhrWygERKTU+HldIQW4+RT71wPXOOfaA88CY3ysJWwt27qX21/+km178kkd1JXbOjUJuiQRiTC+XRpyzs02s4RT7J9T4u5cQK9wx5m9OpcHJmZwQbUY3nrgClpfWDPokkQkApWXnsYhwPsn22lmyWaWbmbpubm5ZVhWcN5M38zglDSa1q3O2w/2VAiIiG8C7yw2s2sJBcGVJzvGOTcG79JRYmKiK6PSAuGc48WP1/DCR2u4qmUcL/frTM2qMUGXJSIRLNAgMLPLgdeAbzvndgVZS3lQUFTMT2Ys4c2MbO7o3ITn7mhPjIaHiojPAgsCM4sHZgD9nXOrg6qjvNifX8CDkxbwxZqdfP9bLfnB9S01MkhEyoRvQWBmk4FeQJyZZQPPADEAzrnRwM+AesDL3gteoXMu0a96yrMd+/IZ+Hoaq3fs53d3tOfupPigSxKRKOLnqKG+p9k/FBjq188PF6t37GfguPnsPVzAuIFJXNOqftAliUiUCbyzOJrNWbuT4RMyqBZTkanDe9Cuca2gSxKRKKQgCMg7i7bwwzczubheLCmDkmhSp3rQJYlIlFIQlDHnHK98vpbff7CKbs3qMqZ/IrWqa3ioiARHQVCGCouKeWbWMibN28R3OlzE83ddTpVKFYMuS0SinIKgjBw6WsioNxby8cocRlxzCY/f1JoKFTQ8VESCpyAoA7n7jzA0NY0lW/by7Hfb0b/7xUGXJCLyFQWBz9bvPMiAcfPJ2Z/Pq/0TueGyhkGXJCLyNQoCHy3ctJshqekATB7WnU7xdQKuSETkmxQEPvl4xQ5GvrGABjWrkjq4K83iYoMuSUTkhBQEPpg8fxNPvb2Edo1rMXZAEvVrVgm6JBGRk1IQlCLnHC98tIYXP15Dr9b1+du9nYmtol+xiJRvepUqJYVFxTw9cylT0jZzZ5cm/PZ2TSEtIuFBQVAKSn5G4KFrW/DYja00hbSIhA0FwXnKO3iUwSlpZGbv0WcERCQsKQjOw+a8Q9w/bj5b9xzmlX5duLndhUGXJCJy1hQE52jplr0MfD2NgqJiJg3tRmJC3aBLEhE5JwqCc/DFmlxGTMigVrUYJg/rQcuGNYMuSUTknCkIztLbC7P50ZuLadGgBimDunJhrapBlyQicl4UBGfIOcers9fx3Psr6d68LmPuT+SCqlpHQETCn4LgDBQVO579x3JS5mzg1ssb8cfvddA6AiISMRQEp5FfUMSj0xbx3pLtDLmyGU/1vlTrCIhIRFEQnMLewwUkj09n3vo8nup9KcOubh50SSIipU5BcBLb9h5mwLj5rN95kBfv6Uifjo2DLklExBcKghNYvWM/A8bNZ39+ISmDutKzRVzQJYmI+EZBcJz56/MYmppGlZiKTB3enbYX1Qq6JBERXykISnh/yTYenrqIJnWqkTqoK03rVg+6JBER3ykIPKlzNvDzd5fRqWltxg5Iok5s5aBLEhEpE1EfBM45fv+vVbzy2Vquv7Qhf+3biWqV9RkBEYkeUR0EBUXFPDF9MTMWbKFv13ie7dOWSlpMRkSiTNQGwYEjhTwwMYMv1uzk0RtaMeq6FlpMRkSiUlQGQe7+IwxKmc+Kbfv53R3tuTspPuiSREQC49t1EDMbZ2Y5Zrb0JPvNzP5iZllmttjMOvtVS0nrdx7k9le+ZG3OQf5+fxeFgIhEPT8viKcAN59i/7eBlt5XMvCKj7UAsGjzHu54ZQ4HjxQxObk717Vp6PePFBEp93wLAufcbCDvFIf0Aca7kLlAbTNr5Fc9X6zJpe+YucRWqchbI3rQsWltv36UiEhYCXKITGNgc4n72d62bzCzZDNLN7P03Nzcc/thtauRmFCH6Q9cQfP6Nc7pMUREIlFYjJV0zo1xziU65xLr169/To/RvH4NJgzpRoOaWlFMRKSkIINgC9C0xP0m3jYRESlDQQbBLOB+b/RQd2Cvc25bgPWIiEQl3z5HYGaTgV5AnJllA88AMQDOudHAe0BvIAs4BAzyqxYRETk534LAOdf3NPsdMNKvny8iImcmLDqLRUTEPwoCEZEopyAQEYlyCgIRkShnoT7b8GFmucDGc/znccDOUiwnHKjN0UFtjg7n0+aLnXMn/ERu2AXB+TCzdOdcYtB1lCW1OTqozdHBrzbr0pCISJRTEIiIRLloC4IxQRcQALU5OqjN0cGXNkdVH4GIiHxTtJ0RiIjIcRQEIiJRLmqCwMxuNrNVZpZlZk8GXc/5MLNxZpZjZktLbKtrZh+a2Rrvex1vu5nZX7x2LzazziX+zQDv+DVmNiCItpwJM2tqZp+a2XIzW2ZmD3vbI7nNVc1svpllem3+hbe9mZnN89o21cwqe9urePezvP0JJR7rx972VWZ2U0BNOmNmVtHMFprZP7z7Ed1mM9tgZkvMbJGZpXvbyva57ZyL+C+gIrAWaA5UBjKBy4Ku6zzaczXQGVhaYtvvgSe9208Cv/Nu9wbeBwzoDszzttcF1nnf63i36wTdtpO0txHQ2btdE1gNXBbhbTaghnc7BpjntWUacI+3fTTwgHf7QWC0d/seYKp3+zLv+V4FaOb9HVQMun2nafujwBvAP7z7Ed1mYAMQd9y2Mn1uR8sZQVcgyzm3zjl3FJgC9Am4pnPmnJsN5B23uQ+Q6t1OBb5bYvt4FzIXqG1mjYCbgA+dc3nOud3Ah8DNvhd/Dpxz25xzC7zb+4EVhNa3juQ2O+fcAe9ujPflgOuAt7ztx7f52O/iLeBbZmbe9inOuSPOufWE1v/o6n8Lzo2ZNQFuAV7z7hsR3uaTKNPndrQEQWNgc4n72d62SNLQ/W+Ft+1AQ+/2ydoelr8T7/S/E6F3yBHdZu8SySIgh9Af9lpgj3Ou0DukZP1ftc3bvxeoR5i1GXgBeBwo9u7XI/Lb7IB/m1mGmSV728r0ue3bwjQSHOecM7OIGxdsZjWA6cAjzrl9oTd/IZHYZudcEdDRzGoDbwNtgq3IX2Z2K5DjnMsws14Bl1OWrnTObTGzBsCHZray5M6yeG5HyxnBFqBpiftNvG2RZId3ioj3PcfbfrK2h9XvxMxiCIXAJOfcDG9zRLf5GOfcHuBToAehSwHH3sCVrP+rtnn7awG7CK829wS+Y2YbCF2+vQ54kchuM865Ld73HEKB35Uyfm5HSxCkAS290QeVCXUszQq4ptI2Czg2UmAA8E6J7fd7ow26A3u9U85/ATeaWR1vRMKN3rZyx7vuOxZY4Zz7U4ldkdzm+t6ZAGZWDbiBUN/Ip8Cd3mHHt/nY7+JO4BMX6kWcBdzjjbBpBrQE5pdJI86Sc+7HzrkmzrkEQn+jnzjn+hHBbTazWDOreew2oefkUsr6uR10j3lZfRHqbV9N6DrrU0HXc55tmQxsAwoIXQscQuja6MfAGuAjoK53rAF/89q9BEgs8TiDCXWkZQGDgm7XKdp7JaHrqIuBRd5X7whv8+XAQq/NS4GfedubE3pRywLeBKp426t697O8/c1LPNZT3u9iFfDtoNt2hu3vxf9GDUVsm722ZXpfy469NpX1c1tTTIiIRLlouTQkIiInoSAQEYlyCgIRkSinIBARiXIKAhGRKKcgkKhlZge87wlmdm8pP/ZPjrs/pzQfX6Q0KQhEIAE4qyAo8UnXk/laEDjnrjjLmkTKjIJABJ4DrvLmg/+BN9nb82aW5s35PhzAzHqZ2RdmNgtY7m2b6U0WtuzYhGFm9hxQzXu8Sd62Y2cf5j32Um8O+rtLPPZnZvaWma00s0lWcjIlER9p0jmR0HzvP3TO3QrgvaDvdc4lmVkV4Esz+7d3bGegnQtNbwww2DmX500DkWZm051zT5rZQ865jif4WbcDHYEOQJz3b2Z7+zoBbYGtwJeE5t75T2k3VuR4OiMQ+aYbCc3nsojQdNf1CM1XAzC/RAgAfN/MMoG5hCb9asmpXQlMds4VOed2AJ8DSSUeO9s5V0xoGo2EUmiLyGnpjEDkmwwY5Zz72qRd3tTIB4+7fz3Qwzl3yMw+IzT/zbk6UuJ2Efr7lDKiMwIR2E9oCcxj/gU84E19jZm18maGPF4tYLcXAm0ILR14TMGxf3+cL4C7vX6I+oSWHS2XM2NK9NA7DpHQDJ9F3iWeFEJz4CcAC7wO21z+t1RgSR8AI8xsBaFZLueW2DcGWGxmC1xoKuVj3ia0rkAmoRlVH3fObfeCRCQQmn1URCTK6dKQiEiUUxCIiEQ5BYGISJRTEIiIRDkFgYhIlFMQiIhEOQWBiEiU+385RJ64O+ASigAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot z\n",
    "plt.plot(opt_output['steps_log'], opt_output['z_log'])\n",
    "plt.xlabel('Iteration')\n",
    "plt.ylabel('z')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcba4972",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}