liangyy/minimal_example_08092021.ipynb

## minimal_example_08092021.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run export PYTHONPATH=path-to-transethnic_prs:$PYTHONPATH before running this notebook\n",
    "import time\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.linear_model import lasso_path, enet_path\n",
    "import transethnic_prs.model1.Model1Blk as model1blk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulating data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# X1: n1 x p\n",
    "# y1: n1 \n",
    "# X2: n2 x p\n",
    "# y2: n2 \n",
    "# beta: p x 1\n",
    "# beta ~iid pi0 * delta(0) + (1 - pi0) * N(0, 1)\n",
    "# X1 ~iid N(0, 1)\n",
    "# X2 ~iid N(0, 1)\n",
    "# y1 = X1 * beta + N(0, s2 = 20)\n",
    "# y2 = X2 * beta + N(0, s2 = 20)\n",
    "\n",
    "n1 = 500\n",
    "n2 = 300\n",
    "p = 1020\n",
    "pi0 = 0.9\n",
    "X1 = np.random.normal(size=(n1, p))\n",
    "X2 = np.random.normal(size=(n2, p))\n",
    "beta = np.random.normal(size=(p))\n",
    "zero_ind = np.random.rand(p) < pi0\n",
    "beta[zero_ind] = 0\n",
    "y1 = X1 @ beta + np.random.normal(size=(n1), scale=20)\n",
    "y2 = X2 @ beta + np.random.normal(size=(n2), scale=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Centering data and calculating X'y and X'X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "X1 = X1 - X1.mean(axis=0)\n",
    "y1 = y1 - y1.mean()\n",
    "X2 = X2 - X2.mean(axis=0)\n",
    "y2 = y2 - y2.mean()\n",
    "b1 = X1.T @ y1  # e.g. b1 = (N1 - 1) diag(X1.cov()) bhat1 where bhat1 is the GWAS effect size estimate\n",
    "A1 = X1.T @ X1  # e.g. A1 = (N1 - 1) X1.cov()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving with my solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod1 = model1blk.Model1Blk([A1], [b1], [X2], y2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Run time = 0.9326646327972412 s\n"
     ]
    }
   ],
   "source": [
    "# elastic net with l1_ratio = 0.1\n",
    "# see l1_ratio at https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.enet_path.html\n",
    "l1_ratio = 0.1\n",
    "t = time.time()\n",
    "beta_mat_en, lambda_seq_en, niters_en, tols_en, convs_en = mod1.solve_path(alpha=l1_ratio) \n",
    "print(f'Run time = {time.time()-t} s')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Run time = 13.78531789779663 s\n"
     ]
    }
   ],
   "source": [
    "# lasso \n",
    "t = time.time()\n",
    "beta_mat_lasso, lambda_seq_lasso, niters_lasso, tols_lasso, convs_lasso = mod1.solve_path(alpha=1) \n",
    "print(f'Run time = {time.time()-t} s')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving with sklearn function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Run time = 0.7238945960998535 s\n"
     ]
    }
   ],
   "source": [
    "# elastic net with l1_ratio = 0.1\n",
    "# need eps=0.01 since my solver by default set lambda_min = lambda_max * 0.01\n",
    "t = time.time()\n",
    "alphas_enet, coefs_enet, kk = enet_path(\n",
    "    np.concatenate([X1, X2], axis=0), \n",
    "    np.concatenate([y1, y2]), \n",
    "    l1_ratio=l1_ratio, \n",
    "    fit_intercept=False, \n",
    "    eps=0.01)\n",
    "print(f'Run time = {time.time()-t} s')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Run time = 2.4606192111968994 s\n"
     ]
    }
   ],
   "source": [
    "# lasso\n",
    "t = time.time()\n",
    "alphas_lasso, coefs_lasso, kk_lasso = lasso_path(\n",
    "    np.concatenate([X1, X2], axis=0), \n",
    "    np.concatenate([y1, y2]), \n",
    "    fit_intercept=False, \n",
    "    eps=0.01)\n",
    "print(f'Run time = {time.time()-t} s')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparing results\n",
    "\n",
    "The function solve for the whole regularization path from lambda_max (min lambda such that the solution is all-zero) to lambda = lambda_max / 100.\n",
    "Here we compare at lambda values along the path: last one (idx = 99) and the one right in the middle (idx = 50)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f9901e5e8d0>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAdR0lEQVR4nO3df5xVdb3v8ddnNgMNaoxEKjMwYhyumQdTmwCTOh6TVOoEeRQxvHmzIs9Db8fyaCDIL0EMjxzjqCkW9+SJVErccZLiYuVD8QgXdBAkJNFQ2UMiKihKwMx87h97Dw7D/r3X/v1+Ph7zmL3X+u69Po8lvvnyXd/1XebuiIhI5aspdgEiIlIYCnwRkSqhwBcRqRIKfBGRKqHAFxGpEj2KXUAy/fr180GDBhW7DBGRsvHMM8/scvePxttX0oE/aNAg1q1bV+wyRETKhpm9kmifhnRERKqEAl9EpEoo8EVEqoQCX0SkSijwRUSqREnP0hERqSbhlgi3rdhC6+59NNTXcf35JzP2jMbAvl+BLyJSAsItESYv3ci+g+0ARHbvY/LSjQCBhb6GdERESsBtK7YcCvtO+w62c9uKLYEdI+fAN7OBZvYHM9tsZpvM7J/jtDEzW2BmW81sg5mdmetxRUQqSevufRltz0YQPfw24Dp3PwUYAVxtZp/o1uZCYEjsZyLwowCOKyJSMRrq6zLano2cA9/dd7j7s7HX7wKbge4DTmOA+z1qNVBvZv1zPbaISKW4/vyTqasNHbatrjbE9eefHNgxAr1oa2aDgDOANd12NQKvdXm/PbZtR5zvmEj0XwE0NTUFWZ6ISFElm4XT+bssZumY2dHAw8C17v5O991xPhL3YbruvhBYCNDc3KwH7opIRUhnFs7YMxoDDfjuApmlY2a1RMN+sbsvjdNkOzCwy/sBQGsQxxYRKQeFmIWTShCzdAz4CbDZ3ecnaLYM+Fpsts4IYI+7HzGcIyJSqQoxCyeVIIZ0zgb+J7DRzNbHtt0INAG4+z3AcmA0sBV4H/h6AMcVESkbDfV1ROKEe5CzcFLJOfDdfRXxx+i7tnHg6lyPJSJSrq4//+TDxvAh+Fk4qWhpBRGRAijELJxUFPgiIgWS71k4qWgtHRGRKqHAFxGpEgp8EZEqocAXEakSCnwRkSqhwBcRqRIKfBGRKqHAFxGpEgp8EZEqocAXEakSWlpBRKpSsqdPVSoFvohUnXSePlWJNKQjIlWnFJ4+VQzq4YtIxes+fBPvQSRQ2KdPFYMCX0QqWrzhGwM8TttCPn2qGDSkIyIVLd7wjXPkY/oK/fSpYlDgi0hFSzRM40BjfR0W+z33oqEVfcEWAhrSMbNFwJeAne7+t3H2nwP8CvhzbNNSd58VxLFFRJJJNGbfWF/HU5POLUJFxRNUD/8/gAtStHnS3U+P/SjsRaQgrj//ZOpqQ4dtq4bhm3gC6eG7+xNmNiiI7xIRCVIpPDy8VBRyls5ZZvYc0Ar8i7tvitfIzCYCEwGampoKWJ6IVKpiPzy8VBQq8J8FTnT3vWY2GggDQ+I1dPeFwEKA5ubmeDOnRKTKVeOyCEEoyCwdd3/H3ffGXi8Has2sXyGOLSKVpXNefWT3PpwPlkUIt0SKXVrJK0jgm9kJZmax18Nix32zEMcWkcpSrcsiBCGoaZkPAOcA/cxsOzAdqAVw93uAi4F/MrM2YB8w3t01XCMiGUs0r77Sl0UIQlCzdC5Lsf9O4M4gjiUi1S3RvPpKXxYhCLrTVkTKiubVZ0+Lp4lIWdG8+uwp8EWk7GhefXY0pCMiUiUU+CIiVUKBLyJSJTSGLyIFpWURikeBLyIFE+9xg5OXbgRQ6BeAAl9E8q6zVx/vhqnOZREU+PmnwBeRvJoa3sji1a/GfWh4Jy2LUBgKfBHJi6nhjfxs9atptdWyCIWhwBeRQIVbInz/4Q3sb+tIq72WRSgcBb6IBCaTXj1EHySuWTqFo8AXkUCEWyIsTjPsDfi3S09X0BeYAl9EcpJsBk4iE0Y0KeyLQIEvIlnrPq8+HUOOO4rZY4fmsSpJREsriEjW4j1uMJmzB/dl5ffOyV9BkpR6+CKStu7LIqQzjNO7toZbLjpNQzglQIEvImnpPgMnVdhrBk7pUeCLSEoT7nuap156K622dbUh5l40VEFfggIZwzezRWa208yeT7DfzGyBmW01sw1mdmYQxxWR/Jsa3pgy7Bvr67DYb4V96Qqqh/8fwJ3A/Qn2XwgMif0MB34U+y0iJSqT6ZZPTTq3ABVJrgIJfHd/wswGJWkyBrjf3R1YbWb1Ztbf3XcEcXwRCU64JcKMZZvYve9gWu2P7V2b54okKIWaltkIvNbl/fbYtiOY2UQzW2dm6954442CFCciUZ3z6tMN+xqD6f9wap6rkqAU6qKtxdkWd7VUd18ILARobm5OtqKqiAQg0x59p54hY97Fn9R4fRkpVOBvBwZ2eT8AaC3QsUUkgeFzVvL6uwcy+kzIjMuGD9TdsmWoUIG/DLjGzB4kerF2j8bvRYon01UtQdMtK0EggW9mDwDnAP3MbDswHagFcPd7gOXAaGAr8D7w9SCOKyKZyybsj+1dy/R/OFVhX+aCmqVzWYr9DlwdxLFEJDvZrGqpu2Uri+60FakC2YzV36H16iuOAl+kgo2a/zgv7nwv48+dPbivwr4CKfBFKlQ2vXqIrle/+Ftn5aEiKTYFvkiFybZXD3D5iCZNt6xgCnyRCvLxKcv5a3tm9ytqumX1UOCLVIBsploChAyFfRVR4IuUuVzG6vW4weqiwBcpU5k8lKSrD/cKsWHmBXmoSEqdAl+kzIRbIlz70PqsPquLstVNgS9SRrIdvlHQCyjwRcpCtr36HgZb534x+IKkLCnwRUrcadN/yzv72zP+3PHH9GTNlFF5qEjKlQJfpIRlE/aafSOJKPBFStBdd8xhzFuLWG+7aO3Zj3lt41jWMTLl57TgmSSjwBcpMd+dMpk5PX5M75roxdkBtotba38MB0kY+hqrl3Qo8EVKQLglwnVL1tPusKrnEnrb4TNxetsBbuixhGUHjgx8DeFIuhT4IkV21x1zGP/W3WztuTdpuwZ787D3IYPbx2kIR9KnwBcpknBLhN//4k5uq72XXjWpL8y2+kcOvdZYvWRDgS9SBJ03UK3quYReljrs3/eezGsbx9mD+2qteslaTRBfYmYXmNkWM9tqZpPi7D/HzPaY2frYz7QgjitSbibc9zSDJj166G7ZBtuVsK07dLixvaMfkw5+kw8P+6rCXnKScw/fzELAXcAoYDuw1syWufsfuzV90t2/lOvxRMrVz266mJ/W/I5Qrw7aqWFx+7m0ej8GJAj9iPdj5IEF0cXObtFiZ5K7IIZ0hgFb3f1lADN7EBgDdA98kaoUbolw8iOjmFATwSy6rQcdfC30GE92nMpH/e0jhnUOeA/mtY3TWL0EKoghnUbgtS7vt8e2dXeWmT1nZr8xs1MTfZmZTTSzdWa27o033gigPJHiCbdE6Lv0Ej5uH4R9JzP4TM1mrj/4bd7sOBr36DDOW340/1r3HRbcMldhL4EKoodvcbZ1f8bas8CJ7r7XzEYDYWBIvC9z94XAQoDm5ubMntUmUiI616r/cs0qfli76Yiw7xSig2UdIw+bX3/Hpadzo4Je8iCIwN8ODOzyfgDQ2rWBu7/T5fVyM7vbzPq5e+IrViJlatCkR/lyzSpW9VxCo+1KGPYA7V3+kf2hkPHCnNEFqFCqVRCBvxYYYmYnARFgPPDVrg3M7ATgdXd3MxtGdCjpzSO+SaSMTbjvab796nX8udcmgKRBD9Hhm8Xt5wJar14KI+fAd/c2M7sGWAGEgEXuvsnMrortvwe4GPgnM2sD9gHj3V3DNVIxHph+Mf/JSqwmddBDNOxf8EburruKbVrCWArESjl3m5ubfd26dcUuQyShtcvu5ZRnbuIo9qcV9BAN+yc7TuVzN/93fouTqmRmz7h7c7x9utNWJAtLFt3O2Ffm0Ex7RkEf8X78u43nBzfPyW+BInEo8EUy9J0bJ3NH7d3UpBn0EF0aYdLBb7Lglrn8IH+liSSlwBfJQKZh7w77qeXXJ05iwZXX5bc4kRQU+CJpeGD6xVzKSn5Ym95FWfjgwuwps/7IuPyWJ5KWQBZPE6lUU8MbiUw7kfGspMbSm2rpDnv9Q/xL+9WcMksrjEjpUA9fJIHvTpnMv/a4O62ghw969JfafDbMvIDb81+iSEYU+CLd7LjldE7Y/2fm98h8+Gb2wJ+wQUsYS4lS4It08dfpx3ICHWkHPcTumO0YxeU3/5LF+StNJGcKfBFgw5y/Y+iB9fQis179fkJ86B/v4fLTdFlWSp8CX6peZNqJDLXdGQU9wP3t53HF7IfzV5hIwBT4UrVen3ESx/lbNKR5URaiYf+219F31l+4Ir/liQRO0zKlKnVM68Nx/haWwQwcd1jVcSp9Z/0l/wWK5IF6+FJV9s5s5KiOvWkHPXwQ9tPOXKUljKWsKfClOvz6e/jan3AUmQU9QKvX0zjrFWbnrTiRwlDgS8XrmNEH8/SDHtSrl8qkwJeK9c6MBo7x9zAy79V3XphVr14qiQJfKlLHtD4ck8E4PXzQq6+ZtYe++StNpGg0S0cqyuoZI/HpfbK6KPsiA6iZtSe/BYoUkXr4UjE6pvVheA69+v+Rv9JESoICX8petlMtAdbYUEbMXJW/4kRKSCBDOmZ2gZltMbOtZjYpzn4zswWx/RvM7MwgjivSMa1P1vPqbeYeRsxQ2Ev1yLmHb2Yh4C5gFLAdWGtmy9y965MfLgSGxH6GAz+K/RbJStuMekLuWfXq91lPes98I3/FiZSoIIZ0hgFb3f1lADN7EBgDdA38McD97u7AajOrN7P+7r4jgONLlemY1odQDmP1vfNXmkhJCyLwG4HXurzfzpG993htGoEjAt/MJgITAZqamgIoTypF+/Q+1EBWvfq/9DqJ/jeuz1dpImUhiMCP97+eZ9EmutF9IbAQoLm5OW4bqTIbltDx8LcOhX26uvbq++etOJHyEUTgbwcGdnk/AGjNoo3IEdpm9CHkUJNh0AO8V3M0R0+P5KcwkTIUROCvBYaY2UlABBgPfLVbm2XANbHx/eHAHo3fSyq5jtUfnb/SRMpSzoHv7m1mdg2wAggBi9x9k5ldFdt/D7AcGA1sBd4Hvp7rcaVyvX3zYOrbdmU1Vm+f/gb2pfn5K06kjAVy45W7Lyca6l233dPltQNXB3EsqVy33DqTSfvmU0/mQd/ZqxeRxHSnrZSEyLQTmZzBc2Xh8OGbDD4mUrW0eJoU1dML/hc+vQ8NGT5EXIudiWROPXwpmsi0ExmRQ69ei52JZEaBLwX3p9s+z5C962jI4qJshxmhmbvzVptIJVPgS0HtndaPIXYw6159KH+liVQ8jeFLQaxddi/t0/twVAZh3xn0L9tAjdWLBEA9fMm7J276DJ+t2ZRxr77docesPQzOX2kiVUU9fMmbqeGNHJjWJ6OwP7RW/cf+jh7q1YsESj18CVy4JcKeh7/DzaHHIMMLs0/53zJy1lP5LVCkSinwJVBrZ36WMR0bIJT5E6h+1PcGrr52Sn4LFKliCnwJRLglQt+ll2Q1Vv+yDWTwzOe19oZIninwJWfh2eP50sEVhGo6MurVtwNLT5zGuCuvy2t9IhKlwJesLVl0Oxdu+wFjbH/GN1Dt7DeC4//3CsblrzwR6UaBL1m5adZNTG+/kx41HWl/putUy+PzWJuIxKfAl4xMuO9pnnrpLZ7vtZAell7Yd/bq19WcxqenP5nH6kQkGQW+pO1vJj9KWyy8j+KvSdt2hvx+ejCl/dvcPvsWPp3n+kQkOQW+pOU7N07m8dolNNguWr1f0rbucH/7eUxvu5Ihxx3Fyu+dU5giRSQpBb4kNWr+49z09o38sPaD6ZYDbNehHnx3nWF/d91VbJsyqnCFikhKCnyJq3OsfmaPRXw2dOTcerPYEghdtneG/YvNM1gzdmhhCxaRlHIKfDPrCzwEDAK2AePc/e047bYB7xKdet3m7s25HFfy6+NTlvPX9mgXfkLo90mnXLZ5DSE6aKeGX/B5rpj9ywJVKSKZynXxtEnA79x9CPC72PtE/t7dT1fYl65wS4RBkx49FPYAIRLPxIl4P/5m/884af/P+fXY57lspsJepJTlOqQzBjgn9vqnwOPA93P8TimCJYtu5zPb7ublXtGLsvPaxrGsYyTt1NAjTui7w7y2cbooK1JGzBNdfUvnw2a73b2+y/u33f3YOO3+DLwNOHCvuy9M8p0TgYkATU1Nn3rllVeyrk9SGzX/cU7ZtYJba39MbztwaPv73pNJB7/Jp2r+xNdCjx0xVv9kx6l87ub/LkLFIpKMmT2TaCQlZQ/fzB4DToizK5NlDc9291YzOw5YaWYvuPsT8RrG/jJYCNDc3Jz930aS0vA5K3n93QP8n55LDgt7gN52gBt6LGHkgQVAdCy/c6z+971H84XvLy5GySKSg5SB7+7nJdpnZq+bWX9332Fm/YGdCb6jNfZ7p5k9AgwD4ga+5F9n0HdqsF1x2zXYmwBMb7uS6W1X0sNg69wv8oWCVCkiQcv1ou0y4IrY6yuAX3VvYGZHmdkxna+BLwDP53hcycLU8EYGTXr0sLAHEt5I1eofOfT68hFNbJ37xbzWJyL5letF21uBJWb2DeBV4BIAM2sAfuzuo4HjgUcsOgjcA/i5u/82x+NKhkbNf5wXd74Xd9+8tnFxx/DntUXXstx2q4JepBLkFPju/ibw+TjbW4HRsdcvA5/M5TiSvXBLhCmPbOS9A+0J2yzrGAkH4YYeS2iwN2n1jzCvbRxvnvRltn3rrAJWKyL5pDttK1j3sfpklnWMZNmBkQCHxupFpLIo8CtQuCXCtQ+tz+qzl49oYraWRRCpSAr8CjM1vJGfrX4148/pBiqRyqfAryCZDOF0+nCvEBtmXpCnikSklCjwK8Bp03/LO/sTX5RNRL16keqiwC9j2Q7fNNbXcf35JzP2jMY8VCUipUqBX4bCLRG+99D6JOtYxvehkPHCnNF5qUlESl+ud9pKgXXOwMk07M8e3FdhL1Ll1MMvI9lMtzz+mJ6s0aMGRQQFflkIt0S4/hfrOZhht14XZUWkKwV+CQu3RJi8dAP7Mk16tP6NiBxJgV+isp2BoyEcEUlEgV+Ckq1smYhuoBKRVBT4JSTcEuGGXz7HgfbMHvSlsXoRSYcCvwSEWyJct2Q9GeY8Zw/uy2ItXywiaVLgF9mE+57mqZfeyugzBkzQqpYikiEFfhFNDW/MOOy1fLGIZEt32hZJuCWS8Swchb2I5EI9/AIKt0S4bcUWIrv3YRl8rrYGbrvkdC12JiI5UeAXQLzZN+lcn9XsGxEJUk5DOmZ2iZltMrMOM2tO0u4CM9tiZlvNbFIuxyw34ZYI331ofcZTLT/cK6SwF5FA5TqG/zxwEfBEogZmFgLuAi4EPgFcZmafyPG4ZaFzsbMMZ1ty9uC+uolKRAKX05COu28GMEs6Ij0M2OruL8faPgiMAf6Yy7FLXXTBs+fSbl9XG2LuRUM1Ti8ieVOIMfxG4LUu77cDwxM1NrOJwESApqam/FaWB1PDG3lgzWu0e/r9ej2BSkQKIWXgm9ljwAlxdk1x91+lcYx43f+EaejuC4GFAM3NzZmOhhRNuCXC9x/ewP629Fe21J2yIlJIKQPf3c/L8RjbgYFd3g8AWnP8zpIxNbyRn695lQ4tiyAiJa4QN16tBYaY2Ulm1hMYDywrwHHzrnMJ40zD/vIRTQp7ESm4XKdlfsXMtgNnAY+a2YrY9gYzWw7g7m3ANcAKYDOwxN035VZ28WWzXv2xvWu549LTdbesiBRFrrN0HgEeibO9FRjd5f1yYHkuxyolmYa9lkQQkVKgtXSy8MCa11I3ijl7cF+FvYiUBC2tkIV0plzW19Uy48unaqqliJQMBX4SXefUh8y4bPhAZo8dSsgsYej3rq3hjzdfWOBKRURSU+An0P3BJO3uh8btLxs+MO4Yfo3BLRedVrAaRUQyocDvJtwSYeZ/beLt9w/G3f/Amtd4aW70evTiNa/S2dHvXVvDLRedpiEcESlZCvwuwi0RJi/dyL6D7QnbdA7lzB47VBdjRaSsaJZOF7et2JI07AFCyReKExEpWQr8Llp370vZ5rLhA1O2EREpRVU7pNP5uMHW3ftoiK1W2VBfRyRJ6GtOvYiUs6oM/KnhjSxe/eqhJTsju/cxeelG/vFTjTz8TOSIYR3NqReRSlB1gR9uiRwW9p32HWznDy+8wdyLhh7R81fQi0glqJrA7xzCSTZk07p7H2PPaFTAi0hFqorA7z6Ek0hDfV1B6hERKYaKn6WTaAinOwOuP//kQpQkIlIUFR/4t63YklbYTxjRpKEcEaloFTek0326ZbIxe9ADxEWkelRU4HdfGiGyex9G/CemG/Bvl56uoBeRqlFRQzrxlkZwouHelYZwRKQaVVQPP9HSCE506EZz60WkmuUU+GZ2CTADOAUY5u7rErTbBrwLtANt7t6cy3ETSTRm31hfx1OTzs3HIUVEykauQzrPAxcBT6TR9u/d/fR8hT1Ep1XW1YYO21ZXG9J0SxERcuzhu/tmACuRJYM7h2m0NIKIyJEKNYbvwP81MwfudfeFiRqa2URgIkBTU1PGB9LSCCIi8aUMfDN7DDghzq4p7v6rNI9ztru3mtlxwEoze8Hd4w4Dxf4yWAjQ3Nyc6p4pERFJU8rAd/fzcj2Iu7fGfu80s0eAYaQ37i8iIgHJ+zx8MzvKzI7pfA18gejFXhERKaCcAt/MvmJm24GzgEfNbEVse4OZLY81Ox5YZWbPAf8PeNTdf5vLcUVEJHO5ztJ5BHgkzvZWYHTs9cvAJ3M5joiI5M7cS/e6qJm9AbySpEk/YFeBygmC6s0v1Zs/5VQrVHe9J7r7R+PtKOnAT8XM1uXzRq6gqd78Ur35U061gupNpKIWTxMRkcQU+CIiVaLcAz/hHbslSvXml+rNn3KqFVRvXGU9hi8iIukr9x6+iIikSYEvIlIlyirwzewSM9tkZh1mlnAKk5ltM7ONZrbezOI+lKUQMqj3AjPbYmZbzWxSIWvsVkdfM1tpZi/Gfh+boF3Rzm+qc2VRC2L7N5jZmYWsL049qeo9x8z2xM7lejObVow6u9SzyMx2mlnc5U9K6fymUWupnduBZvYHM9scy4V/jtMmv+fX3cvmh+iTtU4GHgeak7TbBvQrh3qBEPAS8DGgJ/Ac8Iki1TsPmBR7PQn4QSmd33TOFdE7vH9D9NHFI4A1Rfzvn0695wC/LlaNcWr+HHAm8HyC/aV0flPVWmrntj9wZuz1McCfCv3nt6x6+O6+2d23FLuOdKVZ7zBgq7u/7O4HgAeBMfmvLq4xwE9jr38KjC1SHYmkc67GAPd71Gqg3sz6F7rQmFL6b5sWjy5b/laSJiVzftOotaS4+w53fzb2+l1gM9D94R15Pb9lFfgZ6HzgyjOxB6qUskbgtS7vt3PkH4JCOd7dd0D0DydwXIJ2xTq/6ZyrUjqf6dZylpk9Z2a/MbNTC1Na1krp/KajJM+tmQ0CzgDWdNuV1/NbqCdepa3QD1zJVQD1xns+ZN7myiarN4OvKdj57Sadc1XQ85lCOrU8S3Ttk71mNhoIA0PyXVgOSun8plKS59bMjgYeBq5193e6747zkcDOb8kFvpfZA1cCqHc7MLDL+wFAa47fmVCyes3sdTPr7+47Yv+M3JngO4r1QJt0zlVBz2cKKWvp+j+8uy83s7vNrJ+7l+rCX6V0fpMqxXNrZrVEw36xuy+N0ySv57fihnSs/B64shYYYmYnmVlPYDywrEi1LAOuiL2+AjjiXyhFPr/pnKtlwNdisx1GAHs6h6mKIGW9ZnaCmVns9TCi/0++WfBK01dK5zepUju3sVp+Amx29/kJmuX3/Bb7ynWGV7m/QvRvwP3A68CK2PYGYHns9ceIzoZ4DthEdGilZOv1D67M/4nojI5i1vsR4HfAi7HffUvt/MY7V8BVwFWx1wbcFdu/kSSzuUqk3mti5/E5YDXwmSLX+wCwAzgY+7P7jVI9v2nUWmrndiTR4ZkNwPrYz+hCnl8trSAiUiUqbkhHRETiU+CLiFQJBb6ISJVQ4IuIVAkFvohIlVDgi4hUCQW+iEiV+P/cNVFFdopN5gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# my solver vs sklearn elastic net\n",
    "plt.scatter(coefs_enet[:, -1], beta_mat_en[:, -1])\n",
    "plt.scatter(coefs_enet[:, 50], beta_mat_en[:, 50])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f98f0a58310>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAaOUlEQVR4nO3de3RU5bkG8OedyYVENAFCjyaAWO1Si3BwEZQWbLvUiK0UQSuo0XqKLfaoRylK5GYSg1SEmlIqHrmorQtajQXHFK2AVbxQQEJHyaFILd5ggkvCRYqJ5DLv+WMyMQkzk9mZvWfPnnl+a7kkO+PMu2a1z3r59vt9W1QVRETkXC67CyAiotgwyImIHI5BTkTkcAxyIiKHY5ATETlcmh0fmpeXp4MHD7bjo4mIHGvHjh31qtq/63Vbgnzw4MGoqamx46OJiBxLRD4OdZ1LK0REDscgJyJyOAY5EZHDMciJiByOQU5E5HC2TK0QEaUaj9eHRev3oO5oI/JzszBj7LmYcGGBKe/NICcispjH68OstbVobG4FAPiONmLW2loAMCXMubRCRGSxRev3tId4UGNzKxat32PK+zPIiYgsVne00dB1oxjkREQWy8/NMnTdKAY5EZHFZow9F1np7k7XstLdmDH2XFPenzc7iYgsFryhyakVIiIHm3BhgWnB3RWXVoiIHI5BTkTkcAxyIiKHY5ATETkcb3YSEUVg5RkpZmGQExGFYfUZKWZhkBMRdRHswn0httAHz0hhkBMRJaiuXXgoZp2RYhbe7CQi6iDUSYVdmXVGilkY5EREHXTXbZt5RopZTAtyEXGLiFdE1pn1nkRE8Rap2y7IzcJD1wxNqPVxwNyO/G4Au018PyKiuAt3UuHiycOxeealCRfigElBLiIDAFwFYKUZ70dEZJcJFxbgoWuGoiA3C4LE7cI7MmtqZTGAEgCnhnuBiEwFMBUABg0aZNLHEhGZz8qTCq0Qc0cuIuMAfKaqOyK9TlWXq2qhqhb2798/1o8lIqI2ZnTkowGMF5EfAOgF4DQRWaWqN5nw3kREPeaE7fVmiLkjV9VZqjpAVQcDuB7AqwxxIrJbcGOP72gjFF9tr/d4fXaXZjrOkRNRUgq1sSe4vT7ZmLpFX1U3Adhk5nsSEfVEuI09iba93gzsyIkoKYXb2JNo2+vNwCAnoqQUbmNPom2vNwNPPyQix4lmGiX4cypMrTDIichRjDzswWkbe3qKQU5EjuC0hz3EE4OciBKeEx/2EE+82UlECc+JD3uIJwY5ESU8Jz7sIZ64tEJECSXUREp+blbItXEgcMxssk6jRItBTkQJI9xEyrUjCrBmh6/T8kpWujvhzwmPFwY5Edmuu4mU1947iIeuGZoSM+E9wSAnIltFO5GSKjPhPcGbnURkK06kxI5BTkS24kRK7BjkRGSrSN22Ex58nAgY5ERkq3CnFC6ePBybZ17KEI8Cb3YSka1S6ZRCqzDIicgSRh58zImU2DDIich0cz21WL31E2jbz5GOmqXYMciJyBSRNvUAPGrWSgxyIopZ1w48nFQ+atZKDHIiiknxii3o92E13syoQr7Uo07zsLBlEqr9Y056LTf2WINBTkQ9FgzxBekrkS1NAIABUo8F6SuBZnQKcwG4sccinCMnoh7xeH3YvPcwStKq2kM8KFuaUJJW1ela8ahBXB+3CDtyIjKk603NfKkP+bp8OQQg0IkXjxqEBycMjVeJKYdBTkTR2VmFhr+UYnzDpyjUfljoCqyD12keBoQI8zrtx4c+xAmDnIi6t7MK+PNdyG5uBKTzOvjClkmd1sgBoEEzcKCwBJvHX2pj0amDQU5EYc311OKP2/bh9fRZGODqPDoYXAcf07QEaAZK0qqQL4dwAP1QN6IEI8ffZlPVqYdBTkQhFa/Ygs17DwPofh282j8GO7KL2pdRuJASXwxyIjpJcCIlKNI6OJ+daT+OHxJRO4/Xh9ELXsW0Z9/pdH1hyyQ0aEanaw2agZUZNzHEEwA7ciICEPnZmdX+MZ3WwV25A5B9WSnKh02yoVLqikFOlMLmemqxausnUb222j8G1U1jcBNnwhMOg5woBXm8Psx5vhZfNEV+6HFXo8/uyxBPQAxyohQz11OLY2//AevTqpCfGfmQqyBu7ElsDHKiFOHx+lBevQvfOfFaVIdcAeBEikMwyImSXDDAjzY2AwBKMsIfclXd9FWQswt3DgY5URLzeH14fc1SrHM9276M0t3mHgC8oekwMc+Ri8hAEXlNRHaLyC4RuduMwogodtteeBzz3SswwFUPlwADXKFDHAhs7hEwxJ3IjI68BcA9qvp3ETkVwA4R2aiq/zDhvYkoBnf4/4BsV+dlFJcAfg38Oyh4yNWH46+Kc4VkhpiDXFUPADjQ9ud/i8huAAUAGOREcdR1pFAA7M0M3YELgP3+POTLIXwq/eDjIVeOZuoauYgMBnAhgG0hfjcVwFQAGDRokJkfS5Tyiio34f3Pvuh0TRH+jBQf8lAz8Q1MuLAA+QDy41MmWcS0s1ZEpDeANQCmqeqxrr9X1eWqWqiqhf379zfrY4lSWvGKLRg888WTQjwo3BkpB0aUcBoliZjSkYtIOgIhvlpV15rxnkQUWaguvKvgGSmzM57D17Qen0ke9o2YwWWUJBNzkIuIAHgCwG5VrYy9JCKKpOrJR3D5x5XYgONAJnBYe+OBlh+H3Zn5ol6CJeUPAQBOb/uHkosZHfloADcDqBWRd9quzVbVl0x4byJq4/H68Opzj+JX6cuRIS3t1/vJcSxKXxZyZyYA3HDxwHiWSTYwY2rlLQRughORRYoqN+H8+vWoTH8caeI/6feZ0nrSzkwA+MbXTuFMeArggyWIEtw5s17E+fXrsSB9ZcgQD+q4MxMIbOzZOP17FldHiYBb9IkS1PbqZTijZiH+mVEPP1wRQxwI7MwEAkfNrv7Zt+JRIiUIBjlRotlZhRNrf45CbYW0/Z3ZhcghfkLd+FXrJCyePJxjhSmIQU6UQJYuno/bjixCpmhUd55UA1Mry7KnYvHMMusLpITEICdKEPeVzcF8PIY00ahe36AZmNn8U1x63Z2YzS48pTHIiWw211OLK3bchgWuXZBuuvAWdcEFRZ32wwt9p2DJtDnxKZISGoOcyCbFK7ag34fVmJ/2BHq7TnQb4qrA9Oafo9o/BosnD8cd7MKpDYOcyAa/mDMLv017Gn3Sj3cb4EAgxN/0D8G2Uy7FR3OKrC+QHIVBThRP66ajdfuTqEzTqAPcD8Gq1suQc+0SbGMXTiEwyIniwOP1YeTz30a+HIU7yn3QqsDdzbfj0FnjORdOETHIiSxWVLkJ9x+ZjXzX0ai6cCAQ4ltwAZb88iFri6OkwCAnskhR5SbcdPi3eNn9V7hcxpZS/pp9Fa64b7X1RVJSYJATWeCcWS/i7fRb0cfdGHUX3qqCcvf/YF7pPFxhbXmUZBjkRCbaXr0M5+8oxfsZXwJA1CHepGnIuPZ/MW/YJAuro2TFICcyyfvlQ1Co+6MObyCwlHJMTkXOtZUAQ5x6iMfYEsXI4/XhcOnpOMdAiKsCjX43HutTgpzy/Qxxigk7cqIY3F9xP8pal8At0S+jtKhgYdYvMHtmGe6wtjxKEQxyoh6oevIRTPj4QVTAb6gL9wPtIU5kFgY5kUHHS/NwnTQbXguv01xsn/g3nlRIpuMaOVGUNjxcDC3LwSkGQlwVaFLBY31KUFDxMR/6QJZgR04Uhd2l30SR+Ax34dtkKEaVv8W1cLIUO3KiCJYuno8vy/rhPAMhrhq4oVkzYiFGlb9lbYFEYEdOFNKGh4tR1LAOtyP6aRRte7DP1n4T8a27foeRllVH1BmDnKiLN+7/NoqieFpPR6rAEc3CGxN3cB2c4o5BTtTG4/Wh79rrcImBEA924dtkKEY98BYmWFYdUXgMciIEbmZeLT7AZWwp5S0dgkMTn2MXTrZikFNKK16xBU/s+wHOE2Mbe07AjV7XPo5LuLWeEgCDnFLWG/d/G6tcuwAD2+tVgY3Z43hWOCUUBjmlnOBEyiUGl1EUwJ/OLMWkKfdYWh+RUQxySilflvZBkYFlFKAtxBVwVXwOLqRQIuKGIEoJ26uXwV+Wg0yDa+HB3Zmuis+tLZAoBuzIKem1luag0MA6OBA85KoPCio+wijrSiMyBTtySlpby8dAy3LgMngz069AzYiFKKj4yNL6iMzCjpySUlN5Di5W4134CbjQ64Ej3F5PjsIgp6RyuPQM9JEGpMP4GSnS/zz0unObZbURWYVBTknDX5qDPj1YC69xDcPIsjetK4zIYgxycrxPKoZiYOsnEINr4QAgI2/FyHGV1hVHFAemBLmIXAngNwDcAFaq6gIz3pcoonXTodufwEAY78KbxI3M8sOWlUYUTzEHuYi4ASwFUARgP4DtIlKtqv+I9b2JwjlSdgZy0WA4wIFAF57JLpySiBnjhxcB+JeqfqCqTQCeAXC1Ce9LdJK9T90GLcsxFOLBkcKZuBPywOcAQ5ySjBlLKwUA9nX4eT+Ai014X6JOvizrg6/D+Pb6Os1FQcXHeNi60ohsZUZHHur/VnrSi0SmikiNiNQcPHjQhI+lVOF58HpoWQ4yDYR4cHv9an8RCio+trZAIpuZ0ZHvBzCww88DANR1fZGqLgewHAAKCwtPCnqiUI6X5uFqaTa+sUddeHni/+EmPvCBUoAZHfl2AN8QkbNEJAPA9QCqTXhfSmGr7v8RtCwHpxgI8WAXXpsxHL0qjvCpPZQyYu7IVbVFRO4EsB6B8cMnVXVXzJVRyvKVnoli19EeHzU7zLrSiBKSKXPkqvoSgJfMeC9KXcGNPfk92NhTM2IhRo6/zbriiBIYd3aS7YpXbMFT+67EwB5srz+hLvSq4CFXlNp4jC3ZqqhyE57edyXSDXbhqsBzZ5aiV8URawskcgB25GSbprIcbAAMP/z4iGahb8WnfOwaURsGOcXdkXlnI7elvkdHzX4w+Hqc/ZNlltVG5EQMcoqrltIc5PbgZuYh6Yu88g9xtnWlETkWg5zi4kj5GcjVBrgNhnizAhkVnyPP2vKIHI03O8laO6vgL8tBrjZEfV548Gbmsd5nI4NPryfqFjtyssz75UNwju6Hy+BIYbALz7GuNKKkwo6cTDfXU4svS/vgHN1veKSwxjWMXTiRQezIyVTDyl7G23oDMsXYSYVHNBt9Kw5wYw9RD7AjJ1Nsr16GlrJcvIvJUYd4sAvfmD0OfSsOWF8kUZJiR04xe3rutbjZ/YrhkcJnUIQbHvgTrrCuNKKUwCCnHtv71G0466NncLPb2Ejhe1qA8yv+gRusLY8oZTDIqUeOlefj6/qF4S58Y/Y4XHHfausKI0pBDHIypHjFFszddyvOE2MhHnxiDx/2QGQ+BjlFxeP14fM1d+H37lfhNjiRss89CINKazHB0gqJUhenVqhbRZWb0Hftdfix+xWkGZhI8WvggQ+DSmutL5IohbEjp7CCXfh69ysQl4Ht9QCq07+PCXOf4Vw4URwwyCmkXy54AHc1PoZT3F8ankjZM3Ej18KJ4ohBTp0sXTwftx6pxCy0RN2BBwUnUs63rjwiCoFBTu2eK5+E23W9oQ786dbLUdYyBYsnD2cXTmQTBjkBO6twYu1/40caXRcOBEL8Tf8QvF9Yjo8mDLW2PiKKiEGewjxeH754/m7cKBuRKQAMduK3PLgG37G0QiKKBoM8RW14uBg/bFgHl4GHPQDAF8jES2feh1um3GNtgUQUNQZ5CtpaPgZFWmt4GWXDiGV4cMJQPr2eKMEwyFOIx+tD37XX4RLXLkMh/pZ/CL4z729cRiFKUNzZmSKCuzOjDXFVoNHvxmN9SnDJvL9ZXyAR9Rg78iTn8fpwbM1deNn9V7hc2m2IB3dmrmq9HKdduwR3cKSQKOExyJOUx+vD5ucfQ5msjGp3pirwBXphdvMU7M4bi43TvxeXOokodgzyJFRUuQnn16/HgvSVyJambl8fvJn5X81zUDl5OJawCydyFAZ5khlW9jKOnWjFUxlVhkJ8w4hl+IAbe4gciUGeLH4/Hvrh63gXADKj+09Ugc16Ab4zbzMnUogcjEGeBPaWX4Cv677AOnjbWnjHw6y6UgWOaybeK5yHMeNvi0uNRGQdBrmDebw+vPrco/hN+r6TbmaKBAK743VVwA/gzzwrnCipMMgdaK6nFse3/wH3uqtwdXp9xImU/f485Msh1Gk/LGyZhNMuuhEPci2cKKkwyB3mvrI5KNHfoW/a8ag29oxpWgIAGH12X6z+2bcsro6I7MAgd4qdVfh87XQswL8hUezHDU6juASonMSzwomSGYPcAZYuno+fHFmMHGmK6qjZ4CPXlg16BB+wCydKejGdtSIii0TkPRHZKSLPi0iuSXVRm4vnb8TVh5+MeiZ8vz8Pdzffjj0TN3IphShFxNqRbwQwS1VbRORhALMA3Bd7WVS8Ygs27z0MAMjPrO/29Q2agdktP0XvkTdiCW9mEqWUmIJcVTd0+HErgB/FVg55vD5Me/adTtfqNA8DJHSYqwJH0Bt7R5RiMWfCiVKSmcfYTgHwl3C/FJGpIlIjIjUHDx408WOTg8frw/AHNpwU4gCwsGUSGjSj0zVV4JC/N6a33oE3JryNkQxxopTVbUcuIq8AOD3Er+ao6gttr5kDoAXA6nDvo6rLASwHgMLCwgj7DlOPx+vD62uWYp3rWeRn1qNO87CwZRKq/WMAIPDvZqAkrap9JvyFvlNwx7Q5+LXNtROR/boNclW9PNLvReQWAOMAXKYaaWM4heLx+vDanx7FQ2lfnVQ4QOqxIH0l0IxOYV7dFPjzf5yagW3TimyrmYgSS0xr5CJyJQI3N7+rqg3mlJQ6iio34f3PvsBbIU4qzJYmlKRVtYd30De+dgrPCieiTmKdWnkUgbP2Nkpgm+FWVf15zFUlubmeWqza+kn7z/lhbmTmy6FOP980ahC31xPRSWKdWjnHrEJSgcfrw6y1O9HY7O90PdxUSp32A8AunIgi487OOPF4fZhe9Q78Ie4iLGyZdNLTfBo0AwcKS/DR+KviWCURORGDPA5CzYZ31HUq5VPpB9+IEo4UElFUGOQWCyyn1Hb7uuBUSnAZJT8OtRFRcmCQW8Dj9WHR+j2oO9oIlwhao5zKPC3TzbVwIjKMQW4ij9eH8updONrY3H4t2hDneeFE1FMMcpMEl1Aam1sN/XccKSSiWDHITeDx+nBP1btRd98A4BLgxosZ4kQUOwZ5DwXXwX1HGyEAuotwtwj8qsjPzcKMsefyiT1EZBoGeQ94vD7MeO5dNLcNhXcX4lnpbjx0zVCGNxFZgkHeA+XVu9pDvDt9stNR9sMhDHEisgyDvAc6TqWE4xbBI5P+kwFORJYz88ES1CYr3c0QJ6K4YZD3QJ/s9LC/K8jN4no4EcUVg7wHyn44BOlu6XQt3S1YPHk4Ns+8lCFORHHFNfIu5npq8cdt+9CqCrcIbrh44Emz3sGgDm7D50ghEdmJQd5B8Yot2Lz3cPvPrartD4AIFeYMbiJKBFxawVdPsO8Y4h39cdu+OFdERBS9lO/I53pqsXrrJxE39RjZek9EFG8pGeQdt9dHwy3S/YuIiGySckHek1MKb7h4oIUVERHFJuXWyBet32MoxEef3ZcnFBJRQku5jrwuyuUUnpFCRE6RckGen5sVcW1cABTzYQ9E5CBJG+Qdn5vZccPOjLHnnrRGHjxPvIAbe4jIgZIyyLuOFPqONrY/yZ67Moko2SRdkHu8vpBz4Y3NrVi0fk/7jkwGNxEli6QI8o7LKC6RsJt7or3RSUTkJI4P8q5z4ZF2YebnZsWrLCKiuHH8HHm0c+ECYMbYc60viIgozhwf5NEslwRHCrkuTkTJyPFLK+Hmwt0i8KtyKoWIkp5jgtzIXHhWupuPWyOilOGIIO96Q5Nz4UREX3FEkIe6ocm5cCKiAEfc7Ax3Q5Nz4UREDgnycPPfnAsnInJIkM8Yey6y0t2drmWluzkXTkQEh6yR84YmEVF4pgS5iNwLYBGA/qpab8Z7dsUbmkREocW8tCIiAwEUAfgk9nKIiMgoM9bIfw2gBAh76CAREVkopiAXkfEAfKr6bhSvnSoiNSJSc/DgwVg+loiIOuh2jVxEXgFweohfzQEwG8AV0XyQqi4HsBwACgsL2b0TEZmk2yBX1ctDXReRoQDOAvCuiADAAAB/F5GLVPVTU6skIqKwRCM8iMHQG4l8BKAwmqkVETkI4OMOl/IAWDLtkiT4/UTG7ycyfj/dc8p3dKaq9u960ZY58q6FiEiNqhbaUYsT8PuJjN9PZPx+uuf078i0IFfVwWa9FxERRc8RW/SJiCi8RAny5XYXkOD4/UTG7ycyfj/dc/R3ZNrNTiIiskeidORERNRDDHIiIodLqCAXkXtFREUkz+5aEo2ILBKR90Rkp4g8LyK5dteUCETkShHZIyL/EpGZdteTSERkoIi8JiK7RWSXiNxtd02JSETcIuIVkXV219JTCRPkPEWxWxsBXKCqwwD8E8Asm+uxnYi4ASwF8H0A3wRwg4h8096qEkoLgHtU9XwAowDcwe8npLsB7La7iFgkTJCDpyhGpKobVLWl7cetCByJkOouAvAvVf1AVZsAPAPgaptrShiqekBV/972538jEFY81L8DERkA4CoAK+2uJRYJEeRGTlEkAMAUAH+xu4gEUABgX4ef94NBFZKIDAZwIYBtNpeSaBYj0ED6ba4jJnHbom/WKYrJLNJ3pKovtL1mDgJ/ZV4dz9oSlIS4xr/RdSEivQGsATBNVY/ZXU+iEJFxAD5T1R0i8j2by4lJ3IKcpyh2L9x3FCQitwAYB+Ay5QYAINCBD+zw8wAAdTbVkpBEJB2BEF+tqmvtrifBjAYwXkR+AKAXgNNEZJWq3mRzXYYl3IYgI6cophIRuRJAJYDvqiqfzAFARNIQuPF7GQAfgO0AblTVXbYWliAk0Bn9HsBhVZ1mczkJra0jv1dVx9lcSo8kxBo5ReVRAKcC2Cgi74jI43YXZLe2m793AliPwI28KoZ4J6MB3Azg0rb/zbzT1n1Skkm4jpyIiIxhR05E5HAMciIih2OQExE5HIOciMjhGORERA7HICcicjgGORGRw/0/2iK+3b9Ox7MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# my solver vs sklearn lasso\n",
    "plt.scatter(coefs_lasso[:, -1], beta_mat_lasso[:, -1])\n",
    "plt.scatter(coefs_lasso[:, 50], beta_mat_lasso[:, 50])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f98f0a92a50>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAl7klEQVR4nO3dfZQcdZkv8O/T85J0iEwwM9yYSSAcl8MuC7lMMuPRTe51l2wEFx0jrCNycfW6bPDg3qCHOzFRNi++JSS7BLPq0dzo8XDXdRlMCMMqJ7Li6p0ccTOTSAKyLCKLyUxiXmACMg3z0s/9o7sm/VJVXdVd1fX2/ZyjIdUzVdUveerXTz2/5yeqCiIiiq5U0CdARES1YSAnIoo4BnIioohjICciijgGciKiiGsM4qCtra26aNGiIA5NRBRZQ0NDZ1S1rXR7IIF80aJFGBwcDOLQRESRJSIvmm1naoWIKOIYyImIIo6BnIgo4hjIiYgijoGciCjiAqlaISKysu/wMLbvfxYjoxnMn5NG73VXYFVHe9CnFWoM5EQUGvsOD2P93qPITEwBAIZHM1i/9ygAMJjbYGqFiEJj+/5np4O4ITMxhe37nw3ojKKBgZyIQmNkNONqO+UwkBNRaMyfk3a1nXIYyIkoNHqvuwLppoaibemmBvRed0VAZxQNvNlJRKFh3NBk1Yo7DOREFCqrOtoZuF1iaoWIKOIYyImIIo6BnIgo4hjIiYgijjc7iSix4tLXhYGciBIpTn1dGMgp0eIyIiP37Pq6RO0zwEBOiRWnERm5F6e+LrzZSYnFTnvJFqe+LgzklFhxGpGRe3Hq6+JZIBeRBhE5LCL/7NU+ifwUpxEZubeqox1bbrwa7XPSEADtc9LYcuPVkUyreZkjvxPAMwAu9HCfRL7pve6Kohw5EN0RGVUnLn1dPBmRi8gCADcA2O3F/ojqIU4jMko2r0bk9wFYC+BNVj8gIqsBrAaASy65xKPDEtUmLiMySraaR+Qi8h4Ap1R1yO7nVHWXqnaqamdbW1uthyUiojwvRuTLAHSLyJ8BmAngQhH5B1W91YN9E/mKE4IoDmoekavqelVdoKqLANwM4HEGcYoCY0LQ8GgGivMTgvYdHg761IhcYR05JRYnBFFceDpFX1X/FcC/erlPIr9wQhDFBUfklFicEERxwUBOiRWnKdqUbAzklFirOtpx09J2NIgAABpEcNNS1pVT9DCQU2LtOzyMPUPDmFIFAEypYs/QMKtWKHIYyCmxWLVCccFATok1bFGdYrWdKKwYyCmxjNy40+1EYcVAToll5MadbicKKwZySqx2i3pxq+1EYcVATon1J79v3oXTajtRWHk6RZ8oSn7876ddbY8SdnVMFgZySqy49loxujoapZVGV0cADOYxxdQKJVZce62wPj55GMgpseLaayWu3zTIGgM5JVZcF1+O6zcNssYcOSVaHBdf7r3uiqIcORCPbxpkjYGcKGaMCxOrVpKDgZwSLa5lenH8pkHWGMgpsfYdHkbv957ExFRuSv7waAa933sSAMv0KFp4s5MSa/MjT08HccPElGLzI08HdEZE1WEgp8R6eWzC1XaisGJqhcjEsq2Pxy5vHhlH+oAffQ44dxxoWQCs2AAs7gn6rEKNgZzIhLG4RBynt4f6Bu+RPuCRNcBEfvLSuWO5vwMM5jaYWiGqIE7T240+LMOjGSjOX6hCs07pjz53PogbJjK57WSJgZzIgbhMbw99H5Zzx91tJwAM5JRgbhZ0i8v09tD3YWlZ4G57rY70ATuuAjbNyf15pM+f4/iMgZwSa2aTs49/nKa3h74Py4oNQFPJuTSlc9u9ZuTjzx0DoOfz8REM5gzklFivT2QtH4tbIy2DXx0f9x0exrKtj+Oydd/Hsq2PV59zX9wDvHcn0LIQgOT+fO9Of250xigfz6oVSqyWdBNGM+U143PSTTiw7toAzsh/fvRh8Xwhi8U99alQiVE+noGcEmt8csrV9rjwug+L3Q3UUH+TaVmQT6uYbI8YplYoscYsUitW28lc6G+gWqlnPt5nDOREVJPQ30C1Us98vM+YWqHEEgFUzbeTc5FeyKJe+XifMZBTYpkFcbvtSeN0Kj8XsggeAzlR0jhoSuW2EoULWQSLOXKiJHE4CSb0U/mpCAM5UZI4nAQT2UqUhGIgp0TrTg1goHkNfj3jFgw0r0F3aiDoU/KXw0kwka1ESaiaA7mILBSRH4vIMyLytIjc6cWJEfmtOzWArU27sSB1BikBFqTOYGvT7ngHc4dNqfyayk/+8GJEPgngLlX9AwBvB/AJEbnSg/0S+WptYx9myXjRtlkyjrWN0Wua5JjDSTCrOtqx5carY9tzJm5qrlpR1RMATuT/+1UReQZAO4Bf1rpvIj/NlzMW28/W+UzqyKhOcbCUGitRosPT8kMRWQSgA8DPTR5bDWA1AFxyySVeHpaoKifQinaUB/MTmItYh6+YTIKh8zy72SkiswHsAfBJVX2l9HFV3aWqnara2dbW5tVhiap2z0QPxrS5aNuYNuOeCQa5MPCsNW4CeBLIRaQJuSD+HVXd68U+ifzWn12OdRO34Xi2FVkVHM+2Yt3EbejPLg/61BIv9GuLhkzNqRUREQDfBPCMqt5b+ykRBYHz8svkZ4DqueP4LVqxZfwDGLxwZV2m30e2NW5AvMiRLwPwYQBHReQX+W2fUdUfeLBvoppZ9Qwxyg+NypUFkis/xAQA3BDoOQfOmAE6kYEAmIfT2NK0G+teAdbvzb1efgZUTkhyp+bUiqoOqKqo6mJVvSb/PwZxCgW7r+iJLD90ymQGqPHa1GOqPickucOmWeQ7p130/GD3Ff3/WZYfnsGyrY8nu5OfxQxQozTT75FxpFvjBoCBnHzl+XqOLtl9Rc/OSCGF8tWAskhhOP979T5fTzjobliRxTJoIzoXgP8jY7bGdYe9VshXQXfRs/uK3mASxAGUbY9E178jfcCOq4BNLcDe1RW7G1ZkMgN0TJuxbbKnbiPjVQ0HcGDGGrww83/gwIw1WNVwwPdjnn8d5+T+dPu6BYSBnHwV9E0ru54hw9pq+jtm20N9k62oNS1QVoFj0t2wooJl0BSCk2jD+onbMHThyvpM1XfYbjfyx/QIUyvkq/lz0tNpitLt9WD3FX3Ngz3426ZdaJbJ6Z8f10ZsmyxPQ4T6JptZa9pSVl0P7eRngOaqVoAvV3Nu1bJrt+vXrNQgjukRBnLyVRhuWtn1DBFM2f4diMBNNidBOn2R/+fhJYftdiN/TI8wtUK+CnMXvY2N96NJitMQTaLY2Hh/oOfremq6VWvaQuO/i0SKYJrDdruRP6ZHRANYabazs1MHBwfrflyiQrqxBSIm2xWQzefqf0Ior/IBct8IbC8mBZN3bLUsBD71lIdn6yOz59SUzuXt/UpzBHFMl0RkSFU7S7dzRE4UIlVV+RTcmATEutlABFIE00qeE1oW+h9Q/T6mjxUxzJEThUjVVT75G5P7Dg+ja99/R7vZZKcIpAiKBNFu169jlo72jYoY45g14oicKEQqTU2vlD/fvv9Z0/a8GcwoWwWI6sjhotfV4oicKECl7Qv+5PfbsGdo2LTKp9Is2X2HhzE8msEwlgMTuaXs5stZjOhcbJvswc6Q5HkTyeeKGAZyooAYgXnl1E/wQHMf5mfO4MThViy+7A58+VRHWd37sq2P2+bPjaAO5Hqt94+f76suAK49PByKaqFEsmh54FW6i4GcyAPVNAbbvv9ZrJz6SVEr3XacwXt/cw96bvzKdO7USKeYTawCcvlzs5ukhTR/PAZye741eFuxwbwixqN0FwM5UY2qbQw2MprBA83lrXTTeAN46OPA3tUYS8/DwGs3YXj8jyz3M39O2lELgarbDHjRhCsCfG3w5mLR62rwZidRjaptDDZ/ThrzLVrpQqcAKGZlTuBzsgvdqQHTHzPy505aCFTVZiDC/Ufc8r3B2+KeXB3/ptHcnx5eDBnIiWpUbclg73VX4ATMG3cVKlzsojs1gIHmNfj1jFvwxMw7cX/Xi1jV0W7aHKxQ1W0GfK62CJOgG7zVgqkVohpV2xhsVUc71j98C/5Gv16WXinbl5wtW5puHk5j3tGNwKKLsKojN7oz8rst6SaIAKNjE7XleiPcf8SKVR7c7wZvfi6wwkBOVKNaGoONjU/i9aZmpDUXnLMQNEj53MwRnWu6NF1hdz675mBV87naot7s8uB+Nnjze4EVBnJKlMJR0a9neLPPqlezOdKHrc3fzN3czJvUFKZUilrrTjbMxO7GW7Fh4j7T3ei541hew9J0tiNFn6st6s2oFFrb3If5cgYj2optkz3Yvr8ZB9ZdO/0zXo+a7fLvDORELpg1pLK04ypX1QVVjYZ/9LmiIA4AM2QKL+lszEy/CbMyJ4GWBWhcsQGbFvcAO/ZYLr9W7dJ0FUeKPldb1FvnK49hS0F6aoGcwdam3Vj/CgBc68+3Gviff2cgp8SoVGtdxAiYHvfEKD6GeZ75InkN8mmT1rUmo+MMZuCeieLzcjPSK3xNulMD+dmgZ3Dq4Tag4Uvne49ENHCXWt/8IGahOD01S8axvvlBAFt8O67f+XdWrVBiOB39lLW2NanScN0z3IxFnlms8s8m3fnWjf8l+rPLy37UavJQKeM1MW6kLkidQUpyN1JdlRlGZK3L/wLzck+r7V6xW3LQCwzkFAleBM5aRj/Zc8enj2mkI4ZHM1CcT0e4PieTBY4r5p9LapEHL1xZVJI40LwG3akBSP48KzFeE9sbqZVEqNbc6iJpefH06ALl9wIrTK1Q6JnlcX+y56t428MPYJ6ewSlpw7Elvejqvt12P2ZVCU6NZOdO5449u3HlQf75viufwx8Oled8MQFs35+ueD7Ga2I5MclJmWGU1rp0c/P2SB8mH/5faJx6Pff3c8dyfweqel5+5d8BBnKKgNLA2Z0awBcbdudynfk0QMvQ3TgI2Abz0uoSp8a0Gdsme5DJTtn+rtV226qQGvPPXc//PSDlOd+1jX34b6PlKZdSxnmcergtl04p5aTMMEq15i4unmOPbsAsI4jnNU69ntsesgsUAzmFXmmANEsDpGUcCw9tByqMyguDOSrEclVgOF+eZuShjWDs9MaV3/XDVsFyvpx1nEpa1XAASE9BM7kuiYbJhplodFJmGLVac4cXz5mZk662B4k5cgq90oBklQa4WE1GlCUK89t2MpiBOyfuwPLxnUU3E40RtdMbV77377AIlicw19mNNCO/nXlpOoirAi/pbKybuA37ppZV3kc1uf4IGMnOdbU9SAzkFHqlgXNELfqTCCrejHJSgqgKPLXk83is4Z1F241g7ebGlVm6pTs1gAfG/sr5DTS7G24mQTSDGRhZutbZiN8kvy0CjOlMfG/8j5xdcIJYX7MOdjffWrbS0pg2Y3fzrQGdkTWmVhLKz74PXivNbX81dQu+qDuRKikTTAEV85dOc+PDC9+DLQutZ/k5vXFVmoYp7ZdSsU690lqPJjnf9IoN6HIaRG1SM4CLCSsxqjU3XHPDamx4aBKf1H+aXmnpPtyM5TesDvrUyoiq5Zrbvuns7NTBwcG6H5dyzGY4ppsaPC2H8ptuakFpuTcAZFXQv+ppy+dRuEDDCzNuKa8ZR25EfmX2AU9ej9LXeqB5DRakzBZGXpgrKyy14yqL/LPFz7tlsf/j2VYsH9+J9jnp6anrSRS2AY+IDKlqZ+l2jsgTyO++D2V8WJhAWhZaTle3ex73Xfkc5g9tw1sqTADJTExh8yPWFwSnyvqwpM6a/6Dbyg+z7dW8zibleEaVjpcTVqLKz5JBLzGQh5wfI4K69l2ulBqo1ooNGNvziaLqFSMAWT6PI33oOroREGfP8+WxCezzYJ3LomCww2WFh9OKkGpf54LUjJ47himkMBPj+Ezzg/jwkkXo6rje+ncpNHizM8Q8m0FYwqoszYu+D6UzMMce3eDPwgSLe7Ct6Q4cz7Yiq4Lj2Vasm7gN/dnl1s/DbOJKBVVVl7i8OTmmzdj02k3m76vTipBaFoBY3AOs2ABpSqMR2ekp+l1HN4ZydiaV44g8xPxKgfjVd9msZnrmjJMwTWZ7MFnkmhtWY+XedyAz7vB5uDjmQPOaXJvTsVbgyBbn3x4c3pwce3QDZo6dxIjOzdWpv/E2/PlDX8O7frhnuuthUWqkUsqk1kk5UZqdSWUYyEPMrxRI1f2zKzC78IzoXCwwq/v2YLKI6+dhlaYwYdyQXCBn3KWCnATExT1Y+YNWDL9RXM3yOdmNWRmLapZKx651Uk6UZmdSGQbyEPOz9aUfN3HMLjDbJnuKy+0ATyeLuHoeZn02TFh2P3QSyC0D4rFcqiU/oh4ZvaDoYcumVQ99HNi7uvLNy1oXgIja7Ewq4kmOXESuF5FnReRXIrLOi32S/60vvWZ2genPLsfn5eM4iTZkVXASbTh49eZgvq4XTFzJQuCq8tbpyNQ28J3vDPiR2f9W9Ihl0yqdgqOOgrVOyonp7MykqLmOXEQaAPwHgJUAjgM4COBDqvpLq99hHblzfd/6Oyx78Wt4C87gBFpx4NI70POxu0x/9mD/N7Dw0HZcrKfLOgKWPvbapSvw1tEDjkrV7PZrPH75oc+jRV8FACiK0+JvoAFNyCIFxRRS+JleibemTjrqXPiznR9F19mH0YBsboOc3/frksZkNosL8qvsZCFIQfHb0n2WluVd/i6MPf0DpMdOmNaRW0q/GXj3PbnXqWCfbzS1YHIig1n6xvTJVdrtCFqxdbwHGxvvx5vld7nfcXIuLQtzFTuPbsDMzEmMZOdid/OtuOaG1cXfTGxKEa3ez/Pbz+CUtOLYkl4AsH3vzVT6vFD1rOrIvQjk7wCwSVWvy/99PQCoquVyGwzkNooCxIWQ8deK1m/MaDOeWvqFsn8YB/u/gauG7ka64Ou58bMAyh5TLQkcTWnTEZzdfo1//P916DNF51hJ6bGtntPPdn4Ubz/7kLtgW7rPRRc5Sqc41tAMdHwYePIfa9qnKjCOBswQdy11FcBUw8zzrVUBZBX4rq7EBe//ci6Yl95wBabf34P/+bLp+/nIpeuw8YU/LLrH8f7GA/hiw/8pSvlYvVeGSp8Xqo1VIPcitdIOoDC5djy/jdwqadA/Y+JcWYCc7vJXYuGh7UX/eAp/1uwx0zzw3r8qK5ez26/xuJsgbnZsq+fUdfbhqoJ40T6rKDm0NTUODH275n1OIeU6iBu/11jSWjUlwIfkMfzi+7tyG2xuuFq9n8te/FrZjeq7Ug9Yd5m0UOnzQv7wIpCb/VMrG+aLyGoRGRSRwdOnK3epSySHQediLc+nWnX+u1jPOOoKOK0kF2u3X7vH3TJ7TtPplFr26UfVhboPwIXGtBmpKp6bKtCg5r+XEuC28X/I/cWmAsXq/XoLymecWneZtJ4VW+nzQv7wIpAfB7Cw4O8LAIyU/pCq7lLVTlXtbGtr8+CwMeQw6JyS8u5/p8T8NT0lrZaPWSqYSGK3X7vH3TJ7TlM1fjxPSas/VRdS3XkpgJNow7qJ26w7ONp4GbPxW5vXe3r6v80sUav36wTKW7NanaPZe3X+MfvPC/nDi0B+EMDlInKZiDQDuBlAvwf7DYd6LirrIOhktHn6JlShY0t6kSlpuWn8rNljFW+N5C8qdvs1Hh9Xd1Wspce2ek4H577PXWVJgazmzs20GqNWjWnX+xzXRgwu2YYn3vcTPNbwTmyb7MEb2lD5Fwt+/0v6URxb0lv+dTfv9fS83H/YVKBYvZ8HLr2jrELq77IfLGvjavVeGSp9XsgfNdeRq+qkiPw1gP0AGgB8S1WfrvnMwsCvPiFWTGqBp6QRr2oaF+rvcpUES80rALq6b8dBoLjqoOBnSx97bZFRtWIxQSZ/Uam0X+Px3zv0ecwxqlYKKktQcmNTAbwy+63IvPaK6f4KvWPNt/GzncDbzj5kOuJ4DTPQiCnMwGTRMbIAXlh0c/E+S6pW8NwPoaPHLLsfqgBSekPYMDEG3LjLsmrldUljciqL2ZKrphnFm/CrpX+Dru7b0ZXfxfb9zVj7CrC5+f+iBa9CAGj+/0sj9cuYjZ1Nt2H5DavR1dEOpJ6FDn6zbDWfWe/OT8e3mQ3atbj8s3BsaS96um9Hc0lfn3de9wk8fexSy/feTKXPC/mDbWzt+N1C1IybDnZedBW0qXBwuq/SqfndqQHc07S77KZXUfleDec3ps3TfVUA4KOz/w2bLtjj+nXQjS2WgVw2nwvm/XfKh46SFH5sY1uNIKYtO23Q79W3BQ9Wci+cmt+dGsC9TV9Ho5jclGu+wH2wMbkBbCwu3D++HOmmBlxzw2qgY7O7/TpR62xJP8VwIQeqHgO5nTBPW/ayyVGNQcGYmm+sfmMaxIHqLoA2K9i0+93o34OLHFE9MJDbCfOIzCLAZc8dx1vXfb+uq5kYPWFM+4UUquYCaHExTc1ZgAOfqsPKNQGNfMO2Mg2FG/uR2wnzorIWQXEkO9fT3uVOGD1hLPuFANVfABPYA8SvPvQUXxyRVxLWXKTNEl0Gp73Lax39GT976uE2zIPJhBBpqP4CmKD0hvE+mHW89HUpPoo8BvKQcRxUSwLc8Wx+gYJ8JYehUu9ys8Ug1u89CgCug/nBY71oMeuzseQLzld1NxPWi6mHzBbELuXLUnwUCwzkIeI6qBYEuA8WrA5fqFLvci9XIfrkLy/H0onbsLaxD/Pl7PTqN0O/vBwHul3tKnHM3odSxnvJ/DmVikwgD+2H18N63lqCarXLt3m5CtHIaAbDWI7+8eJvBcKRZEVmF+FCxnvp1TcoipdIBPIgP7y2FxCPZ37WElSrXb7Ny1WI/FzRKM72HR425nSaKiyzXLb1cV/WcaVoi0Qg92sR4koqXkA8XrC21kBYzfJtXi7E7NeiznG3ff+zlkH81rdfgi+sunr6736t40rRFonyw6A+vHYXEACez/wMYmm3VR3t2HLj1Wifk4YgN/rbcuPVVV0gvdxXkth9jvcMDReVHVpd1PmtJ9kiMSIP6it7xQuIxzM//Vrd3slxvTqG2b5Ce38jJKw+30D5N09+6yEzkQjkQX14K15AfJj56cfq9kHizbnKzD7fhQoHFEFd7CncIhHIg/rwVryARGGySsBd8oK6vxElxutwV9+TmDLpRlr6zTNuF3uqXSQCORDMh9fRBSTMk1Xq3U/dBG/OOWN8ppg2oWpEJpAHJdKjH4+rapwqzImnRByNMolpE6oeA3mcBdBPvTQnbhbEOcq0FumBAwWGgTzsqs1xH+nLLRJstuK7j/3UraaaN4ggq8pRJpEPGMjDrNoct/F7ZkHc5xawVrnvrCpe2HqDb8clSrJITAhKLLsct9vfA2prJ+tQ3Cas7Ds8jGVbH8dl676PZVsfZ09wCiUG8jCrNsdt9bhmfa9WCWJ2ql+4wANFBQN5mFnlsivluKv9PQ/EaZp+xRYNRCHBHHmYVTtzNOC1Rn2rvKjz5CbWwFNUcEQeZtWuGRrmtUarZdzAPXcMgJ6/8Xukz7dDxi3fT/HFEXnYVTtzNMwzTqsRwOQmNqiiqGAgT7BIdSUMYHITZ1pSVDCQJ1TkuhJ63DLYKc60pChgjjxOjvQBO64CNs3J/WmTP45cRcaKDbkbtoXqeAOXKMw4Io8Ll7NAI1eREYWWwUQBYSCPC5c3A/1edcmP/Pu+qWXY/sZOjLyewfyZafROXYFVnpwtUbQxtRIXLm8G+jkD048ZkZxlSWQtOoHcLv/rIjccWy5nc67qaMf9XS/iiZl34tczbsETM+/E/V0venJjz4/8e1A5ffZaoSiIRmrFLv8LBL4KTqGD/d/AwkPbcbGexilpw7Elvejqvt3/A7udzXmkD11HNwLIAALMw2nMO7oRWHRRza+bH/l3u31Wm8Z5GbPxZvzOYnsEK3sosaIxIrfL/1bbIdAHB/u/gauG7sY8nEYqHxyvGrobB/u/4f/B3c7m9PF182NGpNXvzpnVVHXKZdPEX2Bci8cy49qITRN/AaA+3wI44icvRCOQ2+V/A5goYmXhoe1Iy3jRtrSMY+Gh7fU5gcU9wKeeAjaN5v60G1n7+Lr5kX+32qcqqg62/dnl+O7UH2NSU1AFJjWF7079MfqzywH4X9nDvD95JRqB3C7/G2Cnv1IX62mL7WfCl8f38XXzowPiqo523LS0HQ0iAHIrDt20tB3nMhOmP+8k2HanBvCBhp+iUbIQARoliw80/BTdqQEA/vdaiVwtP4VWNAK53WSQEE0UOSVtpttfkdl1b/hUkc+v26qOdhxYdy1e2HoDDqy7tvbSw8PD2DM0PL0G6JQq9gwNoyXdZPrzToLt2sY+zCr5BjVLxrG2Mfe++N1bPXK1/BRa0QjkdvnfxT04ePVmnEQbsio4iTYcvHpzIDc6jy3pRUabi7ZltDkXDEKSx58WsQ6JVqNXEVQdbOfLGdvtfvdWZ3dF8kpNVSsish3AewGMA3gewP9U1VEPzqucRTe/fYeHsf7gpchMfHl6W/pgA7YsHK57ZUFX9+04COSrVs7glLTi2NJedB36tPkvBJDHLxKhDolWo9TRsQns+OA1VVWtZJFCClmL7Tl+9lphd0XySq3lh48BWK+qkyJyD4D1ACyilj/s8oxBlIh1dd8O5MsN5+X/h+f/PpCGT0DEOhzasJuJWm2wbTAJ4nbbvcbuiuSVmgK5qv6w4K9PAPjz2k7HvUjkGQNasSdOddB+jF6HtRULTNIrw9qKet0qZ3dF8oKXOfKPAXjU6kERWS0igyIyePq0eXVHNSKRZwwoHx2nqgg/8tXbJnswVnJPY0ybsW0yGukmIkPFEbmI/AvyGYISn1XVh/M/81kAkwC+Y7UfVd0FYBcAdHZ2alVna8JqpHbflc8BO9aEp1NeAPnoSHxbccHr0Wt/djkwkatemS9nMaJzsW2yB/3Z5djp2VGI/FcxkKvqn9o9LiIfAfAeACtU1bMA7ZRZnvG+K5/LTT8PybT9oPjd4TAO+rPL0T++POjTIKpJrVUr1yN3c/OdqjrmzSm5VzZS27Gm7us7hhGrIoiSodaqla8AmAHgMcnNuHtCVT9e81nVKkTT9oPEqgiiZKi1auX3vDoRTwW0vmMYsSqCKP6iMbPTrRBN2/cTO+cRERCVfuRuJWB9xzjViBNRbeIZyIFITT+vRthmtBJRcOKZWkmAuNWIE1H1Yjsij0uPESusESciQyxH5ElYecXvXtlJkF+jwvF2orCKZSCPU48RK373yk4Cq3nI9Z+fTFSbWKZWwpg/9iPVwxrx2lw0qwkvj5UvFXfRLPNVh4jCKpYj8rB1RExCqieKOCKnuIhlIA9b/jgJqZ4oGrVYuNlqO1FYxTK1ErYeI2FM9RDQIDK9mHPpdqIoiWUgB8KVP2apYDiZBXG77URhFcvUStiELdVDOe0WF1Kr7URhxUBeBywVDCdeYCkuYptaCZswpXooZ1VHOx4c/A0OPP/S9LYll7TwfaLI4YicEuvufUeLgjgAHHj+Jdy972hAZ0RUHQbyGrEneHR95+e/cbWdKKyYWqkBe4LXh18N0DghiOKCgbwGYe4JHpfuj7xYElXG1EoNwjrRJ04tAfycFTuryfzjb7WdKKz4ia1B2Hq6GOLUEsDPi+WXblyMVMkkzpTkthNFCQN5DcJahxzWbwrV8PNiuaqjHff2XFNU339vzzVM2VDkMEdeg7D1dDHEqSVA73VXFOXIAW8vlqzvpzhgIK9RGAOB38Gvnvy+WMblpjAlGwN5DIX1m0K1/LpY7js8jN4Hn8RENldvODyaQe+DT04fkygqRAMomu3s7NTBwcG6H5eo0DWbf2jae3xOugm/2PiuAM6IyJ6IDKlqZ+l23uykxOLCEhQXDORERBHHQE6JVVpDXmk7UVgxkFNiZS1uD1ltJworBnJKLK4QRHHBQE6J1XvdFWhqKM6jNDVIJOvtKdkYyCnZStMoTKtQBDGQU2Jt3//s9GQgw0RWI9lcjJKNgZwSK07NxSjZGMgpscLahpjILQZySqywtiEmcotNsyix4tZcjJLLk0AuIv8bwHYAbap6xot9EtVDGNsQE7lVc2pFRBYCWAngN7WfDhERueVFjnwHgLVgBS4RUSBqCuQi0g1gWFWfdPCzq0VkUEQGT58+XcthiYioQMUcuYj8C4B5Jg99FsBnADjqwK+quwDsAnILS7g4RyIislExkKvqn5ptF5GrAVwG4EkRAYAFAA6JyNtU9aSnZ0lERJY8W+pNRP4TQKeTqhUROQ3gRU8O7I9WAEmsvknq8waS+9yT+ryBaD73S1W1rXRjIHXkZicSJiIyaLYuXtwl9XkDyX3uSX3eQLyeu2eBXFUXebUvIiJyjlP0iYgijoHc3K6gTyAgSX3eQHKfe1KfNxCj5+7ZzU4iIgoGR+RERBHHQE5EFHEM5CZEZLuI/LuIHBGRh0RkTtDnVC8i8gEReVpEsiISi9IsOyJyvYg8KyK/EpF1QZ9PvYjIt0TklIg8FfS51JOILBSRH4vIM/nP+Z1Bn5MXGMjNPQbgKlVdDOA/AKwP+Hzq6SkANwL4adAn4jcRaQDwVQDvBnAlgA+JyJXBnlXdfBvA9UGfRAAmAdylqn8A4O0APhGH95yB3ISq/lBVJ/N/fQK59gOJoKrPqGpSVh9+G4BfqeqvVXUcwD8BeF/A51QXqvpTAC8FfR71pqonVPVQ/r9fBfAMgMg3pGcgr+xjAB4N+iTIF+0AjhX8/Thi8I+anBGRRQA6APw84FOpWWKXerPr6qiqD+d/5rPIfRX7Tj3PzW9OnntCiMk21uMmgIjMBrAHwCdV9ZWgz6dWiQ3kVl0dDSLyEQDvAbBCY1ZsX+m5J8hxAAsL/r4AwEhA50J1IiJNyAXx76jq3qDPxwtMrZgQkesBfBpAt6qOBX0+5JuDAC4XkctEpBnAzQD6Az4n8pHkem5/E8Azqnpv0OfjFQZyc18B8CYAj4nIL0Tk60GfUL2IyPtF5DiAdwD4vojsD/qc/JK/of3XAPYjd9OrT1WfDvas6kNEvgvgZwCuEJHjIvKXQZ9TnSwD8GEA1+b/bf9CRP4s6JOqFafoExFFHEfkREQRx0BORBRxDORERBHHQE5EFHEM5EREEcdATkQUcQzkREQR9/8BJvGUVgDbQoMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# my solver vs true beta\n",
    "plt.scatter(beta, beta_mat_lasso[:, -1])\n",
    "plt.scatter(beta, beta_mat_lasso[:, 50])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}