regression/regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![MOSEK ApS](https://www.mosek.com/static/images/branding/webgraphmoseklogocolor.png )\n",
    "\n",
    "# Linear Regression techniques using MOSEK Fusion API\n",
    "\n",
    "\n",
    "Regression is one of the most used techniques in finance, statistic, biology and in general any application where it is necessary to construct models that approximate input data.\n",
    "\n",
    "The aim of this tutorial is two-fold:\n",
    "\n",
    "1. to show some modeling techniques to define regression problems;\n",
    "2. to show how to efficiently implement those problems using MOSEK Fusion API.\n",
    "\n",
    "\n",
    "This tutorial is largerly based on:\n",
    "\n",
    "    [1] Schmelzer, T., Hauser, R., Andersen, E., & Dahl, J. (2013). Regression techniques for Portfolio Optimisation using MOSEK. arXiv preprint arXiv:1310.3397.\n",
    "    [2] Boyd, S., & Vandenberghe, L. (2004). Convex optimization. Cambridge university press.\n",
    "\n",
    "\n",
    "The code is written in Python3 and uses the [MOSEK Fusion API](https://mosek.com/documentation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MOSEK version (9, 0, 40)\n"
     ]
    }
   ],
   "source": [
    "import mosek\n",
    "from mosek.fusion import *\n",
    "import mosek.fusion.pythonic\n",
    "print(\"MOSEK version\", Model.getVersion())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The basics of linear regression\n",
    "\n",
    "A linear model is used to explain the relation between $m$ input variables $x_1,\\ldots,x_m$ and a dependent variable $y$ (the output). To this extent we assume that there exists a set of weights $w\\in \\mathbb{R}^{m}$ such that ideally\n",
    "\n",
    "$$\n",
    "y = x^T w,\n",
    "$$\n",
    "\n",
    "where $x=(x_1,\\ldots,x_m)$.\n",
    "\n",
    "In practice, we are given $n$ observations of the relation that links $x$ and $y$. We store them row-wise in a matrix $X\\in \\mathbb{R}^{n\\times m}$. Corresponding outputs $y=(y_i,\\ldots,y_n)$ are known as well. If the relation we want to describe is truly linear we can hope to solve the linear system\n",
    "\n",
    "$$\n",
    "X w =y.\n",
    "$$\n",
    "\n",
    "But usually this is not the case and we can only settle for an approximation which minimizes (in some sense) the residual vector\n",
    "\n",
    "$$\n",
    "r = Xw - y,\n",
    "$$\n",
    "\n",
    "and we obtain an optimization problem\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    " \\mbox{minimize} & \\phi(r)\\\\ \n",
    "  s.t. & r = Xw - y, \\\\\n",
    "       & w\\in W.\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "The choice of the cost function $\\phi(\\cdot)$ determines the way errors on the residuals are weighted. Here we demonstrate the following variants:\n",
    "\n",
    "* $\\ell_2$-norm,\n",
    "* $\\ell_1$-norm,\n",
    "* $\\ell_p-$norm ,\n",
    "* *Chebyshev* norm (also known as *Min-max*),\n",
    "* *Deadzone-linear*,\n",
    "* *$k$-largest residuals*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create simple input data sets which attempt to lead to the solution $w=(1,\\ldots,1)$ and we add small normally distributed error as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sys\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def dataset(m, n):\n",
    "    X = np.random.uniform(-1.0, 1.0, (n,m))\n",
    "    y = X.dot(np.ones(m)) + np.random.normal(0.0, 0.2, n)\n",
    "    return X, y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\ell_2$-norm\n",
    "\n",
    "The regression problem with respect to the $\\ell_2$-norm (i.e. linear least squares) has the conic quadratic formulation\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize} & t\\\\\n",
    "\\mbox{s.t.} & t\\geq \\|Xw-y\\|_2 \\quad \\mathrm{i.e.}\\ (t,Xw-y)\\in \\mathcal{Q}^{n+1},\\\\\n",
    " & t\\in\\mathbb{R}, w\\in\\mathbb{R}^m.\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "This model models the maximum likelihood of the fit precisely under the assumption that residuals have normal distribution.\n",
    "\n",
    "Below is the model in Fusion.            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def l2norm(X,y):\n",
    "    n,m=X.shape\n",
    "    \n",
    "    M = Model(\"l2norm\")\n",
    "    w = M.variable(m)\n",
    "    t = M.variable()\n",
    "    \n",
    "    M.objective(ObjectiveSense.Minimize, t)\n",
    "    \n",
    "    res = X @ w - y\n",
    "    M.constraint(Expr.vstack(t, res), Domain.inQCone())\n",
    "    \n",
    "    #M.setLogHandler(sys.stdout)\n",
    "    M.solve()\n",
    "    \n",
    "    # Return the weights vector and the residuals\n",
    "    w = w.level()\n",
    "    return w, X.dot(w)-y\n",
    "\n",
    "X,y = dataset(10,1000)\n",
    "w, res = l2norm(X,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If everything went well, we should have recovered the weight vector of all ones and normally distributed residuals, as below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|w-1|_2: 0.025785996710141198\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAF71JREFUeJzt3X20ZXV93/H3R3B8AgRkwAmIo3VQ0UbQG8Ra6wOShaIyTdBg0Q6W1VlKYmMS206euky07ZA0GtfSmlA1jtYHkChQzIM4QtQsUAZRFFBHcYCRcWYURhjxAeTbP/a+erice8++z/fueb/Wuuvsh985+3P2Pfd7f+e3z94nVYUkafl70GIHkCTNDQu6JPWEBV2SesKCLkk9YUGXpJ6woEtST1jQl4kk1yd53mLnWExJ/m2SW5PsTXL8PG7nOUm+PsX69yV5yxxsZ3WSSrL/DO57bJItA/PbkrxwtpkWWpKXJfnIYufoCwv6EjDsjzHJWUk+Nz5fVU+pqitGPM6MC8Qy8b+A36qqA6rq2vnaSFV9tqqeOF+PP0feTLM/HiDJuiTXJLkzyfYkf7ZUXxNVdQnw1CS/vNhZ+sCCrs6WQFF4LHB9l4ZLIOu8SbIKeD5w0SRNHg68ATgMeCZwEvDGBcg1033+YWD9XGbZV1nQl4nBXnySE5JsaXtgO5O8tW32mfZ2Tzss8awkD0ryR0luTrIryfuTPHLgcf99u+77Sf54wnbelOTCJP83yZ3AWe22r0yyJ8mOJO9IsmLg8SrJOUm2JrkryZuT/Iv2PncmuWCw/YTnODRrkock2QvsB3w5ybcmuX8l+c0kW4Gt7bInJbksye1Jvp7kFQPtX5zkhjbnd5K8sV3+vCTbB9odn+SLbbvzgYcOrLvfO6mBHE9op09Ncm373G9N8qYpfsdnJbmp3c63k5w5SdOTgS9W1Y+Hrayqd7XvMn5aVd8BPgg8e4rtVpLXtr+zO5K8M0nadZO+fgbeEZ6d5Bbg0wPLXtM+3zvax/6VJNe1r5t3TIhwBXDqZPk0DVXlzyL/ANuAF05YdhbwuWFtgCuBV7fTBwAnttOrgQL2H7jffwC+CTy+bfsx4APtumOBvcC/BlbQvIW/Z2A7b2rn19L8838Y8AzgRGD/dns3Am8Y2F4BlwAHAU8BfgJsbrf/SOAGYN0k+2HSrAOP/YQp9mMBlwGHtlkfAdwKvKbN+3Tge8BT2vY7gOe004cAT2+nnwdsb6dXADcDvwM8GDi93SdvGfZ7mpizfax/2e6/XwZ2Amsn/r7arHcCT2zXrRrPOeR5/jnwzlGvoYF1FwEbR+y3S4GDgaOB3cApHV4/4/nf3+Z/2MCyv6L5x/erwI/bDIcDRwK7gOcObP/Q9j4HLfbf4nL/sYe+dFzU9l72JNkD/O8p2t4DPCHJYVW1t6qumqLtmcBbq+qmqtoL/D5wRvv2+HTg/1XV56rqp8B/o/nDGnRlVV1UVfdV1Y+q6pqquqqq7q2qbcBfA8+dcJ9zq+rOqroe+CrwyXb7PwD+HpjsgOZUWbv6n1V1e1X9CHgJsK2q/qbN+0Xgb9vnDc1+PDbJQVV1R7t+ohNpCvlfVtU9VXUhcHXXMFV1RVV9pd1/19EML0zcX+PuoxlPflhV7Wj33zAHA3d12X6S1wBjTDLePmBjVe2pqluAy4Hj2uVdfidvqqoftvt83Jur6sdV9Ungh8CHq2pXNe8YPsv9XwPjz+XgLs9Jk7OgLx1rq+rg8R/gnCnang0cA3wtydVJXjJF21+i6WGOu5mmR3hEu+7W8RVVdTfw/Qn3v3VwJskxSS5N8t12GOZ/0IzVDto5MP2jIfMHzCBrV4N5Hws8c8I/yjOBR7frfx14MXBzkn9K8qxJMn2nqgb/0d08pN1QSZ6Z5PIku5P8AHgtD9xfVNUPgd9o1+9I8okkT5rkYe8ADuyw7bXARuBFVfW9Ec2/OzB9N7/4HXX5ndzvNdKazmtg/LnsGZFRI1jQl6Gq2lpVr6R5C3sucGGSR/DA3jXAbTSFbdzRwL00f2A7gKPGVyR5GPCoiZubMP8u4GvAmqo6CPgDIDN/Np2zdjWY91bgnwb/UVbzCZnXAVTV1VV1Gs1+vAi4YMjj7QCOHB9THsg17oc0ByEBSPJo7u9DNENQj6mqR9IMRQzdX1X1j1V1Ms1wy9eA/zPJc7yO5h/6pJKc0t7/pVX1lanajtDldzLbS7Y+mead1J2zfJx9ngV9GUryqiQrq+o+ftGr+RnN2Od9NOOd4z4M/E6SxyU5gKZHfX5V3QtcCLw0yb9qD1T+CaOL84E0Y7172x7k6+bsiU2ddSYuBY5J8uokD25/fiXJk5OsSHJmkkdW1T00z+lnQx7jSpoC9p+S7J/k14ATBtZ/GXhKkuOSPJTmuMOgA4Hbq+rHSU4A/t2woEmOSPOZ7EfQHHfYO0keaI4TPL3d3rDHegHNgdBfr6ovTPIYXc3172SY59IMxWmWLOjL0ynA9Wk++fF24Ix2vPJu4L8D/9wOMZwIvBf4AM0nYL5Nc4Dq9QDtGO3rgY/Q9ETvojlg9ZMptv1GmqJ0F00P8Pw5fF6TZp2JqrqL5qDcGTQ9ze/SvKN5SNvk1cC2dujotcCrhjzGT4Ffozn4eQfNsMjHBtZ/A/hT4FM0n6z53ISHOAf40yR30RyjGPYuAJq/xd9rc95OU+SGDrtV1U7g08BpkzzWH9McgP67NJ922ptkpgVzTn8nk3glzbEYzVLuPzSofVnbA9tDM5zy7cXOo8klORbYBJxQy/iPOMlLaT6x9YqRjTWSBX0f1/5BbaYZavkLmhNRnr6ci4S0r3LIRafRvM2/DVhDM3xjMZeWIXvoktQT9tAlqScW9AJGhx12WK1evXohNylJy94111zzvapaOardghb01atXs2XLltENJUk/l6TT2ckOuUhST1jQJaknLOiS1BMWdEnqCQu6JPWEBV2SesKCLkk9YUGXpJ6woEtSTyzomaLSYlq94ROTrtu28dQFTCLND3voktQTFnRJ6gkLuiT1hAVdknrCgi5JPWFBl6SesKBLUk9Y0CWpJyzoktQTIwt6kicm+dLAz51J3pDk0CSXJdna3h6yEIElScONLOhV9fWqOq6qjgOeAdwNfBzYAGyuqjXA5nZekrRIpjvkchLwraq6GTgN2NQu3wSsnctgkqTpmW5BPwP4cDt9RFXtAGhvD5/LYJKk6elc0JOsAF4GfHQ6G0iyPsmWJFt279493XySpI6m00N/EfDFqtrZzu9Msgqgvd017E5VdV5VjVXV2MqVK2eXVpI0qekU9Ffyi+EWgEuAde30OuDiuQolSZq+TgU9ycOBk4GPDSzeCJycZGu7buPcx5MkddXpG4uq6m7gUROWfZ/mUy+SpCXAM0UlqScs6JLUExZ0SeoJC7ok9YQFXZJ6woIuST1hQZeknrCgS1JPdDqxSFoKVm/4xJTrt208dYGSSEuTPXRJ6gkLuiT1hEMu6o1RQzKzua/DOVoO7KFLUk9Y0CWpJyzoktQTFnRJ6gkLuiT1hAVdknrCgi5JPdH1S6IPTnJhkq8luTHJs5IcmuSyJFvb20PmO6wkaXJde+hvB/6hqp4EPA24EdgAbK6qNcDmdl6StEhGFvQkBwH/BngPQFX9tKr2AKcBm9pmm4C18xVSkjRalx7644HdwN8kuTbJu5M8AjiiqnYAtLeHD7tzkvVJtiTZsnv37jkLLkm6vy4FfX/g6cC7qup44IdMY3ilqs6rqrGqGlu5cuUMY0qSRulS0LcD26vq8+38hTQFfmeSVQDt7a75iShJ6mJkQa+q7wK3Jnliu+gk4AbgEmBdu2wdcPG8JJQkddL18rmvBz6YZAVwE/Aamn8GFyQ5G7gFePn8RJQkddGpoFfVl4CxIatOmts4kqSZ8kxRSeoJC7ok9YRfQSd1MNVX1Pn1dFoq7KFLUk9Y0CWpJyzoktQTFnRJ6gkLuiT1hAVdknrCgi5JPWFBl6SesKBLUk9Y0CWpJzz1X0vKVKfY95WXFdBcsYcuST1hQZeknnDIRQtqXxxSkRaKPXRJ6gkLuiT1RKchlyTbgLuAnwH3VtVYkkOB84HVwDbgFVV1x/zE1HLisIq0OKbTQ39+VR1XVeNfFr0B2FxVa4DN7bwkaZHMZsjlNGBTO70JWDv7OJKkmepa0Av4ZJJrkqxvlx1RVTsA2tvDh90xyfokW5Js2b179+wTS5KG6vqxxWdX1W1JDgcuS/K1rhuoqvOA8wDGxsZqBhklSR106qFX1W3t7S7g48AJwM4kqwDa213zFVKSNNrIgp7kEUkOHJ8GfhX4KnAJsK5ttg64eL5CSpJG6zLkcgTw8STj7T9UVf+Q5GrggiRnA7cAL5+/mJKkUUYW9Kq6CXjakOXfB06aj1CSpOnzTFFJ6gkLuiT1hAVdknrCgi5JPWFBl6Se8AsuNG1eTXF63F9aKPbQJaknLOiS1BMOuUiz5JCKlgp76JLUExZ0SeoJC7ok9YQFXZJ6woIuST1hQZeknrCgS1JPWNAlqScs6JLUE50LepL9klyb5NJ2/nFJPp9ka5Lzk6yYv5iSpFGm00P/beDGgflzgbdV1RrgDuDsuQwmSZqeTgU9yVHAqcC72/kALwAubJtsAtbOR0BJUjdde+h/CfwX4L52/lHAnqq6t53fDhw57I5J1ifZkmTL7t27ZxVWkjS5kQU9yUuAXVV1zeDiIU1r2P2r6ryqGquqsZUrV84wpiRplC6Xz3028LIkLwYeChxE02M/OMn+bS/9KOC2+YspSRplZA+9qn6/qo6qqtXAGcCnq+pM4HLg9LbZOuDieUspSRppNp9D/6/A7yb5Js2Y+nvmJpIkaSam9Y1FVXUFcEU7fRNwwtxHkiTNhGeKSlJPWNAlqScs6JLUExZ0SeoJC7ok9YQFXZJ6woIuST1hQZeknrCgS1JPWNAlqScs6JLUE9O6loukpWX1hk9Mum7bxlMXMImWAnvoktQTFnRJ6gkLuiT1hAVdknrCgi5JPWFBl6SesKBLUk+MLOhJHprkC0m+nOT6JH/SLn9cks8n2Zrk/CQr5j+uJGkyXXroPwFeUFVPA44DTklyInAu8LaqWgPcAZw9fzElSaOMLOjV2NvOPrj9KeAFwIXt8k3A2nlJKEnqpNOp/0n2A64BngC8E/gWsKeq7m2bbAeOnOS+64H1AEcfffRs86qjqU4JB08LXy5G/R5nc19fA/3T6aBoVf2sqo4DjgJOAJ48rNkk9z2vqsaqamzlypUzTypJmtK0PuVSVXuAK4ATgYOTjPfwjwJum9tokqTp6PIpl5VJDm6nHwa8ELgRuBw4vW22Drh4vkJKkkbrMoa+CtjUjqM/CLigqi5NcgPwkSRvAa4F3jOPOSVJI4ws6FV1HXD8kOU30YynS5KWAM8UlaSesKBLUk/4FXT7KL+6TOofe+iS1BMWdEnqCQu6JPWEBV2SesKCLkk94adc9ACzucKfpMVjD12SesKCLkk94ZDLMubQiKRB9tAlqSfsoUv7KL+irn/soUtST1jQJaknLOiS1BMWdEnqiS5fEv2YJJcnuTHJ9Ul+u11+aJLLkmxtbw+Z/7iSpMl06aHfC/xeVT0ZOBH4zSTHAhuAzVW1BtjczkuSFsnIgl5VO6rqi+30XcCNwJHAacCmttkmYO18hZQkjTatMfQkq4Hjgc8DR1TVDmiKPnD4XIeTJHXXuaAnOQD4W+ANVXXnNO63PsmWJFt27949k4ySpA46FfQkD6Yp5h+sqo+1i3cmWdWuXwXsGnbfqjqvqsaqamzlypVzkVmSNESXT7kEeA9wY1W9dWDVJcC6dnodcPHcx5MkddXlWi7PBl4NfCXJl9plfwBsBC5IcjZwC/Dy+YkoSepiZEGvqs8BmWT1SXMbR5I0U54pKkk9YUGXpJ6woEtST1jQJaknLOiS1BMWdEnqCQu6JPWEBV2SeqLLmaJaJKO+lV2SBtlDl6SesKBLUk845CJpqNkM+W3beOocJlFX9tAlqScs6JLUEw65LCI/xSJpLtlDl6SesKBLUk845CJpzo0aTvRTMPPDHrok9cTIgp7kvUl2JfnqwLJDk1yWZGt7e8j8xpQkjdKlh/4+4JQJyzYAm6tqDbC5nZckLaKRBb2qPgPcPmHxacCmdnoTsHaOc0mSpmmmY+hHVNUOgPb28MkaJlmfZEuSLbt3757h5iRJo8z7QdGqOq+qxqpqbOXKlfO9OUnaZ820oO9Msgqgvd01d5EkSTMx04J+CbCunV4HXDw3cSRJM9XlY4sfBq4Enphke5KzgY3AyUm2Aie385KkRTTyTNGqeuUkq06a4yySpFnw1P955hUVpenxsgEz56n/ktQTFnRJ6gkLuiT1hAVdknrCg6KSemNfP6BqD12SesKCLkk94ZCLpAXn+Rnzwx66JPWEBV2SesKCLkk9YUGXpJ6woEtST/gpF0nLip+QmZw9dEnqiWXTQ1+sU3r39VOJJS0f9tAlqScs6JLUE7MacklyCvB2YD/g3VW1JL8sej4PoniARhJMXQsWamh2xj30JPsB7wReBBwLvDLJsXMVTJI0PbMZcjkB+GZV3VRVPwU+Apw2N7EkSdM1myGXI4FbB+a3A8+c2CjJemB9O7s3yddnsc2JDgO+B5Bz5/BRF8bPsy9Tyzn/cs4Oyzv/omafZZ2YcfY5qE+P7dJoNgU9Q5bVAxZUnQecN4vtTB4g2VJVY/Px2PNtOWeH5Z1/OWeH5Z3f7PNrNkMu24HHDMwfBdw2uziSpJmaTUG/GliT5HFJVgBnAJfMTSxJ0nTNeMilqu5N8lvAP9J8bPG9VXX9nCXrZl6GchbIcs4Oyzv/cs4Oyzu/2edRqh4w7C1JWoY8U1SSesKCLkk9sawKepJDk1yWZGt7e8gk7Y5O8skkNya5IcnqhU06NFOn7G3bg5J8J8k7FjLjVLrkT3JckiuTXJ/kuiS/sRhZB/KckuTrSb6ZZMOQ9Q9Jcn67/vNL4XUyrkP2321f29cl2Zyk0+eUF8qo/APtTk9SSZbMxwG7ZE/yinb/X5/kQwudcVJVtWx+gD8DNrTTG4BzJ2l3BXByO30A8PDlkr1d/3bgQ8A7Fjv3dPIDxwBr2ulfAnYABy9S3v2AbwGPB1YAXwaOndDmHOCv2ukzgPMXez9PI/vzx1/XwOuWSvau+dt2BwKfAa4CxhY79zT2/RrgWuCQdv7wxc49/rOseug0lxbY1E5vAtZObNBeT2b/qroMoKr2VtXdCxdxUiOzAyR5BnAE8MkFytXVyPxV9Y2q2tpO3wbsAlYuWML763JpisHndCFwUpJhJ8wttJHZq+rygdf1VTTngSwVXS8L8maajsKPFzLcCF2y/0fgnVV1B0BV7VrgjJNabgX9iKraAdDeHj6kzTHAniQfS3Jtkj9vLyS22EZmT/Ig4C+A/7zA2brosu9/LskJND2cby1AtmGGXZriyMnaVNW9wA+ARy1Iuql1yT7obODv5zXR9IzMn+R44DFVdelCBuugy74/BjgmyT8nuaq96uySsOS+sSjJp4BHD1n1hx0fYn/gOcDxwC3A+cBZwHvmIt9U5iD7OcDfVdWti9FRnIP844+zCvgAsK6q7puLbDPQ5dIUnS5fsQg650ryKmAMeO68JpqeKfO3HZe30fxdLjVd9v3+NMMuz6N5Z/TZJE+tqj3znG2kJVfQq+qFk61LsjPJqqra0RaNYW91tgPXVtVN7X0uAk5kAQr6HGR/FvCcJOfQjP2vSLK3qiY9qDSX5iA/SQ4CPgH8UVVdNU9Ru+hyaYrxNtuT7A88Erh9YeJNqdNlNZK8kOaf7XOr6icLlK2LUfkPBJ4KXNF2XB4NXJLkZVW1ZcFSDtf1dXNVVd0DfLu94OAamrPnF9VyG3K5BFjXTq8DLh7S5mrgkCTjY7cvAG5YgGyjjMxeVWdW1dFVtRp4I/D+hSrmHYzM314C4uM0uT+6gNmG6XJpisHndDrw6WqPci2ykdnbIYu/Bl62lMZwW1Pmr6ofVNVhVbW6fa1fRfM8FruYQ7fXzUU0B6VJchjNEMxNC5pyMot9VHY6PzTjm5uBre3toe3yMZpvTBpvdzJwHfAV4H3AiuWSfaD9WSytT7mMzA+8CrgH+NLAz3GLmPnFwDdoxvH/sF32pzTFA+ChwEeBbwJfAB6/2Pt5Gtk/Bewc2M+XLHbm6eSf0PYKlsinXDru+wBvpekofgU4Y7Ezj/946r8k9cRyG3KRJE3Cgi5JPWFBl6SesKBLUk9Y0CWpJyzoktQTFnRJ6on/DzEGHLd0X6TOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f8235d157b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"|w-1|_2:\", np.linalg.norm(w-np.ones(10)))\n",
    "plt.hist(res,40)\n",
    "plt.title(\"Histogram of residuals (l2 norm)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\ell_1$-norm\n",
    "\n",
    "Under the $\\ell_1$-norm we minimize $\\sum_i |r_i|$. The corresponding linear formulation is\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize}  & \\sum_i t_i \\\\\n",
    "\\mbox{s.t.} & t\\geq Xw-y,\\\\\n",
    " & t \\geq -(Xw-y), \\\\\n",
    " & w\\in\\mathbb{R}^m, t\\in\\mathbb{R}^n,\n",
    "\\end{array}\n",
    "$$ \n",
    "\n",
    "which simply expresses the bound $t_i\\geq |r_i|$. Below is a corresponding Fusion model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|w-1|_2: 0.034229179355081116\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAGIhJREFUeJzt3XuU3GV9x/H3R0LkThKyiRHEYAlotAq4xqC1IBcPGiWpRYUiDTanOYpary3x0mrVtqEq6jlSNRV1tQLBiCTFGzESFQ8gy91wMRACicRkEQKEe+TbP37PlmEzs/Ob3bnsPPm8ztkzv8szM5/57ex3n3l+l1FEYGZm3e9ZnQ5gZmbN4YJuZpYJF3Qzs0y4oJuZZcIF3cwsEy7oZmaZcEHvEpLWSDq60zk6SdJfSdogaZukw1v4PK+RdNsw678l6TNNeJ7pkkLSuBHcd6ak/or59ZKOG22mdpN0oqQLOp0jFy7oY0C1P0ZJp0u6fHA+Il4cEavrPM6IC0SX+BzwnojYKyKua9WTRMSvIuLQVj1+k3yaYnvsQNJLJP1U0r2SxvSJJhGxAniJpJd2OksOXNCttDHwj+L5wJoyDcdA1paRNA14LXBxjSZPAhcCC9oWilFt8/OBhc3MsrNyQe8Slb14SbMk9Ut6UNJmSWenZr9Mt1vTsMSRkp4l6eOS7pK0RdK3Je1b8bh/m9b9UdI/D3meT0paJul/JD0InJ6e+wpJWyVtkvRlSeMrHi8knSFpraSHJH1a0p+l+zwo6cLK9kNeY9Wskp4taRuwC3CDpDtq3D8kvVvSWmBtWvZCSSsl3SfpNklvrWj/Bkk3p5y/l/ThtPxoSRsr2h0u6drUbimwW8W6Z3ySqshxcJqeI+m69No3SPrkML/j0yWtS89zp6RTazQ9Hrg2Ih6rtjIibouIcyn/zy8kvTP9zu6XdI4kpXU13z8VnwgXSLob+HnFsnek13t/euxXSLoxvW++PCTCamBOmaxWR0T4p8M/wHrguCHLTgcur9YGuAI4LU3vBcxO09OBAMZV3O/vgNuBF6S2FwHfSetmAtuAvwDGU3yEf7LieT6Z5udR/PPfHXg5MBsYl57vFuD9Fc8XwApgH+DFwOPAqvT8+wI3A/NrbIeaWSse++BhtmMAK4FJKeuewAbgHSnvEcC9wItT+03Aa9L0ROCINH00sDFNjwfuAj4A7AqclLbJZ6r9nobmTI/152n7vRTYDMwb+vtKWR8EDk3rpg3mrPI6PwucU+I9dDAQJd5/AVwCTAAOBAaAE0q8fwbzfzvl371i2Vcp/vG9DniM4tPEFGB/YAtwVMXzT0r32afTf4vd/uMe+thxceq9bJW0FfivYdo+CRwsaXJEbIuIK4dpeypwdkSsi4htwEeAk9PH45OA/42IyyPiCeBfKP6wKl0RERdHxFMR8WhEXBMRV0bE9ohYD3wNOGrIfc6KiAcjYg3wW+DS9PwPAD8Gau3QHC5rWf8REfdFxKPAG4H1EfHNlPda4PvpdUOxHWdK2ici7k/rh5pNUci/GBFPRsQy4OqyYSJidUTclLbfjRTDC0O316CnKMaTd4+ITWn7VTMBeKhshpIWR8TWiLgbuAw4LC0v8zv5ZEQ8nLb5oE9HxGMRcSnwMHB+RGyJiN8Dv+KZ74HB1zKhya9pp+OCPnbMi4gJgz/AGcO0XQAcAtwq6WpJbxym7XMpepiD7qLoEU5N6zYMroiIR4A/Drn/hsoZSYdIukTSH9IwzL8Dk4fcZ3PF9KNV5vcaQdayKvM+H3jlkH+UpwLPSev/GngDcJekX0g6skam30dE5T+6u6q0q0rSKyVdJmlA0gPAO9lxexERDwNvS+s3SfqhpBfWeNj7gb3LZijpDxXTj/D076jM7+QZ75GkkffA4GvZ2kBeq8IFvQtFxNqIOIXiI+xZwDJJe7Jj7xrgHorCNuhAYDvFH9gm4IDBFZJ2B/Yb+nRD5r8C3ArMiIh9gI8CGvmrKZ21rMq8G4BfVP6jjOIImXcBRMTVETGXYjteTLEjcahNwP6DY8oVuQY9DOwxOCPpOTzTeRRDUM+LiH0phiKqbq+I+GlEHE8x3HIr8N81XuONFP/Q26HM72S0R9K8iOKT1IOjfJydngt6F5L0dkk9EfEUT/dq/kQx9vkUxXjnoPOBD0g6SNJeFD3qpRGxHVgGvEnSq9KOyn+lfnHem2Ksd1vqQb6raS9s+KwjcQlwiKTTJO2afl4h6UWSxks6VdK+EfEkxWv6U5XHuIKigP2DpHGS3gzMqlh/A/BiSYdJ2o1iv0OlvYH7IuIxSbOAv6kWVNJUFcdk70mx32FbjTxQ7Cc4Ij1ftcdSWjc+ze8m6dk1HqueZv9OqjmKYijORskFvTudAKxRceTHl4CT03jlI8C/Ab9OQwyzgW8A36E4AuZOih1U7wVIY7TvBS6g6Ik+RLHD6vFhnvvDFEXpIYoe5NImvq6aWUciIh6i2Cl3MkVP8w8Un2gGi9tpwPo0dPRO4O1VHuMJ4M0UOz/vpxgWuahi/e+ATwE/oziy5vIhD3EG8ClJD1Hso6j2KQCKv8UPpZz3URS5qsNuEbEZ+Dkwt8ZjPZ9iWGNwDP5RoOaJUnU09XdSwykU+2JslPTMoUHbmaUe2FaK4ZQ7O53HapM0E+gDZkUX/xFLehPFEVtvrdvY6nJB38mlP6hVFEMtnwdeSXH4nt8YZl3GQy42l+Jj/j3ADIrhGxdzsy7kHrqZWSZK9dAlfUDF1f5+K+n8tNf8IElXpdOFl6rG6dxmZtYedXvokvan2HM/MyIelXQh8COKEzIuiogLJH0VuCEivjLcY02ePDmmT5/enORmZjuJa6655t6I6KnXruwp1eOA3SU9SXESxSbgGJ4+praP4vjbYQv69OnT6e/vH66JmZkNIanU2cl1h1zStRc+B9xNUcgfAK4BtlacXLCR4qI7ZmbWIXULuqSJFEdCHERxXYc9gddXaVp17EbSQhWXeu0fGBgYTVYzMxtGmZ2ixwF3RsRAOkX6IuBVwISKK64dQHHY2w4iYklE9EZEb09P3SEgMzMboTIF/W5gtqQ90gWKjqW4pvVlPH0Z0vnA8tZENDOzMsqMoV9FcRGna4Gb0n2WAGcCH5R0O8UV+s5tYU4zM6uj1FEuEfEJ4BNDFq/jmVedMzOzDvKp/2ZmmXBBNzPLhAu6mVkmGvnyXbOuNn3RD2uuW794ThuTmLWGe+hmZplwQTczy4QLuplZJlzQzcwy4YJuZpYJF3Qzs0y4oJuZZcIF3cwsEy7oZmaZcEE3M8uEC7qZWSZc0M3MMuGCbmaWCRd0M7NM1C3okg6VdH3Fz4OS3i9pkqSVktam24ntCGxmZtWV+ZLo2yLisIg4DHg58AjwA2ARsCoiZgCr0ryZmXVIo0MuxwJ3RMRdwFygLy3vA+Y1M5iZmTWm0YJ+MnB+mp4aEZsA0u2UZgYzM7PGlC7oksYDJwLfa+QJJC2U1C+pf2BgoNF8ZmZWUiM99NcD10bE5jS/WdI0gHS7pdqdImJJRPRGRG9PT8/o0pqZWU2NFPRTeHq4BWAFMD9NzweWNyuUmZk1rlRBl7QHcDxwUcXixcDxktamdYubH8/MzMoaV6ZRRDwC7Ddk2R8pjnoxM7MxwGeKmpllwgXdzCwTLuhmZplwQTczy4QLuplZJkod5WKWu+mLfjjs+vWL57QpidnIuYduZpYJF3Qzs0y4oJuZZcIF3cwsEy7oZmaZcEE3M8uEC7qZWSZc0M3MMuGCbmaWCRd0M7NMuKCbmWXCBd3MLBNlv1N0gqRlkm6VdIukIyVNkrRS0tp0O7HVYc3MrLayV1v8EvCTiDhJ0nhgD+CjwKqIWCxpEbAIOLNFOc3qXhHRbGdXt4cuaR/gL4FzASLiiYjYCswF+lKzPmBeq0KamVl9ZYZcXgAMAN+UdJ2kr0vaE5gaEZsA0u2UaneWtFBSv6T+gYGBpgU3M7NnKlPQxwFHAF+JiMOBhymGV0qJiCUR0RsRvT09PSOMaWZm9ZQp6BuBjRFxVZpfRlHgN0uaBpBut7QmopmZlVF3p2hE/EHSBkmHRsRtwLHAzelnPrA43S5vaVLbKXjHp9nIlT3K5b3Ad9MRLuuAd1D07i+UtAC4G3hLayKamVkZpQp6RFwP9FZZdWxz45iZ2Uj5TFEzs0y4oJuZZcIF3cwsEy7oZmaZcEE3M8uEC7qZWSZc0M3MMuGCbmaWibJnippZiwx3uYP1i+e0MYl1O/fQzcwy4YJuZpYJF3Qzs0y4oJuZZcI7Rc1K8I5L6wbuoZuZZcIF3cwsEy7oZmaZcEE3M8tEqZ2iktYDDwF/ArZHRK+kScBSYDqwHnhrRNzfmphmZlZPI0e5vDYi7q2YXwSsiojFkhal+TObms6yM9zRIt2q3mvyUTDWLqMZcpkL9KXpPmDe6OOYmdlIlS3oAVwq6RpJC9OyqRGxCSDdTql2R0kLJfVL6h8YGBh9YjMzq6rskMurI+IeSVOAlZJuLfsEEbEEWALQ29sbI8hoZmYllOqhR8Q96XYL8ANgFrBZ0jSAdLulVSHNzKy+ugVd0p6S9h6cBl4H/BZYAcxPzeYDy1sV0szM6isz5DIV+IGkwfbnRcRPJF0NXChpAXA38JbWxTQzs3rqFvSIWAe8rMryPwLHtiKUmZk1zmeKmpllwgXdzCwTLuhmZplwQTczy4QLuplZJlzQzcwy4YJuZpYJF3Qzs0w0cj10s1JyvOa5WTdwD93MLBPuoVvD3AM3G5vcQzczy4QLuplZJlzQzcwy4YJuZpYJF3Qzs0y4oJuZZcIF3cwsE6ULuqRdJF0n6ZI0f5CkqyStlbRU0vjWxTQzs3oaObHofcAtwD5p/izgCxFxgaSvAguArzQ5n1nX84lY1i6leuiSDgDmAF9P8wKOAZalJn3AvFYENDOzcsoOuXwR+CfgqTS/H7A1Iran+Y3A/tXuKGmhpH5J/QMDA6MKa2ZmtdUt6JLeCGyJiGsqF1dpGtXuHxFLIqI3Inp7enpGGNPMzOopM4b+auBESW8AdqMYQ/8iMEHSuNRLPwC4p3Uxzcysnro99Ij4SEQcEBHTgZOBn0fEqcBlwEmp2XxgectSmplZXaM5Dv1M4IOSbqcYUz+3OZHMzGwkGroeekSsBlan6XXArOZHMjOzkfAXXJh1seGOcV+/eE4bk9hY4FP/zcwy4YJuZpYJF3Qzs0y4oJuZZcIF3cwsEy7oZmaZcEE3M8uEC7qZWSZc0M3MMuGCbmaWCRd0M7NMuKCbmWXCBd3MLBMu6GZmmXBBNzPLhAu6mVkmXNDNzDJRt6BL2k3SbyTdIGmNpH9Nyw+SdJWktZKWShrf+rhmZlZLmR7648AxEfEy4DDgBEmzgbOAL0TEDOB+YEHrYpqZWT11C3oUtqXZXdNPAMcAy9LyPmBeSxKamVkppb4kWtIuwDXAwcA5wB3A1ojYnppsBPavcd+FwEKAAw88cLR5rQ2G++Jhay//LqwRpXaKRsSfIuIw4ABgFvCias1q3HdJRPRGRG9PT8/Ik5qZ2bAaOsolIrYCq4HZwARJgz38A4B7mhvNzMwaUeYolx5JE9L07sBxwC3AZcBJqdl8YHmrQpqZWX1lxtCnAX1pHP1ZwIURcYmkm4ELJH0GuA44t4U5zcysjroFPSJuBA6vsnwdxXi6mXWhejtc1y+e06Yk1iw+U9TMLBMu6GZmmXBBNzPLhAu6mVkmXNDNzDJR6tR/6z4+gsF82YCdj3voZmaZcEE3M8uEh1y62Gg+UvvjuFl+3EM3M8uEC7qZWSZc0M3MMuGCbmaWCRd0M7NMuKCbmWXCBd3MLBMu6GZmmXBBNzPLRJkviX6epMsk3SJpjaT3peWTJK2UtDbdTmx9XDMzq6VMD3078KGIeBEwG3i3pJnAImBVRMwAVqV5MzPrkLoFPSI2RcS1afoh4BZgf2Au0Jea9QHzWhXSzMzqa2gMXdJ04HDgKmBqRGyCougDU2rcZ6Gkfkn9AwMDo0trZmY1lS7okvYCvg+8PyIeLHu/iFgSEb0R0dvT0zOSjGZmVkKpgi5pV4pi/t2IuCgt3ixpWlo/DdjSmohmZlZGmaNcBJwL3BIRZ1esWgHMT9PzgeXNj2dmZmWV+YKLVwOnATdJuj4t+yiwGLhQ0gLgbuAtrYloZmZl1C3oEXE5oBqrj21uHDMzGymfKWpmlgkXdDOzTPhLos2squG+SHz94jltTGJluYduZpYJF3Qzs0y4oJuZZcIF3cwsEy7oZmaZcEE3M8uEC7qZWSZ8HPoYNtxxwGadNNr3po9jbw330M3MMuGCbmaWCRd0M7NMuKCbmWXCBd3MLBM+yqWDfBSLmTWTe+hmZpko8yXR35C0RdJvK5ZNkrRS0tp0O7G1Mc3MrJ4yPfRvAScMWbYIWBURM4BVad7MzDqobkGPiF8C9w1ZPBfoS9N9wLwm5zIzswaNdAx9akRsAki3U2o1lLRQUr+k/oGBgRE+nZmZ1dPynaIRsSQieiOit6enp9VPZ2a20xppQd8saRpAut3SvEhmZjYSIy3oK4D5aXo+sLw5cczMbKTKHLZ4PnAFcKikjZIWAIuB4yWtBY5P82Zm1kF1zxSNiFNqrDq2yVnMzOqeQe1rqdfmM0XNzDLhgm5mlgkXdDOzTLigm5llwgXdzCwTLuhmZplwQTczy4S/scjMsrGzH8PuHrqZWSZc0M3MMuEhlxbzF0Gb7Wg0fxf+m6rNPXQzs0y4oJuZZcJDLqPkj39mNla4h25mlomu6aF36vhS98DN8pH7ceruoZuZZcIF3cwsE6MacpF0AvAlYBfg6xExJr9b1MMmZtZqw9WZdg3ljLiHLmkX4Bzg9cBM4BRJM5sVzMzMGjOaIZdZwO0RsS4ingAuAOY2J5aZmTVqNEMu+wMbKuY3Aq8c2kjSQmBhmt0m6bZRPGc1k4F7dVaTH7V9JgP3djrECHVzduju/N2cHcZo/pJ1pOHsTahPzy/TaDQFXVWWxQ4LIpYAS0bxPMOHkPojordVj99q3Zy/m7NDd+fv5uzQ3fnHcvbRDLlsBJ5XMX8AcM/o4piZ2UiNpqBfDcyQdJCk8cDJwIrmxDIzs0aNeMglIrZLeg/wU4rDFr8REWualqy8lg3ntEk35+/m7NDd+bs5O3R3/jGbXRE7DHubmVkX8pmiZmaZcEE3M8tE1xV0SZMkrZS0Nt1OrNHuQEmXSrpF0s2Sprc3aXVl86e2+0j6vaQvtzNjLWWySzpM0hWS1ki6UdLbOpG1Is8Jkm6TdLukRVXWP1vS0rT+qrHyPhlUIv8H0/v7RkmrJJU6Xrld6uWvaHeSpJA0Zg4HLJNd0lvT9l8j6bx2Z9xBRHTVD/CfwKI0vQg4q0a71cDxaXovYI9OZ28kf1r/JeA84Mudzl02O3AIMCNNPxfYBEzoUN5dgDuAFwDjgRuAmUPanAF8NU2fDCzt9HZuMP9rB9/bwLu6LX9qtzfwS+BKoLfTuRvY9jOA64CJaX5Kp3N3XQ+d4vICfWm6D5g3tEG6psy4iFgJEBHbIuKR9kUcVt38AJJeDkwFLm1TrjLqZo+I30XE2jR9D7AF6Glbwmcqc3mKyte0DDhWUrWT5jqhbv6IuKzivX0lxfkgY0XZy4N8mqKz8Fg7w9VRJvvfA+dExP0AEbGlzRl30I0FfWpEbAJIt1OqtDkE2CrpIknXSfpsupjYWFA3v6RnAZ8H/rHN2eops+3/n6RZFL2bO9qQrZpql6fYv1abiNgOPADs15Z09ZXJX2kB8OOWJmpM3fySDgeeFxGXtDNYCWW2/SHAIZJ+LenKdPXZjhqT31gk6WfAc6qs+ljJhxgHvAY4HLgbWAqcDpzbjHz1NCH/GcCPImJDuzuLTcg++DjTgO8A8yPiqWZkG4Eyl6codQmLDimdTdLbgV7gqJYmasyw+VPH5QsUf5tjTZltP45i2OVoik9Gv5L0kojY2uJsNY3Jgh4Rx9VaJ2mzpGkRsSkVjWofczYC10XEunSfi4HZtKmgNyH/kcBrJJ1BMf4/XtK2iKi5U6lZmpAdSfsAPwQ+HhFXtihqGWUuTzHYZqOkccC+wH3tiVdXqctrSDqO4h/uURHxeJuylVEv/97AS4DVqePyHGCFpBMjor9tKasr+965MiKeBO5MFx6cQXEWfUd045DLCmB+mp4PLK/S5mpgoqTBsdtjgJvbkK2Muvkj4tSIODAipgMfBr7djmJeQt3s6TIQP6DI/L02ZqumzOUpKl/TScDPI+3hGgPq5k9DFl8DThwLY7hDDJs/Ih6IiMkRMT2916+keB2dLuZQ7r1zMcVOaSRNphiCWdfWlEN1eq9soz8U45urgLXpdlJa3kvxrUmD7Y4HbgRuAr4FjO909kbyV7Q/nbFzlEvd7MDbgSeB6yt+Dutg5jcAv6MYx/9YWvYpisIBsBvwPeB24DfACzq9nRvM/zNgc8W2XtHpzI3kH9J2NWPkKJeS217A2RSdxZuAkzud2af+m5llohuHXMzMrAoXdDOzTLigm5llwgXdzCwTLuhmZplwQTczy4QLuplZJv4PHYpdXs+2RQoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f81cab32cf8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def l1norm(X,y):\n",
    "    n,m=X.shape\n",
    "    \n",
    "    M = Model(\"l1norm\")\n",
    "    w = M.variable(m)\n",
    "    t = M.variable(n)\n",
    "    \n",
    "    M.objective(ObjectiveSense.Minimize, Expr.sum(t))\n",
    "    \n",
    "    res = X @ w - y\n",
    "    M.constraint(t >= res)\n",
    "    M.constraint(t >= -res)\n",
    "    \n",
    "    #M.setLogHandler(sys.stdout)\n",
    "    M.solve()\n",
    "    \n",
    "    # Return the weights vector and the residuals\n",
    "    w = w.level()\n",
    "    return w, X.dot(w)-y\n",
    "\n",
    "w, res = l1norm(X,y)\n",
    "print(\"|w-1|_2:\", np.linalg.norm(w-np.ones(10)))\n",
    "plt.hist(res,40)\n",
    "plt.title(\"Histogram of residuals (l1 norm)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\ell_p$-norm\n",
    "\n",
    "The $\\ell_p$ norm\n",
    "\n",
    "$$\\phi(r) = \\| r \\|_p = \\left( \\sum_i |r_i|^p \\right)^{1/p}$$\n",
    "\n",
    "For a general $p>1$ we can model it using the power cone. The standard representation of the epigraph $t\\geq\\|r\\|_p$ as a conic model can be found in https://docs.mosek.com/modeling-cookbook/powo.html#norm-cones\n",
    "\n",
    "$$\n",
    "\\begin{array}{l}\n",
    "\\sum_i s_i = t\\\\\n",
    "s_it^{p-1} \\geq |r_i|^p, \\ i=1,\\ldots,n\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "and the last constraint can be expressed using a power cone. In Fusion all of these individual conic constraints can be specified at once in matrix for, where, as always in Fusion, \"matrix in a cone\" is a shorthand for \"each row of the matrix is in the cone\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|w-1|_2: 0.025786431917413328\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAF79JREFUeJzt3XuUnHV9x/H3B0IEuYXIJkQurpSARqsB1hhqLSrgQVCSWlQo2mBzmqOo9UbbeGkPXtqGalHPkaqpqNEqBiOSCFWJkah4AFmIgCFgIAQSE5IVEkJALpFv/3h+a4dlZueZ3bns/PJ5nTNnnstv5vnMM7Pf/c1v5nlGEYGZmXW/PTodwMzMmsMF3cwsEy7oZmaZcEE3M8uEC7qZWSZc0M3MMuGC3iUkrZb0qk7n6CRJfylpg6Sdko5t4XZeKenOYdZ/TdInm7CdXkkhadwIbjtNUn/F/HpJJ482U7tJOkPStzudIxcu6GNAtT9GSedKunZwPiJeFBEr69zPiAtEl/g08O6I2C8iVrVqIxHx84g4plX33ySfoNgfXS0ilgEvlvSSTmfJgQu6lTYG/lE8D1hdpuEYyNoykqYArwau6HSWSqPY55cC85qZZXflgt4lKnvxkmZI6pe0Q9IWSRelZj9L19vTsMQJkvaQ9FFJ90raKunrkg6suN+/SesekPTPQ7ZzgaQlkv5H0g7g3LTt6yRtl7RZ0uclja+4v5B0nqS1kh6W9AlJf5Jus0PSZZXthzzGqlklPUvSTmBP4BZJd9e4fUh6l6S1wNq07AWSlkt6UNKdkt5c0f40SbennL+VdH5a/ipJGyvaHSvp5tRuMbB3xbqnvZOqyHFUmj5d0qr02DdIumCY5/hcSevSdu6RdE6NpqcAN0fEYzXuZ/B5W5zu62ZJLx1muyHpHek52ybpYklK62q+fireEc6VdB/wk4plb0+Pd1u675dJujW9bj4/JMJK4PRa+awBEeFLhy/AeuDkIcvOBa6t1ga4Dnhbmt4PmJmme4EAxlXc7m+Bu4AjU9vLgW+kddOAncCfA+Mp3sI/WbGdC9L8bIp//vsAxwMzgXFpe2uA91VsL4BlwAHAi4DHgRVp+wcCtwNzauyHmlkr7vuoYfZjAMuBiSnrvsAG4O0p73HA74AXpfabgVem6YOA49L0q4CNaXo8cC/wfmAv4My0Tz5Z7XkamjPd15+m/fcSYAswe+jzlbLuAI5J66YM5qzyOD8FXFzrNVTxvJ2ZMp8P3APsNcx+uxKYABwBDACnlnj9DOb/esq/T8WyL1L843st8BjFu4lJwKHAVuDEiu1PTLc5oNN/i91+cQ997Lgi9V62S9oO/NcwbZ8EjpJ0cETsjIjrh2l7DnBRRKyLiJ3Ah4Cz0tvjM4HvR8S1EfEE8C8Uf1iVrouIKyLiqYj4fUTcFBHXR8SuiFgPfAk4cchtLoyIHRGxGvg1cHXa/kPAD4BaH2gOl7Wsf4+IByPi98DrgfUR8dWU92bgu+lxQ7Efp0k6ICK2pfVDzaQoip+NiCcjYglwY9kwEbEyIm5L++9WiuGFoftr0FMU48n7RMTmtP+qmQA8XGfTN0XEkoh4EriIorjOHKb9gojYHhH3AdcA09PyMs/JBRHxSNrngz4REY9FxNXAI8ClEbE1In4L/JynvwYGH8uEOo/J6nBBHztmR8SEwQtw3jBt5wJHA3dIulHS64dp+1yKHuageyl6hJPTug2DKyLiUeCBIbffUDkj6WhJV0q6Pw3D/Btw8JDbbKmY/n2V+f1GkLWsyrzPA14+5B/lOcAhaf1fAacB90r6qaQTamT6bURU/qO7t0q7qiS9XNI1kgYkPQS8g2fuLyLiEeAtaf1mSVdJekGNu90G7F9n05XP61PAxvRYarm/YvpR/v85KvOcPO01kjTyGhh8LNuHyWcluKB3oYhYGxFnU7yFvRBYImlfntm7BthEUdgGHQHsovgD2wwcNrhC0j7Ac4Zubsj8F4A7gKkRcQDwYUAjfzSls5ZVmXcD8NPKf5RRfEPmnQARcWNEzKLYj1cAl1W5v83AoYNjyhW5Bj0CPHtwRtIhPN23KIagDo+IAymGIqrur4j4UUScQjHccgfw3zUe460U/9CHc3hFpj0onudNdW5TTZnnZLSnbH0hxTupHaO8n92eC3oXkvRWST2p5zXYq/kDxdjnUxTjnYMuBd4v6fmS9qPoUS+OiF3AEuANkv4sfVD5MeoX5/0pxnp3ph7kO5v2wIbPOhJXAkdLepukvdLlZZJeKGm8pHMkHZiGJXZQ7MOhrqMoYH8vaZykNwIzKtbfArxI0nRJe1OMX1faH3gwIh6TNAP462pBJU1W8Z3sfSk+d9hZIw8UnxMcl7ZXy/GS3piGRt6X7nO4oblamv2cVHMixVCcjZILenc6FVit4psfnwPOSuOVjwL/CvwiDTHMBL4CfIPiGzD3UHxA9R6ANEb7HuDbFD3Rhyk+sHp8mG2fT1GUHqboQS5u4uOqmXUkIuJhig/lzqLoad5P8Y7mWanJ24D1aejoHcBbq9zHE8AbKT783EYxLHJ5xfrfAB8HfkzxzZprh9zFecDHJT1M8RlFtXcBUPwtfjDlfJCiyFUddouILcBPgFm1HjuwNGXdlh7nG9M/rkY19Tmp4WyKz2JslPT0oUHbnaUe2HaK4ZR7Op3HapM0DVgEzBgyvk/6auRREfGMf1BjjaQ3UHxj6811G1td7qHv5iS9QdKz01v9TwO3UXwFzsawiLg9Il42tJh3m4j4vot587ig2yyKt/mbgKkUwzddXSTMdlcecjEzy4R76GZmmWjrCYwOPvjg6O3tbecmzcy63k033fS7iOip166tBb23t5f+/v76Dc3M7I8klTo62UMuZmaZcEE3M8uEC7qZWSZc0M3MMuGCbmaWCRd0M7NMuKCbmWXCBd3MLBMu6GZmmWjrkaJmndQ7/6qa69YvOL2NScxawz10M7NMuKCbmWXCBd3MLBMu6GZmmXBBNzPLhAu6mVkmXNDNzDLhgm5mlgkXdDOzTNQt6JKOkfSrissOSe+TNFHScklr0/VB7QhsZmbV1S3oEXFnREyPiOnA8cCjwPeA+cCKiJgKrEjzZmbWIY0OuZwE3B0R9wKzgEVp+SJgdjODmZlZYxot6GcBl6bpyRGxGSBdT2pmMDMza0zpgi5pPHAG8J1GNiBpnqR+Sf0DAwON5jMzs5Ia6aG/Drg5Irak+S2SpgCk663VbhQRCyOiLyL6enp6RpfWzMxqaqSgn83/D7cALAPmpOk5wNJmhTIzs8aVKuiSng2cAlxesXgBcIqktWndgubHMzOzskr9YlFEPAo8Z8iyByi+9WJmZmOAjxQ1M8uEC7qZWSZc0M3MMuGCbmaWCRd0M7NMuKCbmWXCBd3MLBMu6GZmmSh1YJHZWNA7/6ph169fcHqbkpiNTe6hm5llwgXdzCwTHnKxbNQbkhnNbT2cY93APXQzs0y4oJuZZcIF3cwsEy7oZmaZcEE3M8uEC7qZWSZc0M3MMlH2R6InSFoi6Q5JaySdIGmipOWS1qbrg1od1szMaivbQ/8c8MOIeAHwUmANMB9YERFTgRVp3szMOqRuQZd0APAXwCUAEfFERGwHZgGLUrNFwOxWhTQzs/rK9NCPBAaAr0paJenLkvYFJkfEZoB0PanajSXNk9QvqX9gYKBpwc3M7OnKFPRxwHHAFyLiWOARGhheiYiFEdEXEX09PT0jjGlmZvWUKegbgY0RcUOaX0JR4LdImgKQrre2JqKZmZVRt6BHxP3ABknHpEUnAbcDy4A5adkcYGlLEpqZWSllT5/7HuCbksYD64C3U/wzuEzSXOA+4E2tiWhmZmWUKugR8Sugr8qqk5obx8zMRspHipqZZcIF3cwsE/4JOrMShvuJOv88nY0V7qGbmWXCBd3MLBMu6GZmmXBBNzPLhAu6mVkmXNDNzDLhgm5mlgkXdDOzTLigm5llwgXdzCwTPvTfxpThDrHPlU8rYM3iHrqZWSZc0M3MMuEhF2ur3XFIxaxd3EM3M8uEC7qZWSZKDblIWg88DPwB2BURfZImAouBXmA98OaI2NaamNZNPKxi1hmN9NBfHRHTI2Lwx6LnAysiYiqwIs2bmVmHjGbIZRawKE0vAmaPPo6ZmY1U2YIewNWSbpI0Ly2bHBGbAdL1pGo3lDRPUr+k/oGBgdEnNjOzqsp+bfEVEbFJ0iRguaQ7ym4gIhYCCwH6+vpiBBnNzKyEUj30iNiUrrcC3wNmAFskTQFI11tbFdLMzOqrW9Al7Stp/8Fp4LXAr4FlwJzUbA6wtFUhzcysvjJDLpOB70kabP+tiPihpBuByyTNBe4D3tS6mGZmVk/dgh4R64CXVln+AHBSK0KZmVnjfKSomVkmXNDNzDLhgm5mlgkXdDOzTLigm5llwj9wYQ3z2RQb4/1l7eIeuplZJlzQzcwy4SEXs1HykIqNFe6hm5llwgXdzCwTLuhmZplwQTczy4QLuplZJlzQzcwy4YJuZpYJF3Qzs0y4oJuZZaJ0QZe0p6RVkq5M88+XdIOktZIWSxrfuphmZlZPIz309wJrKuYvBD4TEVOBbcDcZgYzM7PGlCrokg4DTge+nOYFvAZYkposAma3IqCZmZVTtof+WeAfgafS/HOA7RGxK81vBA6tdkNJ8yT1S+ofGBgYVVgzM6utbkGX9Hpga0TcVLm4StOodvuIWBgRfRHR19PTM8KYZmZWT5nT574COEPSacDewAEUPfYJksalXvphwKbWxTQzs3rq9tAj4kMRcVhE9AJnAT+JiHOAa4AzU7M5wNKWpTQzs7pG8z30fwI+IOkuijH1S5oTyczMRqKhXyyKiJXAyjS9DpjR/EhmZjYSPlLUzCwTLuhmZplwQTczy4QLuplZJlzQzcwy4YJuZpYJF3Qzs0y4oJuZZcIF3cwsEy7oZmaZcEE3M8tEQ+dyMbOxpXf+VTXXrV9wehuT2FjgHrqZWSZc0M3MMuGCbmaWCRd0M7NMuKCbmWXCBd3MLBMu6GZmmahb0CXtLemXkm6RtFrSx9Ly50u6QdJaSYsljW99XDMzq6VMD/1x4DUR8VJgOnCqpJnAhcBnImIqsA2Y27qYZmZWT92CHoWdaXavdAngNcCStHwRMLslCc3MrJRSh/5L2hO4CTgKuBi4G9geEbtSk43AoTVuOw+YB3DEEUeMNq+VNNwh4eDDwrtFvedxNLf1ayA/pT4UjYg/RMR04DBgBvDCas1q3HZhRPRFRF9PT8/Ik5qZ2bAa+pZLRGwHVgIzgQmSBnv4hwGbmhvNzMwaUeZbLj2SJqTpfYCTgTXANcCZqdkcYGmrQpqZWX1lxtCnAIvSOPoewGURcaWk24FvS/oksAq4pIU5zcysjroFPSJuBY6tsnwdxXi6mZmNAT5S1MwsEy7oZmaZ8E/Q7ab802Vm+XEP3cwsEy7oZmaZcEE3M8uEC7qZWSZc0M3MMuFvudgzjOYMf2bWOe6hm5llwgXdzCwTHnLpYh4aMbNK7qGbmWXCPXSz3ZR/oi4/7qGbmWXCBd3MLBMu6GZmmXBBNzPLRJkfiT5c0jWS1khaLem9aflEScslrU3XB7U+rpmZ1VKmh74L+GBEvBCYCbxL0jRgPrAiIqYCK9K8mZl1SN2CHhGbI+LmNP0wsAY4FJgFLErNFgGzWxXSzMzqa2gMXVIvcCxwAzA5IjZDUfSBSc0OZ2Zm5ZUu6JL2A74LvC8idjRwu3mS+iX1DwwMjCSjmZmVUKqgS9qLoph/MyIuT4u3SJqS1k8Btla7bUQsjIi+iOjr6elpRmYzM6uizLdcBFwCrImIiypWLQPmpOk5wNLmxzMzs7LKnMvlFcDbgNsk/Sot+zCwALhM0lzgPuBNrYloZmZl1C3oEXEtoBqrT2puHDMzGykfKWpmlgkXdDOzTLigm5llwgXdzCwTLuhmZplwQTczy4QLuplZJlzQzcwyUeZIUeuQer/KbmZWyT10M7NMuKCbmWXCQy5mVtVohvzWLzi9iUmsLPfQzcwy4YJuZpYJD7l0kL/FYmbN5B66mVkmXNDNzDLhIRcza7p6w4n+FkxruIduZpaJugVd0lckbZX064plEyUtl7Q2XR/U2phmZlZPmR7614BThyybD6yIiKnAijRvZmYdVLegR8TPgAeHLJ4FLErTi4DZTc5lZmYNGukY+uSI2AyQrifVaihpnqR+Sf0DAwMj3JyZmdXT8g9FI2JhRPRFRF9PT0+rN2dmttsaaUHfImkKQLre2rxIZmY2EiMt6MuAOWl6DrC0OXHMzGykynxt8VLgOuAYSRslzQUWAKdIWguckubNzKyD6h4pGhFn11h1UpOzmJnZKPjQ/xbzGRXNGuPTBoycD/03M8uEC7qZWSZc0M3MMuGCbmaWCX8oambZ2N0/UHUP3cwsEy7oZmaZ8JCLmbWdj89oDffQzcwy4YJuZpYJF3Qzs0y4oJuZZcIF3cwsE/6Wi5l1FX9Dpjb30M3MMtE1PfROHdK7ux9KbGbdwz10M7NMuKCbmWViVEMukk4FPgfsCXw5Isbkj0W38kMUf0BjZjB8LWjX0OyIe+iS9gQuBl4HTAPOljStWcHMzKwxoxlymQHcFRHrIuIJ4NvArObEMjOzRo1myOVQYEPF/Ebg5UMbSZoHzEuzOyXdOYptDnUw8DsAXdjEe22PP2bvUt2cv5uzQ3fn72j2UdaJEWdvQn16XplGoynoqrIsnrEgYiGwcBTbqR1A6o+Ivlbcd6t1c3bo7vzdnB26O7+zt9Zohlw2AodXzB8GbBpdHDMzG6nRFPQbgamSni9pPHAWsKw5sczMrFEjHnKJiF2S3g38iOJri1+JiNVNS1ZOS4Zy2qSbs0N35+/m7NDd+Z29hRTxjGFvMzPrQj5S1MwsEy7oZmaZ6KqCLmmipOWS1qbrg2q0O0LS1ZLWSLpdUm97k1bNVCp7anuApN9K+nw7Mw6nTH5J0yVdJ2m1pFslvaUTWSvynCrpTkl3SZpfZf2zJC1O628YC6+TQSWyfyC9tm+VtEJSqe8pt0u9/BXtzpQUksbM1wHLZJf05rT/V0v6Vrsz1hQRXXMB/gOYn6bnAxfWaLcSOCVN7wc8u1uyp/WfA74FfL7TuRvJDxwNTE3TzwU2AxM6lHdP4G7gSGA8cAswbUib84AvpumzgMWd3s8NZH/14OsaeOdYyV42f2q3P/Az4Hqgr9O5G9j3U4FVwEFpflKncw9euqqHTnFqgUVpehEwe2iDdD6ZcRGxHCAidkbEo+2LWFPd7ACSjgcmA1e3KVdZdfNHxG8iYm2a3gRsBXralvDpypyaovIxLQFOklTtgLl2q5s9Iq6peF1fT3EcyFhR9rQgn6DoKDzWznB1lMn+d8DFEbENICK2tjljTd1W0CdHxGaAdD2pSpujge2SLpe0StKn0onEOq1udkl7AP8J/EObs5VRZt//kaQZFD2cu9uQrZpqp6Y4tFabiNgFPAQ8py3phlcme6W5wA9amqgxdfNLOhY4PCKubGewEsrs+6OBoyX9QtL16ayzY8KY+8UiST8GDqmy6iMl72Ic8ErgWOA+YDFwLnBJM/INpwnZzwP+NyI2dKKj2IT8g/czBfgGMCcinmpGthEoc2qKUqev6IDSuSS9FegDTmxposYMmz91XD5D8Xc51pTZ9+Mohl1eRfHO6OeSXhwR21ucra4xV9Aj4uRa6yRtkTQlIjanolHtrc5GYFVErEu3uQKYSRsKehOynwC8UtJ5FGP/4yXtjIiaHyo1UxPyI+kA4CrgoxFxfYuillHm1BSDbTZKGgccCDzYnnjDKnVaDUknU/yzPTEiHm9TtjLq5d8feDGwMnVcDgGWSTojIvrblrK6sq+b6yPiSeCedMLBqRRHz3dUtw25LAPmpOk5wNIqbW4EDpI0OHb7GuD2NmSrp272iDgnIo6IiF7gfODr7SrmJdTNn04B8T2K3N9pY7ZqypyaovIxnQn8JNKnXB1WN3sasvgScMZYGsNNhs0fEQ9FxMER0Zte69dTPI5OF3Mo97q5guJDaSQdTDEEs66tKWvp9KeyjVwoxjdXAGvT9cS0vI/iF5MG250C3ArcBnwNGN8t2Svan8vY+pZL3fzAW4EngV9VXKZ3MPNpwG8oxvE/kpZ9nKJ4AOwNfAe4C/glcGSn93MD2X8MbKnYz8s6nbmR/EParmSMfMul5L4XcBFFR/E24KxOZx68+NB/M7NMdNuQi5mZ1eCCbmaWCRd0M7NMuKCbmWXCBd3MLBMu6GZmmXBBNzPLxP8BH10TkMSEnFYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f825428cb70>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def lpnorm(X,y,p):\n",
    "    n,m=X.shape\n",
    "    \n",
    "    M = Model(\"lpnorm\")\n",
    "    w = M.variable(m)\n",
    "    t = M.variable()\n",
    "    s = M.variable(n)\n",
    "    \n",
    "    M.objective(ObjectiveSense.Minimize, t)\n",
    "    \n",
    "    M.constraint(t == Expr.sum(s))\n",
    "    \n",
    "    res = X @ w - y\n",
    "    M.constraint(Expr.hstack(s, Var.repeat(t,n), res), Domain.inPPowerCone(1.0/p))\n",
    "    \n",
    "    #M.setLogHandler(sys.stdout)\n",
    "    M.solve()\n",
    "\n",
    "    # Return the weights vector and the residuals\n",
    "    w = w.level()\n",
    "    return w, X.dot(w)-y\n",
    "\n",
    "w, res = lpnorm(X,y,p=2)\n",
    "print(\"|w-1|_2:\", np.linalg.norm(w-np.ones(10)))\n",
    "plt.hist(res,40)\n",
    "plt.title(\"Histogram of residuals (lp norm)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Deadzone-linear\n",
    "\n",
    "This approach steams from the $\\ell_1$ norm approximation, but ignoring small errors. Given a threshold $a\\geq 0$, we define \n",
    "\n",
    "$$ \\phi(r) = \\sum_i \\max(0, |r_i|-a )$$\n",
    "\n",
    "The conic optimization model takes the form:\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize} &  \\sum_i t_i\\\\\n",
    "\\mbox{s.t.} & t \\geq |Xw-y|-a,\\\\\n",
    "& t\\geq 0 \n",
    "\\end{array}\n",
    "$$ \n",
    "\n",
    "Again, removing the absolute value we obtain in compact form:\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize} &  \\sum_i t_i\\\\\n",
    "\\mbox{s.t.} & t \\geq Xw-y-a,\\\\\n",
    "  & t \\geq -(Xw-y)-a,\\\\\n",
    "& t\\geq 0.\n",
    "\\end{array}\n",
    "$$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|w-1|_2: 0.03798164802353423\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAGKBJREFUeJzt3X20XXV95/H3B0JAeQ65CeHxag0ooALeiTAOigRaCpakLSIWabCspkjtjK3OmGmdqVU7E6ZVYC1d2iiO0SoEU5BUZlSMUHwAh8uDKCAGQiAhMbk8RIgUEPnOH/t36eZwzt373nse7vndz2uts85+3p+z77nf8zu/s88+igjMzKz/7dTrAGZm1h4u6GZmmXBBNzPLhAu6mVkmXNDNzDLhgm5mlgkX9ClG0l2STux1jl6S9LuSNkraIemYDu7nBEn3jjH/C5I+1ob9DEoKSTMmsO4RkobHmN+WjE22O+HMbdr/6yT9oBf77mcu6F0kaYOkkxumnSfpe6PjEXFkRNxQsZ2e/rN1wd8D742IPSLi9k7tJCK+GxGHd2r7bfJRiuMxrUTEncB2Sb/T6yz9xAXdXmIKvFAcCtxVZ8EpkLVjJM0D3gp8rddZeuTLwJ/0OkQ/cUGfYsqteEkLJA1LekLSVkmfSIvdmO63p26J4yXtJOlDkh6UtE3SFyXtXdruH6Z5j0r6bw37+bCk1ZL+UdITwHlp3zdJ2i5pi6RPSppZ2l5IulDSOklPSvqopN9I6zwh6cry8g2PsWlWSbtK2gHsDPxI0v0t1g9JfyppHbAuTXu1pOskPSbpXklnlZY/TdLdKefDkj6Qpp8oaVNpuWMk3ZaWWwXsVpr3ondSpRyvSsOnS7o9PfaNkj48xt/4PEnr034ekHROi0VPAW6LiKfrZEzz3ybpjvR3+4Gk15XmLZN0f1r3bkm/W5q3s6S/l/SIpPXA6aV5x6fn2ejtaUkb0rxdJV0iaXO6XSJp1/LxlfT+9HfeIundpe3umvb5UHp+f0bSy0oP5wZg4ej2rIaI8K1LN2ADcHLDtPOA7zVbBrgJODcN7wEcl4YHgQBmlNb7I+A+4JVp2auAL6V5RwA7gP8AzKR4C/+r0n4+nMYXU7zIvwx4A3AcMCPt7x7gfaX9BbAG2As4EngGWJv2vzdwN7CkxXFombW07VeNcRwDuA6YlbLuDmwE3p3yHgs8AhyZlt8CnJCG9wWOTcMnApvS8EzgQeDPgV2AM9Mx+Vizv1NjzrSt16bj9zpgK7C48e+Vsj4BHJ7mzRvN2eRx/h3wqdJ4VcZjgW3AGyleFJdQPJ92TfPfDhyQMr4D+CUwL827APgpcHA6rtfT8BxLy+1CUWj/Zxr/CHAzMAcYAH4AfLR0TJ5Ly+wCnAY8Beyb5l9C8RyaBewJ/PPodkv7ewJ4Xa//d/vl1vMA0+mW/rl2ANtLt6doXdBvBP4GmN2wnRcKRGnaWuDC0vjh6Z99BvDfgctL814OPMuLC/qNFdnfB1xdGg/gTaXxW4EPlsY/DlzSYlsts5a2XVXQTyqNvwP4bsMy/wD8dRp+iOKt+14Ny5zIvxX0NwObAZXm/4CaBb1JxkuAixv/XhQFfTvw+8DLKo75Z4HlpfGqjJ8mFdPS/HuBt7TY/h3AojT8HeCC0rzfbHyOlfZxLbBTGr8fOK00/7eADaXj+68Nz9NtFA0FUbyg/EZp3vHAAw37exh4c7f+R/v95i6X7lscEfuM3oALx1j2fOAw4KeSbpH0tjGWPYCi9TbqQYoCMjfN2zg6IyKeAh5tWH9jeUTSYZK+LunnqRvmfwCzG9bZWhr+1ybje0wga13lvIcCb0zdDNslbQfOAfZP83+fonX4oKR/kXR8i0wPR6oipVy1SHqjpOsljUj6BUWLt/F4ERG/pHgBugDYIulaSa9usdnHKVqudTMeCry/4TgcnNYb7Xa7ozTvqFLGFz1HaPLYJf0JRZH+g4h4vrRe49/ygNL4oxHxXGn8KYrnxQBFw+LWUp5vpOlle1K8AFoNLuhTWESsi4h3UrydvQhYLWl3ipZTo80U/9CjDqF4u7uVosvhoNEZqZ9yv8bdNYx/muIt+PyI2Av4S4pWVTuMlbWuct6NwL+UXyijOEPmPQARcUtELKI4jl8DrmyyvS3AgZLKj/GQ0vAvKQoQAJL258W+QtF9cHBE7A18hhbHKyK+GRGnUHS3/JSiJd7MnRQv6HUzbgT+tuE4vDwiLpd0aNrPe4H9UmPiJ6WMWyiKf7PtIukEijNuFkXEL0qzmv0tN7d4PGWPULzoH1nKundEvNAIkHQARTdTy1NL7cVc0KcwSe+SNJBaQ6OtlF8DI8DzFH3Qoy4H/lzSKyTtQdGiXpVaR6uB35H071V8UPk3VBfnPSn6L3ekFuR72vbAxs46EV8HDpN0rqRd0u3fSXqNpJmSzpG0d0T8iuIx/brJNm6ieFH5j5JmSPo9YEFp/o+AIyUdLWk3im6qsj2BxyLiaUkLgD9oFlTSXElnpBfmZyi64JrlgeJzgmPT/upk/CxwQXq3IEm7pw9r96To6gmK5w7pw8mjSutembZ7kKR9gWWlzAcDq4A/jIifNWS8HPiQpAFJsym69/6xxeN5QXpOfxa4WNKctJ8DJf1WabETge9ExDNV27OCC/rUdipwl4ozPy4Fzo6Ip1OXyd8C309vV48DPg98iaLf/QHgaeDPACLirjR8BUVL7EmKvsyx/lE+QFGUnqT4x1vVxsfVMutERMSTFH2+Z1O0Dn9O8Y5m9OyIc4ENqevoAuBdTbbxLPB7FH3lj1N0i1xVmv8zig/3vk1xZs33GjZxIfARSU9SFLVm7wKg+J97f8r5GPAWWnS7RcRWir7tRTUzDgN/DHwyzb8vLUtE3E3xucZNFO+EXgt8v7S7zwLfpHjhuq28XWAhRffV6tKZLqOnlX4MGKZ4N/HjtG7dLzp9MGW8Of1tvk3xecqocyje6VhNenF3nE0HqVW8naI75YFe57HWJB0BrAQWxDT6Z5X0WmBFRDT7vMNacEGfJlR8424tRVfLxylObTt2OhUJs9y5y2X6WETxNn8zMJ+i+8bF3CwjbqGbmWXCLXQzs0x09cJGs2fPjsHBwW7u0sys7916662PRETjl65eoqsFfXBwkOHhlpd2NjOzJiTV+tayu1zMzDLhgm5mlgkXdDOzTLigm5llwgXdzCwTLuhmZplwQTczy4QLuplZJlzQzcwy0dVvippNVYPLrh1z/oblp3cpidnEuYVuZpYJF3Qzs0y4oJuZZcIF3cwsEy7oZmaZcEE3M8uEC7qZWSYqC7qkwyXdUbo9Iel9kmZJuk7SunS/bzcCm5lZc5UFPSLujYijI+Jo4A3AU8DVwDJgbUTMB9amcTMz65HxdrksBO6PiAeBRcDKNH0lsLidwczMbHzGW9DPBi5Pw3MjYgtAup/TbAVJSyUNSxoeGRmZeFIzMxtT7YIuaSZwBvDV8ewgIlZExFBEDA0MDIw3n5mZ1TSeFvpvA7dFxNY0vlXSPIB0v63d4czMrL7xFPR38m/dLQBrgCVpeAlwTbtCmZnZ+NUq6JJeDpwCXFWavBw4RdK6NG95++OZmVldta6HHhFPAfs1THuU4qwXMzObAvwDFzaljPVDE/6RCbOx+av/ZmaZcAvdsuGfkbPpzi10M7NMuKCbmWXCBd3MLBMu6GZmmXBBNzPLhAu6mVkmXNDNzDLhgm5mlgl/scj6RtUXh8ymO7fQzcwy4Ra6TRtu4Vvu3EI3M8uEC7qZWSZc0M3MMuGCbmaWCRd0M7NM1DrLRdI+wOeAo4AA/gi4F1gFDAIbgLMi4vGOpDTrMf80nvWDui30S4FvRMSrgdcD9wDLgLURMR9Ym8bNzKxHKgu6pL2ANwOXAUTEsxGxHVgErEyLrQQWdyqkmZlVq9NCfyUwAvxvSbdL+pyk3YG5EbEFIN3PabaypKWShiUNj4yMtC24mZm9WJ2CPgM4Fvh0RBwD/JJxdK9ExIqIGIqIoYGBgQnGNDOzKnUK+iZgU0T8MI2vpijwWyXNA0j32zoT0czM6qgs6BHxc2CjpMPTpIXA3cAaYEmatgS4piMJzcyslroX5/oz4MuSZgLrgXdTvBhcKel84CHg7Z2JaGZmddQq6BFxBzDUZNbC9sYxM7OJ8jdFzcwy4YJuZpYJF3Qzs0y4oJuZZcIF3cwsEy7oZmaZcEE3M8uEC7qZWSZc0M3MMuGCbmaWCRd0M7NMuKCbmWXCBd3MLBMu6GZmmXBBNzPLhAu6mVkmXNDNzDLhgm5mlgkXdDOzTNT6TVFJG4AngV8Dz0XEkKRZwCpgENgAnBURj3cmppmZVRlPC/2tEXF0RIz+WPQyYG1EzAfWpnEzM+uRyXS5LAJWpuGVwOLJxzEzs4mqW9AD+JakWyUtTdPmRsQWgHQ/p9mKkpZKGpY0PDIyMvnEZmbWVK0+dOBNEbFZ0hzgOkk/rbuDiFgBrAAYGhqKCWQ0M7MaarXQI2Jzut8GXA0sALZKmgeQ7rd1KqSZmVWrLOiSdpe05+gw8JvAT4A1wJK02BLgmk6FNDOzanW6XOYCV0saXf4rEfENSbcAV0o6H3gIeHvnYpqZWZXKgh4R64HXN5n+KLCwE6EsX4PLru11BLNs+ZuiZmaZcEE3M8tE3dMWzV7gbhOzqcktdDOzTLigm5llwgXdzCwTLuhmZplwQTczy4TPcjHrsMmcFbRh+eltTGK5cwvdzCwTLuhmZplwQTczy4QLuplZJvyhqNkk+VIINlW4hW5mlgkXdDOzTLigm5llwgXdzCwTLuhmZpmofZaLpJ2BYeDhiHibpFcAVwCzgNuAcyPi2c7ENLNmxjrDxpcNmH7G00L/T8A9pfGLgIsjYj7wOHB+O4OZmdn41Crokg4CTgc+l8YFnASsTousBBZ3IqCZmdVTt4V+CfBfgOfT+H7A9oh4Lo1vAg5stqKkpZKGJQ2PjIxMKqyZmbVWWdAlvQ3YFhG3lic3WTSarR8RKyJiKCKGBgYGJhjTzMyq1PlQ9E3AGZJOA3YD9qJose8jaUZqpR8EbO5cTDMzq1LZQo+I/xoRB0XEIHA28J2IOAe4HjgzLbYEuKZjKc3MrNJkzkP/IPAXku6j6FO/rD2RzMxsIsZ1tcWIuAG4IQ2vBxa0P5KZmU2EvylqZpYJF3Qzs0y4oJuZZcIF3cwsEy7oZmaZcEE3M8uEC7qZWSZc0M3MMuGCbmaWCRd0M7NMuKCbmWXCBd3MLBMu6GZmmXBBNzPLhAu6mVkmXNDNzDLhgm5mlgkXdDOzTLigm5llovI3RSXtBtwI7JqWXx0Rfy3pFcAVwCzgNuDciHi2k2HNppvBZdf2OoL1kTot9GeAkyLi9cDRwKmSjgMuAi6OiPnA48D5nYtpZmZVKgt6FHak0V3SLYCTgNVp+kpgcUcSmplZLZVdLgCSdgZuBV4FfAq4H9geEc+lRTYBB7ZYdymwFOCQQw6ZbF6rqeqt+oblp3cpiZl1S60PRSPi1xFxNHAQsAB4TbPFWqy7IiKGImJoYGBg4knNzGxM4zrLJSK2AzcAxwH7SBpt4R8EbG5vNDMzG4/Kgi5pQNI+afhlwMnAPcD1wJlpsSXANZ0KaWZm1er0oc8DVqZ+9J2AKyPi65LuBq6Q9DHgduCyDuY0M7MKlQU9Iu4EjmkyfT1Ff7qZmU0B/qaomVkmXNDNzDLhgm5mlgkXdDOzTLigm5llotZX/y0/Y10awJcFmB58eYj8uIVuZpYJF3Qzs0y4y8UsU/5xjOnHLXQzs0y4oJuZZcIF3cwsEy7oZmaZ8Iei9hL+MM2sP7mFbmaWCRd0M7NMuKCbmWXCBd3MLBMu6GZmmag8y0XSwcAXgf2B54EVEXGppFnAKmAQ2ACcFRGPdy6qNfLZKGZWVqeF/hzw/oh4DXAc8KeSjgCWAWsjYj6wNo2bmVmPVBb0iNgSEbel4SeBe4ADgUXAyrTYSmBxp0KamVm1cfWhSxoEjgF+CMyNiC1QFH1gTot1lkoaljQ8MjIyubRmZtZS7YIuaQ/gn4D3RcQTddeLiBURMRQRQwMDAxPJaGZmNdQq6JJ2oSjmX46Iq9LkrZLmpfnzgG2diWhmZnVUFnRJAi4D7omIT5RmrQGWpOElwDXtj2dmZnXVuTjXm4BzgR9LuiNN+0tgOXClpPOBh4C3dyaimZnVUVnQI+J7gFrMXtjeOGZmNlH+pqiZWSZc0M3MMuEfuDCzcau67MSG5ad3KYmVuYVuZpYJF3Qzs0y4y8XMmvLVPPuPW+hmZplwQTczy4S7XKYwv+U1s/FwC93MLBMu6GZmmXBBNzPLhAu6mVkmXNDNzDLhgm5mlgkXdDOzTPg89B7yeeZm1k5uoZuZZcIF3cwsE5UFXdLnJW2T9JPStFmSrpO0Lt3v29mYZmZWpU4L/QvAqQ3TlgFrI2I+sDaNm5lZD1UW9Ii4EXisYfIiYGUaXgksbnMuMzMbp4n2oc+NiC0A6X5OqwUlLZU0LGl4ZGRkgrszM7MqHf9QNCJWRMRQRAwNDAx0endmZtPWRAv6VknzANL9tvZFMjOziZjoF4vWAEuA5en+mrYlMrO+V/WluQ3LT+9SkumlzmmLlwM3AYdL2iTpfIpCfoqkdcApadzMzHqosoUeEe9sMWthm7OYmdkk+JuiZmaZcEE3M8uEr7ZoZl031oemVR+YTmbd3LmFbmaWCRd0M7NMuMulw/wjFmbWLW6hm5llwgXdzCwTLuhmZplwQTczy4QLuplZJnyWi5llY7pf5dEtdDOzTLiFPkk+z9ysvfw/NXFuoZuZZcIF3cwsE33T5TLdP+wwM6viFrqZWSZc0M3MMjGpLhdJpwKXAjsDn4uIvvyxaH+qbjY9dLLrdir88MaEW+iSdgY+Bfw2cATwTklHtCuYmZmNz2S6XBYA90XE+oh4FrgCWNSeWGZmNl6T6XI5ENhYGt8EvLFxIUlLgaVpdIekeyexz5Z0EbOBRzqx7S7o1+z9mhucvRemfG5d1HLWpLKPsd26Dq2z0GQKuppMi5dMiFgBrJjEfuqFkYYjYqjT++mEfs3er7nB2XuhX3ND/2SfTJfLJuDg0vhBwObJxTEzs4maTEG/BZgv6RWSZgJnA2vaE8vMzMZrwl0uEfGcpPcC36Q4bfHzEXFX25KNX8e7dTqoX7P3a25w9l7o19zQJ9kV8ZJubzMz60P+pqiZWSZc0M3MMtG3BV3SLEnXSVqX7vdtsdwhkr4l6R5Jd0sa7G7SpplqZU/L7iXpYUmf7GbGFlkqc0s6WtJNku6SdKekd/QiaynPqZLulXSfpGVN5u8qaVWa/8Op8PyAWrn/Ij2f75S0VlKt85S7oSp7abkzJYWkKXM6YJ3sks5Kx/4uSV/pdsYxRURf3oD/BSxLw8uAi1osdwNwShreA3h5v2RP8y8FvgJ8sh9yA4cB89PwAcAWYJ8e5d0ZuB94JTAT+BFwRMMyFwKfScNnA6umwHGuk/uto89l4D1TIXfd7Gm5PYEbgZuBoV7nHsdxnw/cDuybxuf0Onf51rctdIrLDKxMwyuBxY0LpGvLzIiI6wAiYkdEPNW9iC1VZgeQ9AZgLvCtLuWqUpk7In4WEevS8GZgGzDQtYQvVufyFOXHtBpYKKnZl+a6qTJ3RFxfei7fTPE9kKmg7iVBPkrRQHi6m+Eq1Mn+x8CnIuJxgIjY1uWMY+rngj43IrYApPs5TZY5DNgu6SpJt0v6u3RRsV6rzC5pJ+DjwH/ucrax1DnmL5C0gKKlc38XsjXT7PIUB7ZaJiKeA34B7NeVdK3VyV12PvB/O5qovsrsko4BDo6Ir3czWA11jvthwGGSvi/p5nTF2SljSv9ikaRvA/s3mfVXNTcxAzgBOAZ4CFgFnAdc1o58Y2lD9guB/xMRG7vZYGxD7tHtzAO+BCyJiOfbkW0C6lyeotYlLLqsdiZJ7wKGgLd0NFF9Y2ZPDZWLKf4Pp5o6x30GRbfLiRTvir4r6aiI2N7hbLVM6YIeESe3midpq6R5EbElFY9mb302AbdHxPq0zteA4+hCQW9D9uOBEyRdSNH3P1PSjoho+SFTO7QhN5L2Aq4FPhQRN3coah11Lk8xuswmSTOAvYHHuhOvpVqX1ZB0MsUL7Vsi4pkuZatSlX1P4CjghtRQ2R9YI+mMiBjuWsrm6j5fbo6IXwEPpIsNzqf45nzP9XOXyxpgSRpeAlzTZJlbgH0ljfbhngTc3YVsVSqzR8Q5EXFIRAwCHwC+2OliXkNl7nQZiKsp8n61i9maqXN5ivJjOhP4TqRPu3qoMnfqtvgH4Iwp1o87ZvaI+EVEzI6IwfTcvpniMfS6mEO958vXKD6QRtJsii6Y9V1NOZZefyo70RtFP+daYF26n5WmD1H8etLocqcAdwI/Br4AzOyX7KXlz2NqnOVSmRt4F/Ar4I7S7egeZj4N+BlFP/5fpWkfoSgiALsBXwXuA/4f8MpeH+eaub8NbC0d4zW9zlw3e8OyNzBFznKpedwFfIKiYfhj4OxeZy7f/NV/M7NM9HOXi5mZlbigm5llwgXdzCwTLuhmZplwQTczy4QLuplZJlzQzcwy8f8BJDve1TLY/PIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f81bf9b1278>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def deadzone(X,y,a):\n",
    "    n,m=X.shape\n",
    "    \n",
    "    M = Model(\"deadzone\")\n",
    "    w = M.variable(m)\n",
    "    t = M.variable(n, Domain.greaterThan(0))\n",
    "    \n",
    "    M.objective(ObjectiveSense.Minimize, Expr.sum(t))\n",
    "    \n",
    "    res = X @ w - y\n",
    "    M.constraint(t >= res - a)\n",
    "    M.constraint(t >= -res - a)\n",
    "    \n",
    "    #M.setLogHandler(sys.stdout)\n",
    "    M.solve()\n",
    "    \n",
    "    # Return the weights vector and the residuals\n",
    "    w = w.level()\n",
    "    return w, X.dot(w)-y\n",
    "\n",
    "w, res = deadzone(X,y,0.05)\n",
    "print(\"|w-1|_2:\", np.linalg.norm(w-np.ones(10)))\n",
    "plt.hist(res,40)\n",
    "plt.title(\"Histogram of residuals (deadzone)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chebyshev Approximation\n",
    "\n",
    "Sometimes referred as the *minmax* approximation, it uses \n",
    "\n",
    "$$\\phi(r) = \\max_i( |r_i| ) = \\|r\\|_{\\infty}$$\n",
    "\n",
    "as penalty function. It is a limit case of $\\ell_p-$norms as $p\\rightarrow \\infty$. The effect is to only look at the largest error in absolute value. The formulation of our problem is therefore\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize} & t\\\\\n",
    "\\mbox{s.t.} & te \\geq Xw-y,\\\\\n",
    "  & te \\geq -(Xw-y),\n",
    "\\end{array}\n",
    "$$ \n",
    "where $e$ is the all-ones vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|w-1|_2: 0.16455126712218598\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAFmZJREFUeJzt3X+UZGV95/H3RxA08lsGRHAcjZCIbgJmRFzXIyrmsGKEzaLBRXbIspmjJNlNjO6SmF2NuglsEsSzupuQYBxNFAjyayGuIDpRPGD4pSiiIjgIMjIojPxQEeS7f9w7pmy6p6q7q/rHM+/XOX26qu5z636f6u5PPfXUradTVUiS2vK4xS5AkjR+hrskNchwl6QGGe6S1CDDXZIaZLhLUoMM92UgyY1JDlvsOhZTkn+T5PYkDyQ5eILHeXGSr25l+weSvGsMx1mVpJJsP4d9D0xyzRyPu7J/DLeby/6TkuS8JEcsdh0tMdwXWZINSQ6fctsJSa7Ycr2qnlNV64fcz5zDYpn4M+C3qmqnqrp+Ugepqs9U1c9N6v7H5J10j8esVdU3+8fwx2Ouab5OAf7HYhfREsNdI1kCTxpPB24cpeESqHVikuwDvBS4YLFrGaeq+idglySrF7uWVhjuy8Dg6D7JIUmuSXJfkruSnNY3+3T/fXP/svuFSR6X5A+T3JZkU5IPJtl14H7/fb/tu0n+25TjvD3JuUn+Nsl9wAn9sa9MsjnJxiTvTbLDwP1VkpOS3Jzk/iTvTPKz/T73JTlnsP2UPk5ba5IdkzwAbAd8IcktM+xfSX4zyc3Azf1tP5/ksiT3JPlqktcOtH9lki/3dX4ryZv72w9LcsdAu4OTXNe3Oxt4wsC2n3qFNVDHs/rLRya5vu/77UnevpWf8QlJbu2P840kx83Q9BXAdVX1w4F9NyR5S5IbkjyY5Mwkeyf5WH9/n0iye9/2p17hJVnf/5w+27e9NMmeU9r+el//vUnekOT5/bE2J3nvQB0/m+ST/e/Td5L8XZLdBrbdk+R5/fWn9m0OG+jbeuDImR4jzVJV+bWIX8AG4PApt50AXDFdG+BK4Pj+8k7Aof3lVUAB2w/s9x+ArwPP7NueB3yo33Yg8ADwr4Ad6F7mPzxwnLf314+mGwQ8Efgl4FBg+/54NwG/M3C8Ai4CdgGeAzwEXN4ff1fgy8CaGR6HGWsduO9nbeVxLOAyYI++1icBtwO/3tf7POA7wHP69huBF/eXdwee118+DLijv7wDcBvwu8DjgWP6x+Rd0/2cptbZ39e/6B+/XwDuAo6e+vPqa70P+Ll+2z5b6pymn38KvG+a36GrgL2BfYFNwHXAwcCOwCeBt033e0IXqLcAB/SP23rglClt/4LuSe2XgR/SvWrYa+BYL+nbP4vuyWdHYAXdgOP0gTp/g+535meAjwN/NqUfbwLOW+y/yVa+HLkvDRf0o6DNSTYD/3srbR8GnpVkz6p6oKqu2krb44DTqurWqnoA+H3g2H7Udgzwf6vqiqr6EfDf6f6QB11ZVRdU1aNV9YOquraqrqqqR6pqA/CXwEum7HNqVd1XVTcCXwIu7Y//PeBjdIEz21pH9SdVdU9V/QB4FbChqv6mr/c64KN9v6F7HA9MsktV3dtvn+pQulA/vaoerqpzgatHLaaq1lfVF/vH7wbgIzz28driUeC5SZ5YVRv7x286uwH3T3P7/6qqu6rqW8BngM9V1fVV9RBwPjM/7gB/U1Vf6x+3c4CDpmx/Z1X9sKouBR4EPlJVmwaOdXDf369X1WVV9VBV3Q2cNtjfqvoruldVn6N7AnvrlOPc3/dPY2C4Lw1HV9VuW76Ak7bS9kS6UdZXklyd5FVbaftUupHnFrfRjRT37rfdvmVDVX0f+O6U/W8fvJLkgCQXJ/l2P1Xzx8CeU/a5a+DyD6a5vtMcah3VYL1PB14w5UnzOOAp/fZ/C7wSuC3JPyZ54Qw1fauqBp/0bpum3bSSvCDJp5LcneR7wBt47ONFVT0I/Fq/fWOSS5L8/Ax3ey+w8zS3z/VxB/j2wOXvT9N2pPtOsleSs/pprvuAv+Wx/f0r4Ll0T0YPTdm2M7B5K3VqFgz3Zaaqbq6q19G9LD4VODfJk3jsqBvgTrqQ22Il8AjdH+dGYL8tG5I8EXjy1MNNuf5/gK8A+1fVLsAfAJl7b0audVSD9d4O/OPgk2Z1Z4m8EaCqrq6qo+gexwvoRqxTbQT2TTLYx5UDlx+km2IAIMlT+GkfppumelpV7Uo3vTHt41VVH6+qV9CNaL9CF4LTuYHuyX0p+hO6n8Ev9L8fr2egv0l2Ak4HzgTenmSPKfs/G/jCAtXaPMN9mUny+iQrqupR/nmU82PgbrqX9s8caP4R4HeTPKP/w/pj4OyqegQ4F/iVJP8y3Zucf8TwoN6Zbm74gX5k+caxdWzrtc7FxcABSY5P8vj+6/lJnp1khyTHJdm1qh6m69N0pwZeSfcE85+SbJ/kV4FDBrZ/AXhOkoOSPIHufYpBOwP3VNUPkxwC/LvpCu3f/Hx1/yT9EN17ITOdqngZ8Lz+eEvNznS1b06yL/CWKdvfA1xbVf8RuITuyW7QS+im7jQGhvvycwRwY7ozSN4DHNvPh36f7jzhz/bTEIcC7wc+RPfG1jfo3gz7bYB+Tve3gbPoRqj30705NvWl8qA30wXU/XQjy7PH2K8Za52Lqrqf7g3AY+leFXyb7pXOjn2T44EN/fTBG+hGmVPv40fAr9K9cXov3dTJeQPbvwa8A/gE3VzyFVPu4iTgHUnup3tPY7pXB9D9Hf5eX+c9dCE37dRcVd1F9wbpUTP1fRH9Ed0b19+jC++fPFZJjqL73X1Df9Ob6J6kjuu3Px94sLpTIjUG+enpRG2r+tHyZropl28sdj2aWZIDgXXAIdXIH3CSjwJnVtU/LHYtrTDct2FJfoXuVMUAfw68gO6UQH8ppGXOaZlt21F0UwF3AvvTTfEY7FIDHLlLUoMcuUtSgxZ0gaU999yzVq1atZCHlKRl79prr/1OVa2YzT4LGu6rVq3immvmtAy1JG2zkoz8yegtnJaRpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGLegnVKVhVp18yYzbNpxy5Jz3HWV/qSWO3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJatBI57kn2QDcD/wYeKSqVifZAzgbWAVsAF5bVfdOpkxJ0mzMZuT+0qo6qKpW99dPBi6vqv2By/vrkqQlYD7TMkcB6/rL64Cj51+OJGkcRg33Ai5Ncm2Stf1te1fVRoD++17T7ZhkbZJrklxz9913z79iSdJQo64t86KqujPJXsBlSb4y6gGq6gzgDIDVq1fXHGqUJM3SSCP3qrqz/74JOB84BLgryT4A/fdNkypSkjQ7Q8M9yZOS7LzlMvDLwJeAi4A1fbM1wIWTKlKSNDujTMvsDZyfZEv7D1fV/0tyNXBOkhOBbwKvmVyZkqTZGBruVXUr8IvT3P5d4OWTKEqSND9+QlWSGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1aNTlB6SmrTr5kq1u33DKkRPZV5oUR+6S1CDDXZIaZLhLUoOcc9eCGjY/LWk8HLlLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgz3OXFtnWzv13XRrNlSN3SWqQ4S5JDTLcJalBzrlLE+Z6OloMjtwlqUGGuyQ1yHCXpAY5565txnzmvp0313LjyF2SGmS4S1KDDHdJatDI4Z5kuyTXJ7m4v/6MJJ9LcnOSs5PsMLkyJUmzMZuR+38Gbhq4firw7qraH7gXOHGchUmS5m6kcE+yH3Ak8Nf99QAvA87tm6wDjp5EgZKk2Rt15H468F+AR/vrTwY2V9Uj/fU7gH3HXJskaY6Gnuee5FXApqq6NslhW26epmnNsP9aYC3AypUr51imtG0adn69671rJqOM3F8EvDrJBuAsuumY04Hdkmx5ctgPuHO6navqjKpaXVWrV6xYMYaSJUnDDA33qvr9qtqvqlYBxwKfrKrjgE8Bx/TN1gAXTqxKSdKszGf5gf8KnJXkXcD1wJnjKUmanksASKObVbhX1XpgfX/5VuCQ8ZckSZovP6EqSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGuS/2ZO2US5t0DZH7pLUIMNdkhpkuEtSg5xz19i5Boy0+By5S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIM9z30Zt7Vz0YWuKeB67tPQ5cpekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkho0dG2ZJE8APg3s2Lc/t6reluQZwFnAHsB1wPFV9aNJFitpaRi2vtCw9Yk0eaOM3B8CXlZVvwgcBByR5FDgVODdVbU/cC9w4uTKlCTNxtBwr84D/dXH918FvAw4t799HXD0RCqUJM3aSHPuSbZL8nlgE3AZcAuwuaoe6ZvcAew7mRIlSbM10nruVfVj4KAkuwHnA8+ertl0+yZZC6wFWLly5RzL1EJyvfblYz7r8qttszpbpqo2A+uBQ4Hdkmx5ctgPuHOGfc6oqtVVtXrFihXzqVWSNKKh4Z5kRT9iJ8kTgcOBm4BPAcf0zdYAF06qSEnS7IwyLbMPsC7JdnRPBudU1cVJvgycleRdwPXAmROsU5I0C0PDvapuAA6e5vZbgUMmUZQkaX78hKokNchwl6QGjXQqpJYfT2eUtm2O3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAa5tozUKNcX2rY5cpekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBQ9dzT/I04IPAU4BHgTOq6j1J9gDOBlYBG4DXVtW9kyt12zNsPe4Npxy5QJVIWm5GGbk/AvxeVT0bOBT4zSQHAicDl1fV/sDl/XVJ0hIwNNyramNVXddfvh+4CdgXOApY1zdbBxw9qSIlSbMzqzn3JKuAg4HPAXtX1UbongCAvcZdnCRpbkb+H6pJdgI+CvxOVd2XZNT91gJrAVauXDmXGiUtAv8H6/I20sg9yePpgv3vquq8/ua7kuzTb98H2DTdvlV1RlWtrqrVK1asGEfNkqQhhoZ7uiH6mcBNVXXawKaLgDX95TXAheMvT5I0F6NMy7wIOB74YpLP97f9AXAKcE6SE4FvAq+ZTImSpNkaGu5VdQUw0wT7y8dbjmbDOVFJM/ETqpLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaNPKSv5qbrS0R4L/JU6vmszSGfxfj4chdkhpkuEtSgwx3SWqQc+6LyCV7JU2KI3dJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhrkee6SlhTXpRkPR+6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBg0N9yTvT7IpyZcGbtsjyWVJbu6/7z7ZMiVJszHKyP0DwBFTbjsZuLyq9gcu769LkpaIoeFeVZ8G7ply81HAuv7yOuDoMdclSZqHua4ts3dVbQSoqo1J9pqpYZK1wFqAlStXzvFwS5f/B1XSUjTxN1Sr6oyqWl1Vq1esWDHpw0mSmHu435VkH4D++6bxlSRJmq+5hvtFwJr+8hrgwvGUI0kah6Fz7kk+AhwG7JnkDuBtwCnAOUlOBL4JvGaSRcLwuW3XcZY035zY2v7LLWOGhntVvW6GTS8fcy2SpDHxE6qS1CDDXZIa5P9QlaQxWGrz9Y7cJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXI5QdYeh8bljQZ29K/xXTkLkkNMtwlqUGGuyQ1aJuYc5/PPNu2NEcnaWbLLQscuUtSgwx3SWqQ4S5JDWpmzn25zYdJ0iQ5cpekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1aF7hnuSIJF9N8vUkJ4+rKEnS/Mw53JNsB7wP+NfAgcDrkhw4rsIkSXM3n5H7IcDXq+rWqvoRcBZw1HjKkiTNx3zWltkXuH3g+h3AC6Y2SrIWWNtffSDJd4HvzOO4S92etNs/+7Z8tdy/Jd+3nDrnXbf07emz3XE+4Z5pbqvH3FB1BnDGT3ZKrqmq1fM47pLWcv/s2/LVcv/s2/TmMy1zB/C0gev7AXfO4/4kSWMyn3C/Gtg/yTOS7AAcC1w0nrIkSfMx52mZqnokyW8BHwe2A95fVTeOsOsZw5ssay33z74tXy33z75NI1WPmSaXJC1zfkJVkhpkuEtSgyYe7kn2SHJZkpv777vP0G5lkkuT3JTky0lWTbq2cRi1f33bXZJ8K8l7F7LGuRqlb0kOSnJlkhuT3JDk1xaj1lENWzIjyY5Jzu63f265/B7CSH17U/+3dUOSy5PM+tzpxTTqcidJjklSSZbN6ZGj9C3Ja/uf341JPjz0Tqtqol/A/wRO7i+fDJw6Q7v1wCv6yzsBPzPp2hayf/329wAfBt672HWPq2/AAcD+/eWnAhuB3Ra79hn6sx1wC/BMYAfgC8CBU9qcBPxFf/lY4OzFrnuMfXvplr8r4I3LpW+j9q9vtzPwaeAqYPVi1z3Gn93+wPXA7v31vYbd70JMyxwFrOsvrwOOntqgX5Nm+6q6DKCqHqiq7y9AbeMwtH8ASX4J2Bu4dIHqGoehfauqr1XVzf3lO4FNwIoFq3B2RlkyY7DP5wIvTzLdB/aWmqF9q6pPDfxdXUX32ZTlYtTlTt5JNyj54UIWN0+j9O03gPdV1b0AVbVp2J0uRLjvXVUbAfrve03T5gBgc5Lzklyf5E/7hcmWg6H9S/I44M+BtyxwbfM1ys/uJ5IcQjfyuGUBapuL6ZbM2HemNlX1CPA94MkLUt38jNK3QScCH5toReM1tH9JDgaeVlUXL2RhYzDKz+4A4IAkn01yVZIjht3pfJYf+IkknwCeMs2mt454F9sDLwYOBr4JnA2cAJw5jvrmawz9Own4h6q6fakNAsfQty33sw/wIWBNVT06jtomYJQlM0ZaVmMJGrnuJK8HVgMvmWhF47XV/vUDqHfT5cZyM8rPbnu6qZnD6F5xfSbJc6tq80x3OpZwr6rDZ9qW5K4k+1TVxj4Apns5cQdwfVXd2u9zAXAoSyTcx9C/FwIvTnIS3fsJOyR5oKoWfQ38MfSNJLsAlwB/WFVXTajUcRhlyYwtbe5Isj2wK3DPwpQ3LyMtB5LkcLon7pdU1UMLVNs4DOvfzsBzgfX9AOopwEVJXl1V1yxYlXMz6u/lVVX1MPCNJF+lC/urZ7rThZiWuQhY019eA1w4TZurgd2TbJmrfRnw5QWobRyG9q+qjquqlVW1Cngz8MGlEOwjGNq3fumJ8+n69PcLWNtcjLJkxmCfjwE+Wf07WEvc0L710xZ/Cbx6lDnbJWar/auq71XVnlW1qv87u4qun0s92GG038sL6N4QJ8medNM0t271XhfgneAnA5cDN/ff9+hvXw389UC7VwA3AF8EPgDssNjvYo+zfwPtT2D5nC0ztG/A64GHgc8PfB202LVvpU+vBL5G977AW/vb3kEXBABPAP4e+DrwT8AzF7vmMfbtE8BdAz+nixa75nH2b0rb9SyTs2VG/NkFOI1u0PtF4Nhh9+nyA5LUID+hKkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg/4/2eyUIU3uKWsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f81bf9727b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def minmax(X,y):\n",
    "    n,m=X.shape\n",
    "    \n",
    "    M = Model(\"minmax\")\n",
    "    w = M.variable(m)\n",
    "    t = M.variable()\n",
    "    \n",
    "    M.objective(ObjectiveSense.Minimize, t)\n",
    "    \n",
    "    res = X @ w - y\n",
    "    t_repeat = Var.repeat(t,n)\n",
    "    M.constraint(t_repeat >= res)\n",
    "    M.constraint(t_repeat >= -res)\n",
    "    \n",
    "    #M.setLogHandler(sys.stdout)\n",
    "    M.solve()\n",
    "    \n",
    "    # Return the weights vector and the residuals\n",
    "    w = w.level()\n",
    "    return w, X.dot(w)-y\n",
    "\n",
    "w, res = minmax(X,y)\n",
    "print(\"|w-1|_2:\", np.linalg.norm(w-np.ones(10)))\n",
    "plt.hist(res,40)\n",
    "plt.title(\"Histogram of residuals (minmax)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sum of the largest $k$ residuals\n",
    "\n",
    "The penalty function $\\phi$ is defined as \n",
    "\n",
    "$$\\phi(u) = \\sum_{i=1}^k | r_{[i]} |, $$\n",
    "\n",
    "where the notation $r_{[i]}$ indicates the $i$-th element of $r$ sorted. It means that the penalty depends on the $k$ largest (in absolute value) residuals. This is a generalization of the minmax norm, which here corresponds to $k=1$, and of the $\\ell_1$-norm, which corresponds to $k=n$. As shown in https://docs.mosek.com/modeling-cookbook/linear.html#sum-of-largest-elements whis problem can be modeled in a linear fashion:\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize} & ks+\\sum u_i\\\\\n",
    "\\mbox{s.t.} & se + u \\geq Xw -y\\\\\n",
    "            & se + u \\geq -(Xw -y)\\\\\n",
    "            & u\\geq 0,\\\\\n",
    "            & u\\in \\mathbb{R}^n,s\\in \\mathbb{R},w\\in\\mathbb{R}^m.\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|w-1|_2: 0.08698784890647508\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAGKpJREFUeJzt3XmUXGWZx/HvT0JAZQ10YgAxIgGNjBJsI446IIsHQUlGUXGQCQ5jjqKOuI1xmRHF0TAqwhkdNYoSFzYjQgZciJGoeAjSEEBDwEAMJhCTRggkKEvgmT/u21p0qrtuVdfWb/8+5/Spu7z33ue9Vf3UW+/dFBGYmdno95ROB2BmZs3hhG5mlgkndDOzTDihm5llwgndzCwTTuhmZplwQu8yklZIOrzTcXSSpH+UtFbSFknTW7idV0i6fZj550v6VBO2M0VSSBrXwLLTJPWNNIY6t9lwvE2O41JJx3QyhtHGCb2NJK2RdNSgaadIumZgPCKeHxFLa6ynK/7hWuhzwLsiYqeIWN6qjUTELyPiwFatv0nOpNgfY9E84L86HcRo4oRu2+iCL4pnASvKFOyCWFtG0mTglcBlnY6lEyLi18Aukno7Hcto4YTeZSpb8ZJmSOqT9KCkDZLOTsV+kV43pW6Jl0p6iqSPSbpL0kZJ35K0a8V6/znN+5Ok/xi0nTMkLZT0HUkPAqekbV8raZOk9ZK+KGl8xfpC0mmSVknaLOlMSc9Jyzwo6ZLK8oPqWDVWSTtI2gJsB9ws6c4hlg9J75S0CliVpj1X0mJJ90m6XdIbK8ofK+nWFOfdkj6Qph8uaV1FuemSbkzlLgZ2rJj3pF9SFXHsn4aPk7Q81X2tpDOGeY9PkbQ6bef3kk4aoujRwI0R8XDFsh9Kddic6nlkmv6k7qEqdVsj6YOSbpH0kKTzJE2S9KO0rp9K2n2IeF+flj+o4tfhW1M975f0dkkvTuveJOmLFcs+R9LP0ufuXknflbRbxbz7JB2SxvdKZQ6v2PxS4Lih9qUNEhH+a9MfsAY4atC0U4BrqpUBrgVOTsM7AYem4SlAAOMqlvsX4A5gv1T2UuDbad40YAvwcmA8xU/4xyq2c0Yan0XxJf9U4EXAocC4tL2VwOkV2wtgEbAL8HzgEWBJ2v6uwK3A7CH2w5CxVqx7/2H2YwCLgQkp1qcDa4G3pngPAe4Fnp/KrwdekYZ3Bw5Jw4cD69LweOAu4L3A9sAJaZ98qtr7NDjOtK6/S/vvBcAGYNbg9yvF+iBwYJo3eSDOKvX8LPClivEDUz33qljvc9Lw+QOxDq5bxedqGTAJ2BvYCNwITAd2AH4GfLxKvG9N79X+g+Z9heIL71XAwxS/IiZWrPuwVH5/ii+mHYAeisbIORVxvY3is/U04CfA5wbtg/cBl3b6f3e0/LmF3n6XpVbMJkmbgP8dpuxjwP6S9oyILRGxbJiyJwFnR8TqiNgCfBg4UUWXxAnA/0XENRHxKPCfFP+Ula6NiMsi4omI+EtE3BARyyJia0SsAb4KHDZombMi4sGIWAH8Frgqbf8B4EcUyaLeWMv6TETcFxF/AV4DrImIb6Z4bwS+n+oNxX6cJmmXiLg/zR/sUIpEfk5EPBYRC4HrywYTEUsj4jdp/90CXMi2+2vAE8BBkp4aEevT/qtmN2BzxfjjFIlxmqTtI2JNRFT9FTOE/4mIDRFxN/BL4LqIWB4RjwA/YNv363Tgg8DhEXHHoHlnRsTDEXEV8BBwYURsrFj3dICIuCMiFkfEIxHRD5xNxX6JiK9R/Mq6juLL7aODtrM57QcrwQm9/WZFxG4Df8Bpw5Q9FTgAuE3S9ZJeM0zZvShamAPuomhhTUrz1g7MiIg/A38atPzayhFJB0i6QtIfUzfMp4E9By2zoWL4L1XGd2og1rIq430W8JJBX5QnAc9I818PHAvcJennkl46REx3R0TlF91dVcpVJeklkq6W1C/pAeDtbLu/iIiHgDel+eslXSnpuUOs9n5g54pl76BIsmcAGyVdJGmvsjFS//v1QYpfCOvYVql1SZqY4rw7fY6+w7b75WvAQRRfOI8MmrczsGnoKlklJ/QuFhGrIuLNFD9lzwIWSno627auAe6hSGwD9gW2UvyjrQf2GZgh6anAHoM3N2j8y8BtwNSI2AX4CKDGa1M61rIq410L/LzyizKKM2TeARAR10fETIr9eBlwSZX1rQf2llRZx30rhh+i6BYAQNIzeLILKLqgnhkRu1J0SVTdXxHxk4g4mqJFehtFQqvmFoov9MplL4iIl1Psv6D4XGwTH3/7MhuJVwEfk/T6EazjMxRxviB9jt5CxX6RtBNwDnAecIakCYOWfx5w8wi2P6Y4oXcxSW+R1BMRT/C3VsrjQD/Fz/b9KopfCLxX0rPTP8mngYsjYiuwEHitpL9XcaDyE9ROzjtT9PVuSS3IdzStYsPH2ogrgAMknSxp+/T3YknPkzRe0kmSdo2Ixyjq9HiVdVxL8aXyb5LGSXodMKNi/s3A8yUdLGlHilZypZ2B+yLiYUkzgH+qFmg6EHl8+mJ+hOLYRrV4oDhOcEjaHpIOlHSEpB0o+q3/UrHsTcCxkiakL5vTh9xb5a0AjgG+JOn4BtexM0UdN0nam6LVX+lc4IaI+FfgSoovwkqHUXTfWQlO6N3tGGCFijM/zgVOTP2Wf6Y4P/dXqYvhUOAbwLcpDjr9nuIf/t0AqY/23cBFFC3RzRQHrgb/vK30AYqktJmiBXlxE+s1ZKyNiIjNFK3JEyla/3+kaLnukIqcDKxJP/nfTtFKHLyOR4HXURz8vJ+iW+TSivm/Az4J/JSiz/eaQas4DfikpM0Uxyiq/QqA4n/u/SnO+ygSVtVut4jYQHGwcmaatAPFudn3pjpOpPjlBMX+vJni4OdVNOn9ioibKY5RfE3SqxtYxScoDlI/QJGw/7pPJc2k+Iy/PU16H8UX2Elp/ouBh6I4fdFK0JO7DG0sSK3iTRTdKb/vdDw2NEnTgAXAjBhj/6ySvg+cFxE/7HQso4UT+hgh6bUUpxUK+DzwEorT9/wBMMuEu1zGjpkUP/PvAaZSdN84mZtlxC10M7NMuIVuZpaJtt7YaM8994wpU6a0c5NmZqPeDTfccG9E9NQq19aEPmXKFPr62nprZzOzUU9SqauW3eViZpYJJ3Qzs0w4oZuZZcIJ3cwsE07oZmaZcEI3M8tEzYSebtl5U8Xfg5JOT7fpXKzimZKLh3oeoZmZtUfNhB4Rt0fEwRFxMMVzJv9M8biqucCSiJhKcdOnuS2N1MzMhlVvl8uRwJ0RcRfFzZ4WpOkLKB4wbGZmHVLvlaInUjxtBmBSRKwHiIj1kiZWW0DSHGAOwL777lutiFlbTJl75ZDz1sw7ro2RmLVG6RZ6enTZ8cD36tlARMyPiN6I6O3pqXkrAjMza1A9XS6vBm5Mj8UC2CBpMkB63djs4MzMrLx6Evqb+Vt3CxRPOJ+dhmcDlzcrKDMzq1+phC7pacDRVDzgleJhtUdLWpXmzWt+eGZmVlapg6LpKfN7DJr2J4qzXszMrAv4SlEzs0w4oZuZZcIJ3cwsE07oZmaZcEI3M8uEE7qZWSac0M3MMuGEbmaWCSd0M7NMOKGbmWXCCd3MLBNO6GZmmXBCNzPLhBO6mVkm6n2mqFmWhnveKPiZozY6uIVuZpYJJ3Qzs0y4y8VGDXeLmA3PLXQzs0w4oZuZZaJUQpe0m6SFkm6TtFLSSyVNkLRY0qr0unurgzUzs6GVbaGfC/w4Ip4LvBBYCcwFlkTEVGBJGjczsw6pmdAl7QL8A3AeQEQ8GhGbgJnAglRsATCrVUGamVltZVro+wH9wDclLZf0dUlPByZFxHqA9Dqx2sKS5kjqk9TX39/ftMDNzOzJyiT0ccAhwJcjYjrwEHV0r0TE/IjojYjenp6eBsM0M7NayiT0dcC6iLgujS+kSPAbJE0GSK8bWxOimZmVUTOhR8QfgbWSDkyTjgRuBRYBs9O02cDlLYnQzMxKKXul6LuB70oaD6wG3krxZXCJpFOBPwBvaE2IZmZWRqmEHhE3Ab1VZh3Z3HDMzKxRvlLUzCwTTuhmZplwQjczy4QTuplZJpzQzcwy4YRuZpYJJ3Qzs0w4oZuZZcIJ3cwsE07oZmaZcEI3M8uEE7qZWSbK3m3RrC2mzL2y0yFUNVxca+Yd18ZIzIbmFrqZWSac0M3MMuGEbmaWCSd0M7NMOKGbmWXCCd3MLBNO6GZmmXBCNzPLRKkLiyStATYDjwNbI6JX0gTgYmAKsAZ4Y0Tc35owzcyslnpa6K+MiIMjojeNzwWWRMRUYEkaNzOzDhlJl8tMYEEaXgDMGnk4ZmbWqLIJPYCrJN0gaU6aNiki1gOk14nVFpQ0R1KfpL7+/v6RR2xmZlWVvTnXyyLiHkkTgcWSbiu7gYiYD8wH6O3tjQZiNDOzEkq10CPinvS6EfgBMAPYIGkyQHrd2KogzcystpotdElPB54SEZvT8KuATwKLgNnAvPR6eSsDNcuVb81rzVKmy2US8ANJA+UviIgfS7oeuETSqcAfgDe0LkwzM6ulZkKPiNXAC6tM/xNwZCuCMjOz+vlKUTOzTDihm5llwgndzCwTfki0tVW3PgR6tPIZMlbJLXQzs0w4oZuZZcJdLpYNd+fYWOcWuplZJpzQzcwy4S4XsxGq1dXjs02sXdxCNzPLhBO6mVkm3OVidXMXg1l3cgvdzCwTTuhmZplwQjczy4QTuplZJpzQzcwy4YRuZpYJJ3Qzs0yUTuiStpO0XNIVafzZkq6TtErSxZLGty5MMzOrpZ4Li94DrAR2SeNnAV+IiIskfQU4Ffhyk+MzG/V8W19rl1ItdEn7AMcBX0/jAo4AFqYiC4BZrQjQzMzKKdvlcg7w78ATaXwPYFNEbE3j64C9qy0oaY6kPkl9/f39IwrWzMyGVjOhS3oNsDEibqicXKVoVFs+IuZHRG9E9Pb09DQYppmZ1VKmD/1lwPGSjgV2pOhDPwfYTdK41ErfB7indWGamVktNVvoEfHhiNgnIqYAJwI/i4iTgKuBE1Kx2cDlLYvSzMxqGsl56B8C3ifpDoo+9fOaE5KZmTWirvuhR8RSYGkaXg3MaH5IZmbWCF8pamaWCSd0M7NMOKGbmWXCzxQ162K+bYDVwy10M7NMOKGbmWXCXS62jZH+zHc3gVlnuIVuZpYJJ3Qzs0w4oZuZZcIJ3cwsE07oZmaZcEI3M8uEE7qZWSac0M3MMuGEbmaWCSd0M7NMOKGbmWXCCd3MLBNO6GZmmXBCNzPLRM3b50raEfgFsEMqvzAiPi7p2cBFwATgRuDkiHi0lcFa8/gWt2b5KdNCfwQ4IiJeCBwMHCPpUOAs4AsRMRW4Hzi1dWGamVktNRN6FLak0e3TXwBHAAvT9AXArJZEaGZmpZTqQ5e0naSbgI3AYuBOYFNEbE1F1gF7D7HsHEl9kvr6+/ubEbOZmVVRKqFHxOMRcTCwDzADeF61YkMsOz8ieiOit6enp/FIzcxsWHWd5RIRm4ClwKHAbpIGDqruA9zT3NDMzKweNRO6pB5Ju6XhpwJHASuBq4ETUrHZwOWtCtLMzGqredoiMBlYIGk7ii+ASyLiCkm3AhdJ+hSwHDivhXGamVkNNRN6RNwCTK8yfTVFf7qZmXUBXylqZpYJJ3Qzs0w4oZuZZaLMQVEbhXyvFrOxxy10M7NMOKGbmWXCCd3MLBNO6GZmmXBCNzPLhM9yMRujap0JtWbecW2KxJrFLXQzs0w4oZuZZcJdLmaZ8sVlY49b6GZmmXBCNzPLhBO6mVkmnNDNzDLhhG5mlgmf5TKK+SwGM6vkFrqZWSZqJnRJz5R0taSVklZIek+aPkHSYkmr0uvurQ/XzMyGUqaFvhV4f0Q8DzgUeKekacBcYElETAWWpHEzM+uQmgk9ItZHxI1peDOwEtgbmAksSMUWALNaFaSZmdVWVx+6pCnAdOA6YFJErIci6QMTmx2cmZmVVzqhS9oJ+D5wekQ8WMdycyT1Serr7+9vJEYzMyuhVEKXtD1FMv9uRFyaJm+QNDnNnwxsrLZsRMyPiN6I6O3p6WlGzGZmVkWZs1wEnAesjIizK2YtAman4dnA5c0Pz8zMyipzYdHLgJOB30i6KU37CDAPuETSqcAfgDe0JsSxyxcOWScN9/nz04y6U82EHhHXABpi9pHNDcfMzBrlK0XNzDLhhG5mlgkndDOzTDihm5llwrfP7SCfxWJmzeQWuplZJpzQzcwy4YRuZpYJJ3Qzs0z4oKiZNV2tA/6+dUBruIVuZpYJJ3Qzs0y4y8XM2s53cmwNt9DNzDLhhG5mlgkndDOzTDihm5llwgndzCwTPsvFzEYVnyEzNLfQzcwy4YRuZpaJmgld0jckbZT024ppEyQtlrQqve7e2jDNzKyWMi3084FjBk2bCyyJiKnAkjRuZmYdVDOhR8QvgPsGTZ4JLEjDC4BZTY7LzMzq1Ggf+qSIWA+QXicOVVDSHEl9kvr6+/sb3JyZmdXS8oOiETE/Inojorenp6fVmzMzG7MaTegbJE0GSK8bmxeSmZk1otELixYBs4F56fXypkVkZl2v1hOJrDPKnLZ4IXAtcKCkdZJOpUjkR0taBRydxs3MrINqttAj4s1DzDqyybGYmdkI+F4uZtZV3J3TOF/6b2aWCSd0M7NMuMulxfzz0czaxS10M7NMOKGbmWXCCd3MLBNO6GZmmXBCNzPLhM9yqaHWWSpj/aG0ZtY93EI3M8uEE7qZWSbc5TJCvnDIzLqFW+hmZplwQjczy8SY6HLxmSpmNha4hW5mlgkndDOzTIyJLhczGxtG2r063PKjoWvWLXQzs0w4oZuZZWJEXS6SjgHOBbYDvh4R85oSVRWtvIDHFweZ2Uh1Q3dNwy10SdsBXwJeDUwD3ixpWrMCMzOz+oyky2UGcEdErI6IR4GLgJnNCcvMzOo1ki6XvYG1FePrgJcMLiRpDjAnjW6RdHsD29oTuLeB5bpNDvXIoQ7genSbttRDZ7V82ar1GMl2k2eVKTSShK4q02KbCRHzgfkj2A6S+iKidyTr6AY51COHOoDr0W1cj+YYSZfLOuCZFeP7APeMLBwzM2vUSBL69cBUSc+WNB44EVjUnLDMzKxeDXe5RMRWSe8CfkJx2uI3ImJF0yJ7shF12XSRHOqRQx3A9eg2rkcTKGKbbm8zMxuFfKWomVkmnNDNzDLRlQld0gRJiyWtSq+7D1FuX0lXSVop6VZJU9ob6fDK1iOV3UXS3ZK+2M4YaylTB0kHS7pW0gpJt0h6UydirUbSMZJul3SHpLlV5u8g6eI0/7pu+wwNKFGP96X/gVskLZFU6rzldqtVj4pyJ0gKSV13KmOZOkh6Y3o/Vki6oG3BRUTX/QH/DcxNw3OBs4YotxQ4Og3vBDyt07E3Uo80/1zgAuCLnY673joABwBT0/BewHpgty6IfTvgTmA/YDxwMzBtUJnTgK+k4ROBizsdd4P1eOXA5x94x2itRyq3M/ALYBnQ2+m4G3gvpgLLgd3T+MR2xdeVLXSKWwgsSMMLgFmDC6T7xoyLiMUAEbElIv7cvhBLqVkPAEkvAiYBV7UprnrUrENE/C4iVqXhe4CNQE/bIhxamdtTVNZvIXCkpGoXzXVSzXpExNUVn/9lFNeFdJuytws5k6Ih8XA7gyupTB3eBnwpIu4HiIiN7QquWxP6pIhYD5BeJ1YpcwCwSdKlkpZL+my6YVg3qVkPSU8BPg98sM2xlVXmvfgrSTMoWi53tiG2WqrdnmLvocpExFbgAWCPtkRXXpl6VDoV+FFLI2pMzXpImg48MyKuaGdgdSjzXhwAHCDpV5KWpbvStkXHnlgk6afAM6rM+mjJVYwDXgFMB/4AXAycApzXjPjKakI9TgN+GBFrO9UwbEIdBtYzGfg2MDsinmhGbCNU5vYUpW5h0WGlY5T0FqAXOKylETVm2Hqkxs0XKP6Pu1WZ92IcRbfL4RS/lH4p6aCI2NTi2DqX0CPiqKHmSdogaXJErE9JotpPlnXA8ohYnZa5DDiUNif0JtTjpcArJJ1GcRxgvKQtETHkAaNma0IdkLQLcCXwsYhY1qJQ61Xm9hQDZdZJGgfsCtzXnvBKK3WbDUlHUXwJHxYRj7QptnrUqsfOwEHA0tS4eQawSNLxEdHXtiiHV/YztSwiHgN+n25IOJXi6vqW6tYul0XA7DQ8G7i8Spnrgd0lDfTVHgHc2obY6lGzHhFxUkTsGxFTgA8A32pnMi+hZh3SrR9+QBH799oYWy1lbk9RWb8TgJ9FOpLVRWrWI3VVfBU4vp19tnUath4R8UBE7BkRU9L/wzKK+nRLModyn6nLKA5SI2lPii6Y1W2JrtNHjYc4krwHsARYlV4npOm9FE9GGih3NHAL8BvgfGB8p2NvpB4V5U+h+85yqVkH4C3AY8BNFX8Hdzr2FNuxwO8o+vQ/mqZ9kiJRAOwIfA+4A/g1sF+nY26wHj8FNlTs/0WdjrmRegwqu5QuO8ul5Hsh4GyKBuZvgBPbFZsv/Tczy0S3drmYmVmdnNDNzDLhhG5mlgkndDOzTDihm5llwgndzCwTTuhmZpn4f8iMPPvgNHtBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f81b60b82e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def sumkmax(X,y,k):\n",
    "    n,m=X.shape\n",
    "    \n",
    "    M = Model(\"sumkmax\")\n",
    "    w = M.variable(m)\n",
    "    s = M.variable()\n",
    "    u = M.variable(n, Domain.greaterThan(0))\n",
    "    \n",
    "    M.objective(ObjectiveSense.Minimize, Expr.sum(u) + k * s)\n",
    "    \n",
    "    res = X @ w - y\n",
    "    lhs = u + Var.repeat(s,n)\n",
    "    M.constraint(lhs >= res)\n",
    "    M.constraint(lhs >= -res)\n",
    "    \n",
    "    #M.setLogHandler(sys.stdout)\n",
    "    M.solve()\n",
    "    \n",
    "    # Return the weights vector and the residuals\n",
    "    w = w.level()\n",
    "    return w, X.dot(w)-y\n",
    "\n",
    "w, res = sumkmax(X,y,20)\n",
    "print(\"|w-1|_2:\", np.linalg.norm(w-np.ones(10)))\n",
    "plt.hist(res,40)\n",
    "plt.title(\"Histogram of residuals (sumkmax)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some observations\n",
    "\n",
    "All the regression models presented in this tutorial lead to either a linear or a second-order cone or a power cone  optimization problem. That means that\n",
    "\n",
    "* they can be solved in polynomial time by standard optimization algorithms;\n",
    "* add additional (conic)constraints on the weights can be easily included;\n",
    "* sparsity in the input, i.e. in $X$, can be easily handled and usually leads to great savings in time and memory;\n",
    "* the solver will automatically choose between the primal and dual formulation;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\"><img alt=\"Creative Commons License\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by/4.0/80x15.png\" /></a><br />This work is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\">Creative Commons Attribution 4.0 International License</a>. The **MOSEK** logo and name are trademarks of <a href=\"http://mosek.com\">Mosek ApS</a>. The code is provided as-is. Compatibility with future release of **MOSEK** or the `Fusion API` are not guaranteed. For more information contact our [support](mailto:support@mosek.com). "
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}