{ "cells": [ { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "**********\n", "Regression\n", "**********\n", "\n", ":cite:`sullivanDD3364lmkktc` contains an elegant illustration of the KKT\n", "conditions. Given the constrained optimization problem\n", "\n", ".. math::\n", "\n", " \\min_{\\mathbf{x} \\in \\mathbb{R}^n} \\quad f(\\mathbf{x})\\\\\n", " \\text{subject to} \\quad\n", " h_i(\\mathbf{x}) &= 0 \\quad i = 1, \\ldots, l\\\\\n", " g_j(\\mathbf{x}) &\\leq 0 \\quad j = 1, \\ldots, m,\n", "\n", "define the Lagrangian as\n", "\n", ".. math::\n", "\n", " \\mathcal{L}(\\mathbf{x}, \\boldsymbol{\\mu}, \\boldsymbol{\\lambda}) =\n", " f(\\mathbf{x}) + \\boldsymbol{\\mu}^\\top \\mathbf{h}(\\mathbf{x}) +\n", " \\boldsymbol{\\lambda}^\\top \\mathbf{g}(\\mathbf{x}).\n", "\n", "The KKT optimality conditions are\n", "\n", " Lagrangian Stationarity\n", " :math:`\\nabla_{\\mathbf{x}} \\mathcal{L} = \\boldsymbol{0}`.\n", "\n", " Primal Feasibility\n", " :math:`\\nabla_{\\boldsymbol{\\lambda}} \\mathcal{L} =\n", " \\mathbf{g}(\\mathbf{x}^*) \\preceq \\boldsymbol{0}` and\n", " :math:`\\nabla_{\\boldsymbol{\\mu}} \\mathcal{L} =\n", " \\mathbf{h}(\\mathbf{x}^*) = \\boldsymbol{0}`.\n", "\n", " Dual Feasibility\n", " :math:`\\boldsymbol{\\lambda}^* \\succeq \\boldsymbol{0}`.\n", "\n", " Complementary Slackness\n", " :math:`\\lambda_j^* g_j(\\mathbf{x}^*) = 0` for all :math:`j = 1, \\ldots, m`.\n", "\n", " Positive Definiteness\n", " :math:`\\mathbf{w}^\\top\n", " \\left( \\nabla_{\\mathbf{x} \\mathbf{x}} \\mathcal{L} \\right) \\mathbf{w} > 0`\n", " for all :math:`\\mathbf{w}` in the feasible set.\n", "\n", "Note that non-negativity constraints are merely a subset of inequality\n", "constraints :cite:`lusbysdokkt`. For example,\n", ":math:`\\mathbf{x} \\succeq \\boldsymbol{0}` can be incorporated explicitly by\n", "defining the corresponding Lagrangian as\n", "\n", ".. math::\n", "\n", " \\tilde{\\mathcal{L}}(\\mathbf{x}, \\boldsymbol{\\mu},\n", " \\boldsymbol{\\lambda}, \\boldsymbol{\\nu}) =\n", " \\mathcal{L}(\\mathbf{x}, \\boldsymbol{\\mu}, \\boldsymbol{\\lambda}) +\n", " \\boldsymbol{\\nu}^\\top (-\\mathbf{x})\n", " \\quad \\text{where} \\quad\n", " \\boldsymbol{\\nu} \\succeq \\boldsymbol{0}.\n", "\n", "To incorporate it implicitly, :cite:`wildemcnktt` used\n", "\n", ".. math::\n", "\n", " \\boldsymbol{0}\n", " &= \\nabla_{\\mathbf{x}} \\tilde{\\mathcal{L}}\\\\\n", " &= \\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{x}} - \\boldsymbol{\\nu}\\\\\n", " \\boldsymbol{\\nu} &= \\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{x}}\n", "\n", "to simplify the Lagrange multipliers :math:`\\boldsymbol{\\nu}` in the\n", "(additional) complementary slackness conditions to\n", "\n", ".. math::\n", "\n", " 0 &= -x_k^* \\nu_k^*\\\\\n", " &= -x_k^* \\frac{\\partial \\mathcal{L}}{\\partial x_k}\\\\\n", " &= x_k^* \\frac{\\partial \\mathcal{L}}{\\partial x_k}.\n", "\n", "Furthermore, the Lagrangian stationarity conditions are replaced by\n", ":math:`\\nabla_{\\mathbf{x}} \\mathcal{L} \\succeq \\boldsymbol{0}`, and the\n", "(additional) primal feasibility conditions are typically rewritten as\n", ":math:`\\mathbf{x} \\succeq \\boldsymbol{0}`.\n", "\n", "In some cases, it may seem that the appropriate constraint is a strict\n", "inequality. However, a strict inequality constraint may lead to an optimization\n", "problem that is ill-posed in that a minimizer is infeasible but on the boundary\n", "of the feasible set :cite:`gockenbach2003sic`. For this reason, strict\n", "inequality constraints are not used in nonlinear programming. In practice,\n", "if the strict inequality constraint is not active, one can treat it as a\n", "nonstrict inequality. The alternative is to select the smallest value that is\n", "feasible and use that to establish a nonstrict inequality." ] }, { "cell_type": "raw", "metadata": { "collapsed": true, "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 12.1\n", "=============" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy\n", "\n", "b = numpy.asfarray([0, 3, 1, 2])\n", "_ = numpy.asfarray([0, 1, 2, 4])\n", "A = numpy.vstack((_, numpy.ones(_.shape))).T\n", "\n", "lhs = numpy.dot(A.T, b)\n", "rhs = numpy.linalg.inv(numpy.dot(A.T, A))\n", "print('L2-Regression: x = {}'.format(numpy.dot(rhs, lhs)))" ] }, { "cell_type": "raw", "metadata": { "collapsed": true, "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 12.2\n", "=============" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "m = len(b)\n", "n = 2\n", "\n", "from pymprog import model\n", "\n", "lp = model('L1-Regression LP')\n", "\n", "#lp.verbose(True)\n", "x = lp.var('x', n, bounds=(None, None))\n", "t = lp.var('t', m, bounds=(0, None))\n", "\n", "lp.minimize(sum(t))\n", "\n", "for i in range(m):\n", " -t[i] <= b[i] - sum(A[i, j] * x[j] for j in range(n)) <= t[i]\n", "\n", "lp.solve()\n", "lp.sensitivity()\n", "\n", "lp.end()" ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 12.3\n", "=============\n", "\n", "The minimization of the maximum deviation over a sorted set of real numbers\n", ":math:`\\{ b_1, \\ldots, b_m \\}`\n", "\n", ".. math::\n", "\n", " \\min_{x \\in \\mathbb{R}} \\max_i \\left\\vert x - b_i \\right\\vert\n", "\n", "can be recasted as\n", "\n", ".. math::\n", "\n", " \\mathcal{P} \\quad \\text{minimize} \\quad t\\\\\n", " \\text{subject to} \\quad\n", " x - b_i &\\leq t\\\\\n", " -x + b_i &\\leq t\\\\\n", " t &\\geq 0.\n", "\n", "Recall that :math:`\\left\\vert x - b \\right\\vert` is the distance to :math:`b`\n", "from :math:`x` on the number line. The midrange is defined as\n", "\n", ".. math::\n", "\n", " \\tilde{x} = \\frac{b_1 + b_m}{2}\n", "\n", "where :math:`b_1` and :math:`b_m` are the minimum and maximum values in the\n", "data set. Suppose :math:`\\tilde{x}` is not the minimizer. Then there must\n", "exist :math:`\\epsilon \\neq 0` such that\n", "\n", ".. math::\n", "\n", " \\left\\vert \\tilde{x} + \\epsilon - b_i \\right\\vert <\n", " \\left\\vert \\tilde{x} - b_i \\right\\vert\n", " \\qquad \\forall\\, i = 1, \\ldots, m.\n", "\n", "Define the data points in terms of the midrange: :math:`d_i = \\tilde{x} - b_i`.\n", "The preceding observation must still hold, but\n", "\n", ".. math::\n", "\n", " \\left\\vert \\tilde{x} + \\epsilon - b_i \\right\\vert\n", " &< \\left\\vert \\tilde{x} - b_i \\right\\vert\\\\\n", " \\left\\vert \\epsilon - d_i \\right\\vert\n", " &< \\left\\vert d_i \\right\\vert\n", "\n", "is a contradiction because it can only hold for half of the data points.\n", "Therefore, the midrange minimizes the maximum deviation from the set of\n", "observations." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 12.4\n", "=============\n", "\n", "Given a set of points :math:`\\left\\{ b_i \\in \\mathbb{R}^2 \\right\\}_{i = 1}^m`,\n", "the :math:`\\bar{x}` that solves\n", "\n", ".. math::\n", "\n", " L = \\min_x \\sum_{i = 1}^m \\left\\Vert x - b_i \\right\\Vert_2^2\n", "\n", "satisfies\n", "\n", ".. math::\n", "\n", " \\frac{\\partial L}{\\partial x}\n", " &= \\sum_{i = 1}^m \\frac{\\partial}{\\partial x} (x - b_i)^\\top (x - b_i)\\\\\n", " 0 &= \\sum_{i = 1}^m 2x - 2b_i\\\\\n", " x &= \\frac{1}{m} \\sum_{i = 1}^m b_i." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 12.5\n", "=============\n", "\n", "Given a set of points :math:`\\left\\{ b_i \\in \\mathbb{R}^2 \\right\\}_{i = 1}^m`,\n", "the :math:`\\hat{x}` that solves\n", "\n", ".. math::\n", "\n", " L &= \\min_x \\sum_{i = 1}^m \\left\\Vert x - b_i \\right\\Vert_1\\\\\n", " &= \\min_x \\sum_i \\sum_j \\left\\vert x_j - b_{ij} \\right\\vert\n", "\n", "satisfies\n", "\n", ".. math::\n", "\n", " \\frac{\\partial L}{\\partial x_j}\n", " &= \\sum_{i = 1}^m \\sum_{j'}\n", " \\frac{x_{j'} - b_{ij'}}{\\left\\vert x_{j'} - b_{ij'} \\right\\vert}\n", " \\delta[j - j']\\\\\n", " 0 &= \\sum_i\n", " \\frac{x_j - b_{ij}}{\\left\\vert x_j - b_{ij} \\right\\vert}\\\\\n", " \\sum_i \\frac{x_j}{\\epsilon_{ij}(x)}\n", " &= \\sum_i \\frac{b_{ij}}{\\epsilon_{ij}(x)}\n", " & \\quad & \\epsilon_{ij}(x) = \\left\\vert x_j - b_{ij} \\right\\vert\\\\\n", " x_j &= \\left(\n", " \\sum_i \\frac{1}{\\epsilon_{ij}(x)}\n", " \\right)^{-1}\n", " \\sum_i \\frac{b_{ij}}{\\epsilon_{ij}(x)}\n", " & \\quad & \\forall\\, j = 1, \\ldots, n." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 12.6\n", "=============\n", "\n", "Given a set of points :math:`\\left\\{ b_i \\in \\mathbb{R}^2 \\right\\}_{i = 1}^3`,\n", "the :math:`\\hat{x}` that solves\n", "\n", ".. math::\n", "\n", " L = \\min_x \\sum_i \\left\\Vert x - b_i \\right\\Vert_2\n", "\n", "satisfies\n", "\n", ".. math::\n", "\n", " \\frac{\\partial L}{\\partial x}\n", " &= \\sum_i \\frac{2x - 2b_i}{\\left\\Vert x - b_i \\right\\Vert_2}\\\\\n", " 0 &= \\sum_i \\frac{x}{\\left\\Vert x - b_i \\right\\Vert_2} -\n", " \\sum_i \\frac{b_i}{\\left\\Vert x - b_i \\right\\Vert_2}\\\\\n", " x &= \\left(\n", " \\sum_i \\frac{1}{\\epsilon_i(x)}\n", " \\right)^{-1}\n", " \\sum_i \\frac{b_i}{\\epsilon_i(x)}\n", " & \\quad & \\epsilon_i(x) = \\left\\Vert x - b_i \\right\\Vert_2.\n", "\n", "The solution :math:`\\hat{x}` is known as the Fermat point\n", ":cite:`hajja1994advanced`. When one of the angles exceeds :math:`120^\\circ`,\n", "the minimum is at the vertex holding the largest angle, which is opposite the\n", "edge with the greatest length." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 12.7\n", "=============\n", "\n", "Let :math:`a_i` denote the number of hours of labor needed for month :math:`i`,\n", "and let :math:`x` denote the number of hours per month of labor that will be\n", "handled by regular employees. The total labor costs\n", "\n", ".. math::\n", "\n", " \\min_{x > 0} L(x) =\n", " \\min_{x > 0} \\sum_{i = 1}^{12} 25 \\max(a_i - x, 0) + 17.5x\n", "\n", "is a piecewise-linear minimization.\n", "\n", "(i)\n", "---\n", "\n", "The optimization problem can be recasted as a LP\n", "\n", ".. math::\n", "\n", " \\text{minimize} \\quad \\sum_i 17.5x + 25t_i\\\\\n", " \\text{subject to} \\quad\n", " a_i - x &\\leq t_i \\quad i = 1, \\ldots, 12\\\\\n", " 0 &\\leq t_i, x.\n", "\n", "(ii)\n", "----\n", "\n", "Given the projected labor hours by month in Table 12.2, the optimal solution is\n", ":math:`x^* = 340`.\n", "\n", "(iii)\n", "-----\n", "\n", ".. math::\n", "\n", " \\frac{\\partial L}{\\partial x} =\n", " \\sum_{i = 1}^{12} 25 f(x) + 17.5\n", " \\quad \\text{where} \\quad\n", " f(x) =\n", " \\begin{cases}\n", " -1 & \\text{if } a_i > x\\\\\n", " 0 & \\text{otherwise.}\n", " \\end{cases}\n", "\n", "Applying gradient descent with an adaptive learning rate schedule converges to\n", ":math:`x^* \\approx 340`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "a = [390, 420, 340, 320,\n", " 310, 590, 340, 580,\n", " 550, 360, 420, 600]\n", "\n", "from pymprog import model\n", "\n", "lp = model('Piecewise-Linear LP')\n", "\n", "#lp.verbose(True)\n", "x = lp.var('x', bounds=(0, None))\n", "t = lp.var('t', len(a), bounds=(0, None))\n", "\n", "lp.minimize(sum(17.5 * x + 25 * t_i for t_i in t))\n", "\n", "for a_i, t_i in zip(a, t):\n", " a_i - x <= t_i\n", "\n", "lp.solve()\n", "lp.sensitivity()\n", "\n", "lp.end()" ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 12.8\n", "=============\n", "\n", "Minimizing over the square of the :math:`L^2`-norm of\n", "\n", ".. math::\n", "\n", " \\DeclareMathOperator*{\\argmin}{arg\\,min}\n", " \\argmin_{f, g} L(f, g)\n", " &= \\argmin_{f, g}\n", " \\sum_{i = 1}^n \\left( \\frac{1}{f} - \\frac{f}{g} x_i - t_i \\right)^2\\\\\n", " &= \\argmin_{a, b}\n", " \\sum_i \\left( a + b x_i - t_i \\right)^2\n", "\n", "gives\n", "\n", ".. math::\n", "\n", " \\frac{\\partial L}{\\partial a}\n", " &= \\sum_i 2 \\left( a + b x_i - t_i \\right) (1)\\\\\n", " 0 &= \\sum_i a + b x_i - t_i\\\\\n", " a &= \\frac{1}{n} \\sum_i t_i - b x_i\n", "\n", "and\n", "\n", ".. math::\n", "\n", " \\frac{\\partial L}{\\partial b}\n", " &= \\sum_i 2 \\left( a + b x_i - t_i \\right) (x_i)\\\\\n", " 0 &= \\sum_i a x_i + b x_i^2 - t_i x_i\\\\\n", " \\sum_i b x_i^2 &= \\sum_i t_i x_i - a x_i\\\\\n", " b &= \\frac{\\sum_i t_i x_i - a x_i}{\\sum_i x_i^2}.\n", "\n", "These resulting update equations converge to :math:`f = 1 / a^* \\approx 3.00`\n", "and :math:`g = -\\left( a^* b^* \\right)^{-1} \\approx 8.92`.\n", "\n", "The :math:`L^1`-regression converges to :math:`f \\approx 3.125` and\n", ":math:`g \\approx 9.25`. Note that this solver requires the use of an equivalent\n", "formulation where :math:`b = f / g`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "obs = [(-10, 3.72), (-20, 7.06), (-30, 10.46),\n", " (-10, 3.71), (-20, 7.00), (-30, 10.48),\n", " (-10, 3.67), (-20, 7.08), (-30, 10.33)]\n", "\n", "a = 1\n", "b = 1\n", "n = len(obs)\n", "\n", "denom_const = sum(x_i**2 for x_i, _ in obs)\n", "for i in range(256):\n", " a = sum(t_i - b * x_i for x_i, t_i in obs) / n\n", " b = sum(t_i * x_i - a * x_i for x_i, t_i in obs) / denom_const\n", "\n", " f = round(1 / a, 7)\n", " g = round(-1 / (a * b), 7)\n", "\n", "print('L2-Regression: f = {}, g = {}'.format(f, g))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from pymprog import model\n", "\n", "lp = model('L1-Regression LP')\n", "\n", "#lp.verbose(True)\n", "a = lp.var('a', bounds=(0, None))\n", "b = lp.var('b', bounds=(0, None))\n", "s = lp.var('s', n, bounds=(0, None))\n", "\n", "lp.minimize(sum(s))\n", "\n", "for (x_i, t_i), s_i in zip(obs, s):\n", " -s_i <= a - b * x_i - t_i <= s_i\n", "\n", "lp.solve()\n", "lp.sensitivity()\n", "\n", "lp.end()\n", "\n", "f = round(1 / a.primal, 7)\n", "g = round(1 / (a.primal * b.primal), 7)\n", "print('L1-Regression: f = {}, g = {}'.format(f, g))" ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 12.9\n", "=============\n", "\n", "(a)\n", "---\n", "\n", "A minimum of :math:`g_x(z)` is given by\n", "\n", ".. math::\n", "\n", " \\nabla g_x(z)\n", " &= \\frac{\\partial}{\\partial z}\n", " \\left(\n", " b^\\top E_x^{-1} b - 2 b^\\top E_x^{-1} A z + z^\\top A^\\top E_x^{-1} A z\n", " \\right)\n", " & \\quad & E_x = E_x^\\top\\\\\n", " 0 &= 2 A^\\top E_x^{-1} A z - 2 A^\\top E_x^{-1} b\\\\\n", " z &= \\left( A^\\top E_x^{-1} A \\right)^{-1} A^\\top E_x^{-1} b\\\\\n", " &= T(x)\n", "\n", "because\n", "\n", ".. math::\n", "\n", " \\nabla^2 g_x(z)\n", " &= \\frac{\\partial}{\\partial z}\n", " \\left(\n", " 2 A^\\top E_x^{-1} A z - 2 A^\\top E_x^{-1} b\n", " \\right)\\\\\n", " &= 2 A^\\top E_x^{-1} A\n", "\n", "is positive semi-definite. Furthermore, :math:`T(x)` is a global minimum\n", "seeing that :math:`g_x(z)` is convex:\n", "\n", ".. math::\n", "\n", " g_x(\\alpha z_1 + (1 - \\alpha) z_2)\n", " &= b^\\top E_x^{-1} b -\n", " 2 b^\\top E_x^{-1} A \\left( \\alpha z_1 + (1 - \\alpha) z_2 \\right) +\n", " \\left( \\alpha z_1 + (1 - \\alpha) z_2 \\right)^\\top\n", " A^\\top E_x^{-1} A \\left( \\alpha z_1 + (1 - \\alpha) z_2 \\right)\\\\\n", " &= \\alpha g_x(z_1) + (1 - \\alpha) g_x(z_2) -\n", " \\alpha (1 - \\alpha)\n", " \\left\\Vert E_x^{-1 / 2} A (z_1 - z_2) \\right\\Vert_2^2\\\\\n", " &\\leq \\alpha g_x(z_1) + (1 - \\alpha) g_x(z_2).\n", "\n", "(b)\n", "---\n", "\n", ".. math::\n", "\n", " g_x(x)\n", " &= \\sum_{i = 1}^m \\frac{\\epsilon_i^2(x)}{\\epsilon_i(x)}\\\\\n", " &= \\sum_{i = 1}^m \\epsilon_i(x)\\\\\n", " &= f(x).\n", "\n", "The preceding equality implies that\n", "\n", ".. math::\n", "\n", " \\left\\Vert E_x^{-1 / 2} \\left( b - Ax \\right) \\right\\Vert_2^2\n", " &= \\left( b - Ax \\right)^\\top E_x^{-1} \\left( b - Ax \\right)\\\\\n", " &= \\left\\Vert b - Ax \\right\\Vert_1.\n", "\n", "(c)\n", "---\n", "\n", "Let :math:`y = T(x)` and\n", ":math:`\\epsilon_i(y) =\n", "\\epsilon_i(x) + \\left( \\epsilon_i(y) - \\epsilon_i(x) \\right)`.\n", "\n", ".. math::\n", "\n", " g_x(y)\n", " &= \\sum_i \\frac{\\epsilon^2_i(y)}{\\epsilon_i(x)}\\\\\n", " &= \\sum_i \\frac{\n", " \\left[\n", " \\epsilon_i(x) + \\left( \\epsilon_i(y) - \\epsilon_i(x) \\right)\n", " \\right]^2\n", " }{\n", " \\epsilon_i(x)\n", " }\\\\\n", " &= \\sum_i \\frac{\n", " \\epsilon^2_i(x) +\n", " 2 \\epsilon_i(x) \\left( \\epsilon_i(y) - \\epsilon_i(x) \\right) +\n", " \\left( \\epsilon_i(y) - \\epsilon_i(x) \\right)^2\n", " }{\n", " \\epsilon_i(x)\n", " }\\\\\n", " &= \\sum_i 2 \\epsilon_i(y) - \\epsilon_i(x) +\n", " \\frac{\n", " \\left( \\epsilon_i(y) - \\epsilon_i(x) \\right)^2\n", " }{\n", " \\epsilon_i(x)\n", " }\\\\\n", " &= 2 f(y) - f(x) +\n", " \\sum_i \\frac{\n", " \\left( \\epsilon_i(y) - \\epsilon_i(x) \\right)^2\n", " }{\n", " \\epsilon_i(x)\n", " }\\\\\n", " &\\geq 2 f\\left( T(x) \\right) - f(x)\n", "\n", "(d)\n", "---\n", "\n", "Let :math:`y = T(x)`. By inspection, :math:`T(x) \\neq x` for each valid\n", ":math:`x`.\n", "\n", ".. math::\n", "\n", " f(y) &\\leq \\frac{1}{2} \\left[ g_x(y) + f(x) \\right]\n", " & \\quad & \\text{(c)}\\\\\n", " &< \\frac{1}{2} \\left[ g_x(x) + f(x) \\right]\n", " & \\quad & \\text{(a)}\\\\\n", " &< f(x)\n", " & \\quad & \\text{(b)}\n", "\n", "Therefore, the sequence of iterates in the iteratively reweighted least squares\n", "algorithm produces a monotonically decreasing sequence of objective function\n", "values." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 12.10\n", "==============\n", "\n", "Let parameter :math:`\\mu \\in \\mathbb{R}` and the set of numbers\n", ":math:`\\left\\{ b_j \\right\\}_{j = 1}^n`.\n", "\n", "(a)\n", "---\n", "\n", "The optimization problem can be recast as\n", "\n", ".. math::\n", "\n", " \\min_x L\n", " &= \\min_x \\sum_j \\left\\vert x - b_j \\right\\vert + \\mu (x - b_j)\\\\\n", " &= \\min_x \\sum_i^{P} (x - b_i) + \\mu (x - b_i) +\n", " \\sum_k^{N} (b_k - x) + \\mu (x - b_k)\n", " & \\quad & P + N = n \\land x \\geq b_i \\land x < b_k.\n", "\n", "The optimal solution satisfies\n", "\n", ".. math::\n", "\n", " \\frac{\\partial L}{\\partial x}\n", " &= \\left( \\sum_i^{P} 1 + \\mu \\right) + \\left( \\sum_k^{N} -1 + \\mu \\right)\\\\\n", " 0 &= n \\mu + P - N.\n", "\n", "The following analysis shows that this optimization problem can be interpreted\n", "as quantile estimation.\n", "\n", "- :math:`\\mu = 0 \\implies P = N`.\n", "\n", " - This is the median or second quantile.\n", "\n", "- :math:`\\mu = \\frac{1}{2} \\implies P = \\frac{1}{4} n \\land N = \\frac{3}{4} n`.\n", "\n", " - This is the first quantile.\n", "\n", "- :math:`\\mu = -\\frac{1}{2} \\implies P = \\frac{3}{4} n \\land N = \\frac{1}{4} n`.\n", "\n", " - This is the third quantile.\n", "\n", "- :math:`\\mu = 1 \\implies P = 0 \\land N = n`.\n", "\n", " - This is the minimum value of the set or zeroth quantile.\n", "\n", "- :math:`\\mu = -1 \\implies P = n \\land N = 0`.\n", "\n", " - This is the maximum value of the set or fourth quantile.\n", "\n", "By inspection, the minimum does not exist when :math:`\\mu \\not\\in [-1, 1]`.\n", "\n", "(b)\n", "---\n", "\n", ".. math::\n", "\n", " \\DeclareMathOperator*{\\argmax}{arg\\,max}\n", " \\text{minimize} \\quad \\sum_j t_j + \\mu (x - b_j)\\\\\n", " \\text{subject to} \\quad\n", " x - b_j &\\leq t_j\\\\\n", " -x + b_j &\\leq t_j\\\\\n", " t_j &\\geq 0.\n", "\n", "Observe that the constant term :math:`\\sum_j \\mu b_j` in the LP objective does\n", "not affect the :math:`\\argmax L`. Therefore, the LP in standard form is\n", "\n", ".. math::\n", "\n", " \\text{maximize} \\quad -n \\mu x - \\sum_j t_j\\\\\n", " \\text{subject to} \\quad\n", " x - b_j &\\leq t_j\\\\\n", " -x + b_j &\\leq t_j\\\\\n", " t_j &\\geq 0.\n", "\n", "(c)\n", "---\n", "\n", "The LP in standard form can be thought of as a\n", ":ref:`one parameter family of linear programming problems `.\n", "\n", "Suppose :math:`n = 4` and :math:`b_j = 2^{j - 1}`. Rewriting the standard form\n", "LP into matrix form gives\n", "\n", ".. math::\n", "\n", " \\text{maximize} \\quad -n \\mu x - \\sum_j t_j\\\\\n", " \\text{subject to} \\quad\n", " x - b_j &\\leq t_j\\\\\n", " -x + b_j &\\leq t_j\\\\\n", " t_j &\\geq 0.\n", "\n", ".. math::\n", "\n", " \\mathcal{B} &= \\{ 5, 6, 7, 8, 9, 10, 11, 12 \\} \\quad & \\quad\n", " \\mathcal{N} &= \\{ 0, 1, 2, 3, 4 \\}\\\\\n", " B &= \\mathbf{I}_8 \\quad & \\quad\n", " N &= \\begin{bmatrix}\n", " 1 & -1\\\\\n", " -1 & -1\\\\\n", " 1 & & -1\\\\\n", " -1 & & -1\\\\\n", " 1 & & & -1\\\\\n", " -1 & & & -1\\\\\n", " 1 & & & & -1\\\\\n", " -1 & & & & -1\n", " \\end{bmatrix}\\\\\n", " x^*_\\mathcal{B} &= b = \\begin{bmatrix}\n", " 1\\\\ -1\\\\ 2\\\\ -2\\\\ 4\\\\ -4\\\\ 8\\\\ -8\n", " \\end{bmatrix}\n", " \\quad & \\quad\n", " z^*_\\mathcal{N} &= -c_\\mathcal{N} =\n", " \\begin{bmatrix} n\\mu' \\\\ 1\\\\ 1\\\\ 1\\\\ 1 \\end{bmatrix}.\n", "\n", "To analyze the one parameter family of LPs, arbitrarily perturb the initial\n", "dictionary to make it optimal by defining\n", ":math:`\\bar{x}_\\mathcal{B} = \\boldsymbol{1}_8` and\n", ":math:`\\bar{z}_\\mathcal{N} = \\boldsymbol{1}_5`.\n", "\n", "Evaluating :math:`x^*_\\mathcal{B} + \\bar{x}_\\mathcal{B} \\mu \\geq \\boldsymbol{0}`\n", "and :math:`z^*_\\mathcal{N} + \\bar{z}_\\mathcal{N} \\mu \\geq \\boldsymbol{0}`\n", "reveals the current range for :math:`\\mu` is\n", "\n", ".. math::\n", "\n", " 2^{n - 1} \\leq -n\\mu' \\leq \\mu \\leq \\infty.\n", "\n", "By inspection, the initial dictionary is optimal as long as\n", ":math:`\\mu' \\leq -\\frac{8}{n}`.\n", "\n", "Iteration 1\n", "^^^^^^^^^^^\n", "\n", "Step 1\n", " Computing\n", "\n", " .. math::\n", "\n", " \\mu^*\n", " &= \\max\\left\\{\n", " \\max_{j \\in \\mathcal{N}, \\bar{z}_j > 0} -\\frac{z^*_j}{\\bar{z}_j},\n", " \\max_{i \\in \\mathcal{B}, \\bar{x}_i > 0} -\\frac{x^*_i}{\\bar{x}_i},\n", " \\right\\}\\\\\n", " &= \\max\\left\\{\n", " \\max\\left\\{\n", " -\\frac{n\\mu'}{1},\n", " -\\frac{1}{1}, -\\frac{1}{1}, -\\frac{1}{1}, -\\frac{1}{1}\n", " \\right\\},\n", " \\max\\left\\{\n", " -\\frac{1}{1}, -\\frac{-1}{1},\n", " -\\frac{2}{1}, -\\frac{-2}{1},\n", " -\\frac{4}{1}, -\\frac{-4}{1},\n", " -\\frac{8}{1}, -\\frac{-8}{1}\n", " \\right\\}\n", " \\right\\}\\\\\n", " &= -n\\mu'\n", "\n", " shows that :math:`x_0` is the entering variable for\n", " :math:`\\mu' < -\\frac{8}{n}`.\n", "\n", "Step 2\n", " The maximum is achieved by :math:`j = 0 \\in \\mathcal{N}`, so do one step of\n", " the primal simplex method.\n", "\n", " The primal step direction is\n", "\n", " .. math::\n", "\n", " \\Delta x_\\mathcal{B} =\n", " B^{-1} N e_j =\n", " \\begin{bmatrix}\n", " 1 & -1\\\\\n", " -1 & -1\\\\\n", " 1 & & -1\\\\\n", " -1 & & -1\\\\\n", " 1 & & & -1\\\\\n", " -1 & & & -1\\\\\n", " 1 & & & & -1\\\\\n", " -1 & & & & -1\n", " \\end{bmatrix} \\begin{bmatrix} 1\\\\ 0\\\\ 0\\\\ 0\\\\ 0 \\end{bmatrix} =\n", " \\begin{bmatrix} 1\\\\ -1\\\\ 1\\\\ -1\\\\ 1\\\\ -1\\\\ 1\\\\ -1 \\end{bmatrix}.\n", "\n", " The index of the leaving variable is\n", "\n", " .. math::\n", "\n", " \\argmin_{i \\in \\mathcal{B}}\n", " \\left\\{\n", " \\frac{x^*_i + \\mu^* \\bar{x}_i}{\\Delta x_i} \\colon \\Delta x_i > 0\n", " \\right\\}\n", " &= \\argmin_{i \\in \\mathcal{B}}\\left\\{\n", " \\frac{1 + \\mu^* (1)}{1},\n", " \\frac{2 + \\mu^* (1)}{1},\n", " \\frac{4 + \\mu^* (1)}{1},\n", " \\frac{8 + \\mu^* (1)}{1}\n", " \\right\\}\\\\\n", " &= 5.\n", "\n", " The dual step direction is\n", "\n", " .. math::\n", "\n", " \\Delta z_\\mathcal{N} =\n", " -\\left( B^{-1} N \\right)^\\top e_i =\n", " -\\begin{bmatrix}\n", " 1 & -1 & 1 & -1 & 1 & -1 & 1 & -1\\\\\n", " -1 & -1\\\\\n", " & & -1 & -1\\\\\n", " & & & & -1 & -1\\\\\n", " & & & & & &-1 & -1\n", " \\end{bmatrix} \\begin{bmatrix}\n", " 1\\\\ 0\\\\ 0\\\\ 0\\\\ 0\\\\ 0\\\\ 0\\\\ 0\n", " \\end{bmatrix} =\n", " \\begin{bmatrix} -1\\\\ 1\\\\ 0\\\\ 0\\\\ 0 \\end{bmatrix}.\n", "\n", "Step 3\n", " The step length adjustments are given by\n", "\n", " .. math::\n", "\n", " t &= \\frac{x^*_i}{\\Delta x_i} = \\frac{1}{1} = 1\n", " \\quad & \\quad\n", " \\bar{t} = \\frac{\\bar{x}_i}{\\Delta x_i} = \\frac{1}{1} = 1\\\\\n", " s &= \\frac{z^*_j}{\\Delta z_j} = \\frac{n\\mu'}{-1} = -n\\mu'\n", " \\quad & \\quad\n", " \\bar{s} = \\frac{\\bar{z}_j}{\\Delta z_j} = \\frac{1}{-1} = -1.\n", "\n", "Step 4\n", " The current primal and dual solutions are updated to\n", "\n", " .. math::\n", "\n", " x^*_\\mathcal{B} &\\leftarrow x^*_\\mathcal{B} - t \\Delta x_\\mathcal{B} =\n", " \\begin{bmatrix} 1\\\\ -1\\\\ 2\\\\ -2\\\\ 4\\\\ -4\\\\ 8\\\\ -8 \\end{bmatrix} -\n", " (1) \\begin{bmatrix} 1\\\\ -1\\\\ 1\\\\ -1\\\\ 1\\\\ -1\\\\ 1\\\\ -1 \\end{bmatrix} =\n", " \\begin{bmatrix} 0\\\\ 0\\\\ 1\\\\ -1\\\\ 3\\\\ -3\\\\ 7\\\\ -7 \\end{bmatrix}\n", " \\quad & \\quad\n", " x^*_j &\\leftarrow t = 1\\\\\n", " \\bar{x}_\\mathcal{B} &\\leftarrow\n", " \\bar{x}_\\mathcal{B} - \\bar{t} \\Delta x_\\mathcal{B} =\n", " \\begin{bmatrix} 1\\\\ 1\\\\ 1\\\\ 1\\\\ 1\\\\ 1\\\\ 1\\\\ 1 \\end{bmatrix} -\n", " (1) \\begin{bmatrix} 1\\\\ -1\\\\ 1\\\\ -1\\\\ 1\\\\ -1\\\\ 1\\\\ -1 \\end{bmatrix} =\n", " \\begin{bmatrix} 0\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2 \\end{bmatrix}\n", " \\quad & \\quad\n", " \\bar{x}_j &\\leftarrow \\bar{t} = 1\n", "\n", " and\n", "\n", " .. math::\n", "\n", " z^*_\\mathcal{N} &\\leftarrow z^*_\\mathcal{N} - s \\Delta z_\\mathcal{N} =\n", " \\begin{bmatrix} n\\mu' \\\\ 1\\\\ 1\\\\ 1\\\\ 1 \\end{bmatrix} -\n", " (-n\\mu') \\begin{bmatrix} -1\\\\ 1\\\\ 0\\\\ 0\\\\ 0 \\end{bmatrix} =\n", " \\begin{bmatrix} 0\\\\ 1 + n\\mu'\\\\ 1\\\\ 1\\\\ 1 \\end{bmatrix}\n", " \\quad & \\quad\n", " z^*_i &\\leftarrow s = -n\\mu'\\\\\n", " \\bar{z}_\\mathcal{N} &\\leftarrow\n", " \\bar{z}_\\mathcal{N} - \\bar{s} \\Delta z_\\mathcal{N} =\n", " \\begin{bmatrix} 1\\\\ 1\\\\ 1\\\\ 1\\\\ 1 \\end{bmatrix} -\n", " (-1) \\begin{bmatrix} -1\\\\ 1\\\\ 0\\\\ 0\\\\ 0 \\end{bmatrix} =\n", " \\begin{bmatrix} 0\\\\ 2\\\\ 1\\\\ 1\\\\ 1 \\end{bmatrix}\n", " \\quad & \\quad\n", " \\bar{z}_i &\\leftarrow \\bar{s} = -1.\n", "\n", "Step 5\n", " The current basis is updated to\n", "\n", " .. math::\n", "\n", " \\mathcal{B} &= \\{ 0, 6, 7, 8, 9, 10, 11, 12 \\} \\quad & \\quad\n", " \\mathcal{N} &= \\{ 5, 1, 2, 3, 4 \\}\\\\\n", " B &= \\begin{bmatrix}\n", " 1\\\\\n", " -1 & 1\\\\\n", " 1 & & 1\\\\\n", " -1 & & & 1\\\\\n", " 1 & & & & 1\\\\\n", " -1 & & & & & 1\\\\\n", " 1 & & & & & & 1\\\\\n", " -1 & & & & & & & 1\n", " \\end{bmatrix} \\quad & \\quad\n", " N &= \\begin{bmatrix}\n", " 1 & -1\\\\\n", " & -1\\\\\n", " & & -1\\\\\n", " & & -1\\\\\n", " & & & -1\\\\\n", " & & & -1\\\\\n", " & & & & -1\\\\\n", " & & & & -1\n", " \\end{bmatrix}\\\\\n", " x^*_\\mathcal{B} &= \\begin{bmatrix}\n", " 1\\\\ 0\\\\ 1\\\\ -1\\\\ 3\\\\ -3\\\\ 7\\\\ -7\n", " \\end{bmatrix}\n", " \\quad & \\quad\n", " z^*_\\mathcal{N} &= \\begin{bmatrix}\n", " -n\\mu' \\\\ 1 + n\\mu'\\\\ 1\\\\ 1\\\\ 1\n", " \\end{bmatrix}\\\\\n", " \\bar{x}_\\mathcal{B} &= \\begin{bmatrix}\n", " 1\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2\n", " \\end{bmatrix}\n", " \\quad & \\quad\n", " \\bar{z}_\\mathcal{N} &= \\begin{bmatrix} -1\\\\ 2\\\\ 1\\\\ 1\\\\ 1 \\end{bmatrix}.\n", "\n", " The current range for :math:`\\mu` is\n", "\n", " .. math::\n", "\n", " \\frac{-1 - n\\mu'}{2} \\leq \\mu \\leq -n\\mu'.\n", "\n", " By inspection, the dictionary is optimal as long as\n", " :math:`\\mu' \\leq -\\frac{8}{n}`.\n", "\n", "Iteration 2\n", "^^^^^^^^^^^\n", "\n", "Step 1\n", " Computing\n", "\n", " .. math::\n", "\n", " \\mu^*\n", " &= \\max\\left\\{\n", " \\max_{j \\in \\mathcal{N}, \\bar{z}_j > 0} -\\frac{z^*_j}{\\bar{z}_j},\n", " \\max_{i \\in \\mathcal{B}, \\bar{x}_i > 0} -\\frac{x^*_i}{\\bar{x}_i},\n", " \\right\\}\\\\\n", " &= \\max\\left\\{\n", " \\max\\left\\{\n", " -\\frac{1 + n\\mu'}{2},\n", " -\\frac{1}{1}, -\\frac{1}{1}, -\\frac{1}{1}\n", " \\right\\},\n", " \\max\\left\\{\n", " -\\frac{1}{1}, -\\frac{0}{2}, -\\frac{-1}{2},\n", " -\\frac{-3}{2}, -\\frac{-7}{2}\n", " \\right\\}\n", " \\right\\}\\\\\n", " &= -\\frac{1 + n\\mu'}{2}\n", "\n", " shows that :math:`x_1` is the entering variable for\n", " :math:`\\mu' < -\\frac{8}{n}`.\n", "\n", "Step 2\n", " The maximum is achieved by :math:`j = 1 \\in \\mathcal{N}`, so do one step of\n", " the primal simplex method.\n", "\n", " The primal step direction is\n", "\n", " .. math::\n", "\n", " \\Delta x_\\mathcal{B} =\n", " B^{-1} N e_j =\n", " \\begin{bmatrix}\n", " 1 & -1\\\\\n", " 1 & -2\\\\\n", " -1 & 1 & -1\\\\\n", " 1 & -1 & -1\\\\\n", " -1 & 1 & 0 & -1\\\\\n", " 1 & -1 & 0 & -1\\\\\n", " -1 & 1 & 0 & 0 & -1\\\\\n", " 1 & -1 & 0 & 0 & -1\n", " \\end{bmatrix} \\begin{bmatrix} 0\\\\ 1\\\\ 0\\\\ 0\\\\ 0 \\end{bmatrix} =\n", " \\begin{bmatrix} -1\\\\ -2\\\\ 1\\\\ -1\\\\ 1\\\\ -1\\\\ 1\\\\ -1 \\end{bmatrix}.\n", "\n", " The index of the leaving variable is\n", "\n", " .. math::\n", "\n", " \\argmin_{i \\in \\mathcal{B}}\n", " \\left\\{\n", " \\frac{x^*_i + \\mu^* \\bar{x}_i}{\\Delta x_i} \\colon \\Delta x_i > 0\n", " \\right\\}\n", " &= \\argmin_{i \\in \\mathcal{B}}\\left\\{\n", " \\frac{1 + \\mu^* (0)}{1},\n", " \\frac{3 + \\mu^* (0)}{1},\n", " \\frac{7 + \\mu^* (0)}{1}\n", " \\right\\}\\\\\n", " &= 7.\n", "\n", " The dual step direction is\n", "\n", " .. math::\n", "\n", " \\Delta z_\\mathcal{N} =\n", " -\\left( B^{-1} N \\right)^\\top e_i =\n", " -\\begin{bmatrix}\n", " 1 & 1 & -1 & 1 & -1 & 1 & -1 & 1\\\\\n", " -1 & -2 & 1 & -1 & 1 & -1 & 1 & -1\\\\\n", " & & -1 & -1\\\\\n", " & & & & -1 & -1\\\\\n", " & & & & & & -1 & -1\n", " \\end{bmatrix} \\begin{bmatrix}\n", " 0\\\\ 0\\\\ 1\\\\ 0\\\\ 0\\\\ 0\\\\ 0\\\\ 0\n", " \\end{bmatrix} =\n", " \\begin{bmatrix} 1\\\\ -1\\\\ 1\\\\ 0\\\\ 0 \\end{bmatrix}.\n", "\n", "Step 3\n", " The step length adjustments are given by\n", "\n", " .. math::\n", "\n", " t &= \\frac{x^*_i}{\\Delta x_i} = \\frac{1}{1} = 1\n", " \\quad & \\quad\n", " \\bar{t} = \\frac{\\bar{x}_i}{\\Delta x_i} = \\frac{0}{1} = 0\\\\\n", " s &= \\frac{z^*_j}{\\Delta z_j} = \\frac{1 + n\\mu'}{-1} = -1 - n\\mu'\n", " \\quad & \\quad\n", " \\bar{s} = \\frac{\\bar{z}_j}{\\Delta z_j} = \\frac{2}{-1} = -2.\n", "\n", "Step 4\n", " The current primal and dual solutions are updated to\n", "\n", " .. math::\n", "\n", " x^*_\\mathcal{B} &\\leftarrow x^*_\\mathcal{B} - t \\Delta x_\\mathcal{B} =\n", " \\begin{bmatrix} 1\\\\ 0\\\\ 1\\\\ -1\\\\ 3\\\\ -3\\\\ 7\\\\ -7 \\end{bmatrix} -\n", " (1) \\begin{bmatrix}\n", " -1\\\\ -2\\\\ 1\\\\ -1\\\\ 1\\\\ -1\\\\ 1\\\\ -1\n", " \\end{bmatrix} =\n", " \\begin{bmatrix} 2\\\\ 2\\\\ 0\\\\ 0\\\\ 2\\\\ -2\\\\ 6\\\\ -6 \\end{bmatrix}\n", " \\quad & \\quad\n", " x^*_j &\\leftarrow t = 1\\\\\n", " \\bar{x}_\\mathcal{B} &\\leftarrow\n", " \\bar{x}_\\mathcal{B} - \\bar{t} \\Delta x_\\mathcal{B} =\n", " \\begin{bmatrix} 1\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2 \\end{bmatrix} -\n", " (0) \\begin{bmatrix}\n", " -1\\\\ -2\\\\ 1\\\\ -1\\\\ 1\\\\ -1\\\\ 1\\\\ -1\n", " \\end{bmatrix} =\n", " \\begin{bmatrix} 1\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2 \\end{bmatrix}\n", " \\quad & \\quad\n", " \\bar{x}_j &\\leftarrow \\bar{t} = 0\n", "\n", " and\n", "\n", " .. math::\n", "\n", " z^*_\\mathcal{N} &\\leftarrow z^*_\\mathcal{N} - s \\Delta z_\\mathcal{N} =\n", " \\begin{bmatrix} -n\\mu' \\\\ 1 + n\\mu'\\\\ 1\\\\ 1\\\\ 1 \\end{bmatrix} -\n", " (-1 - n\\mu') \\begin{bmatrix} 1\\\\ -1\\\\ 1\\\\ 0\\\\ 0 \\end{bmatrix} =\n", " \\begin{bmatrix} 1\\\\ 0\\\\ 2 + n\\mu'\\\\ 0\\\\ 0 \\end{bmatrix}\n", " \\quad & \\quad\n", " z^*_i &\\leftarrow s = -1 - n\\mu'\\\\\n", " \\bar{z}_\\mathcal{N} &\\leftarrow\n", " \\bar{z}_\\mathcal{N} - \\bar{s} \\Delta z_\\mathcal{N} =\n", " \\begin{bmatrix} -1\\\\ 2\\\\ 1\\\\ 1\\\\ 1 \\end{bmatrix} -\n", " (-2) \\begin{bmatrix} 1\\\\ -1\\\\ 1\\\\ 0\\\\ 0 \\end{bmatrix} =\n", " \\begin{bmatrix} 1\\\\ 0\\\\ 3\\\\ 1\\\\ 1 \\end{bmatrix}\n", " \\quad & \\quad\n", " \\bar{z}_i &\\leftarrow \\bar{s} = -2.\n", "\n", "Step 5\n", " The current basis is updated to\n", "\n", " .. math::\n", "\n", " \\mathcal{B} &= \\{ 0, 6, 1, 8, 9, 10, 11, 12 \\} \\quad & \\quad\n", " \\mathcal{N} &= \\{ 5, 7, 2, 3, 4 \\}\\\\\n", " B &= \\begin{bmatrix}\n", " 1 & & -1\\\\\n", " -1 & 1 & -1\\\\\n", " 1\\\\\n", " -1 & & & 1\\\\\n", " 1 & & & & 1\\\\\n", " -1 & & & & & 1\\\\\n", " 1 & & & & & & 1\\\\\n", " -1 & & & & & & & 1\n", " \\end{bmatrix} \\quad & \\quad\n", " N &= \\begin{bmatrix}\n", " 1\\\\\n", " \\\\\n", " & 1 & -1\\\\\n", " & & -1\\\\\n", " & & & -1\\\\\n", " & & & -1\\\\\n", " & & & & -1\\\\\n", " & & & & -1\n", " \\end{bmatrix}\\\\\n", " x^*_\\mathcal{B} &= \\begin{bmatrix}\n", " 2\\\\ 2\\\\ 1\\\\ 0\\\\ 2\\\\ -2\\\\ 6\\\\ -6\n", " \\end{bmatrix}\n", " \\quad & \\quad\n", " z^*_\\mathcal{N} &= \\begin{bmatrix}\n", " 1\\\\ -1 - n\\mu'\\\\ 2 + n\\mu'\\\\ 0\\\\ 0\n", " \\end{bmatrix}\\\\\n", " \\bar{x}_\\mathcal{B} &= \\begin{bmatrix}\n", " 1\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2\n", " \\end{bmatrix}\n", " \\quad & \\quad\n", " \\bar{z}_\\mathcal{N} &= \\begin{bmatrix} 1\\\\ -2\\\\ 3\\\\ 1\\\\ 1 \\end{bmatrix}.\n", "\n", " The current range for :math:`\\mu` is\n", "\n", " .. math::\n", "\n", " \\frac{-2 - n\\mu'}{3} \\leq \\mu \\leq \\frac{-1 - n\\mu'}{2}.\n", "\n", " By inspection, the dictionary is optimal as long as\n", " :math:`\\mu' \\leq -\\frac{2}{n}`.\n", "\n", "Iteration 3\n", "^^^^^^^^^^^\n", "\n", "Step 1\n", " Computing\n", "\n", " .. math::\n", "\n", " \\mu^*\n", " &= \\max\\left\\{\n", " \\max_{j \\in \\mathcal{N}, \\bar{z}_j > 0} -\\frac{z^*_j}{\\bar{z}_j},\n", " \\max_{i \\in \\mathcal{B}, \\bar{x}_i > 0} -\\frac{x^*_i}{\\bar{x}_i},\n", " \\right\\}\\\\\n", " &= \\max\\left\\{\n", " \\max\\left\\{\n", " -\\frac{1}{1}, -\\frac{2 + n\\mu'}{3}, -\\frac{0}{1}, -\\frac{0}{1}\n", " \\right\\},\n", " \\max\\left\\{\n", " -\\frac{2}{1}, -\\frac{2}{2},\n", " -\\frac{0}{2}, -\\frac{-2}{2}, -\\frac{-6}{2}\n", " \\right\\}\n", " \\right\\}\\\\\n", " &= -\\frac{1 + n\\mu'}{2}\n", "\n", " shows that :math:`x_2` is the entering variable for\n", " :math:`\\mu' < -\\frac{2}{n}`.\n", "\n", "Step 2\n", " The maximum is achieved by :math:`j = 2 \\in \\mathcal{N}`, so do one step of\n", " the primal simplex method.\n", "\n", " The primal step direction is\n", "\n", " .. math::\n", "\n", " \\Delta x_\\mathcal{B} =\n", " B^{-1} N e_j =\n", " \\begin{bmatrix}\n", " & 1 & -1\\\\\n", " -1 & 2 & -2\\\\\n", " -1 & 1 & -1\\\\\n", " & 1 & -2\\\\\n", " & -1 & 1 & -1\\\\\n", " & 1 & -1 & -1\\\\\n", " & -1 & 1 & & -1\\\\\n", " & 1 & -1 & & -1\n", " \\end{bmatrix} \\begin{bmatrix} 0\\\\ 0\\\\ 1\\\\ 0\\\\ 0 \\end{bmatrix} =\n", " \\begin{bmatrix} -1\\\\ -2\\\\ -1\\\\ -2\\\\ 1\\\\ -1\\\\ 1\\\\ -1 \\end{bmatrix}.\n", "\n", " The index of the leaving variable is\n", "\n", " .. math::\n", "\n", " \\argmin_{i \\in \\mathcal{B}}\n", " \\left\\{\n", " \\frac{x^*_i + \\mu^* \\bar{x}_i}{\\Delta x_i} \\colon \\Delta x_i > 0\n", " \\right\\}\n", " &= \\argmin_{i \\in \\mathcal{B}}\\left\\{\n", " \\frac{2 + \\mu^* (0)}{1},\n", " \\frac{6 + \\mu^* (0)}{1}\n", " \\right\\}\\\\\n", " &= 9.\n", "\n", " The dual step direction is\n", "\n", " .. math::\n", "\n", " \\Delta z_\\mathcal{N} =\n", " -\\left( B^{-1} N \\right)^\\top e_i =\n", " -\\begin{bmatrix}\n", " & -1 & -1\\\\\n", " 1 & 2 & 1 & 1 & -1 & 1 & -1 & 1\\\\\n", " -1 & -2 & -1 & -2 & 1 & -1 & 1 & -1\\\\\n", " & & & & -1 & -1\\\\\n", " & & & & & & -1 & -1\n", " \\end{bmatrix} \\begin{bmatrix}\n", " 0\\\\ 0\\\\ 0\\\\ 0\\\\ 1\\\\ 0\\\\ 0\\\\ 0\n", " \\end{bmatrix} =\n", " \\begin{bmatrix} 0\\\\ 1\\\\ -1\\\\ 1\\\\ 0 \\end{bmatrix}.\n", "\n", "Step 3\n", " The step length adjustments are given by\n", "\n", " .. math::\n", "\n", " t &= \\frac{x^*_i}{\\Delta x_i} = \\frac{2}{1} = 2\n", " \\quad & \\quad\n", " \\bar{t} = \\frac{\\bar{x}_i}{\\Delta x_i} = \\frac{0}{1} = 0\\\\\n", " s &= \\frac{z^*_j}{\\Delta z_j} = \\frac{2 + n\\mu'}{-1} = -2 - n\\mu'\n", " \\quad & \\quad\n", " \\bar{s} = \\frac{\\bar{z}_j}{\\Delta z_j} = \\frac{3}{-1} = -3.\n", "\n", "Step 4\n", " The current primal and dual solutions are updated to\n", "\n", " .. math::\n", "\n", " x^*_\\mathcal{B} &\\leftarrow x^*_\\mathcal{B} - t \\Delta x_\\mathcal{B} =\n", " \\begin{bmatrix} 2\\\\ 2\\\\ 1\\\\ 0\\\\ 2\\\\ -2\\\\ 6\\\\ -6 \\end{bmatrix} -\n", " (2) \\begin{bmatrix}\n", " -1\\\\ -2\\\\ -1\\\\ -2\\\\ 1\\\\ -1\\\\ 1\\\\ -1\n", " \\end{bmatrix} =\n", " \\begin{bmatrix} 4\\\\ 6\\\\ 3\\\\ 4\\\\ 0\\\\ 0\\\\ 4\\\\ -4 \\end{bmatrix}\n", " \\quad & \\quad\n", " x^*_j &\\leftarrow t = 2\\\\\n", " \\bar{x}_\\mathcal{B} &\\leftarrow\n", " \\bar{x}_\\mathcal{B} - \\bar{t} \\Delta x_\\mathcal{B} =\n", " \\begin{bmatrix} 1\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2 \\end{bmatrix} -\n", " (0) \\begin{bmatrix}\n", " -1\\\\ -2\\\\ -1\\\\ -2\\\\ 1\\\\ -1\\\\ 1\\\\ -1\n", " \\end{bmatrix} =\n", " \\begin{bmatrix} 1\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2 \\end{bmatrix}\n", " \\quad & \\quad\n", " \\bar{x}_j &\\leftarrow \\bar{t} = 0\n", "\n", " and\n", "\n", " .. math::\n", "\n", " z^*_\\mathcal{N} &\\leftarrow z^*_\\mathcal{N} - s \\Delta z_\\mathcal{N} =\n", " \\begin{bmatrix} 1\\\\ -1 - n\\mu'\\\\ 2 + n\\mu'\\\\ 0\\\\ 0 \\end{bmatrix} -\n", " (-2 - n\\mu') \\begin{bmatrix} 0\\\\ 1\\\\ -1\\\\ 1\\\\ 0 \\end{bmatrix} =\n", " \\begin{bmatrix} 1\\\\ 1\\\\ 0\\\\ 2 + n\\mu'\\\\ 0 \\end{bmatrix}\n", " \\quad & \\quad\n", " z^*_i &\\leftarrow s = -2 - n\\mu'\\\\\n", " \\bar{z}_\\mathcal{N} &\\leftarrow\n", " \\bar{z}_\\mathcal{N} - \\bar{s} \\Delta z_\\mathcal{N} =\n", " \\begin{bmatrix} 1\\\\ -2\\\\ 3\\\\ 1\\\\ 1 \\end{bmatrix} -\n", " (-3) \\begin{bmatrix} 0\\\\ 1\\\\ -1\\\\ 1\\\\ 0 \\end{bmatrix} =\n", " \\begin{bmatrix} 1\\\\ 1\\\\ 0\\\\ 4\\\\ 1 \\end{bmatrix}\n", " \\quad & \\quad\n", " \\bar{z}_i &\\leftarrow \\bar{s} = -3.\n", "\n", "Step 5\n", " The current basis is updated to\n", "\n", " .. math::\n", "\n", " \\mathcal{B} &= \\{ 0, 6, 1, 8, 2, 10, 11, 12 \\} \\quad & \\quad\n", " \\mathcal{N} &= \\{ 5, 7, 9, 3, 4 \\}\\\\\n", " B &= \\begin{bmatrix}\n", " 1 & & -1\\\\\n", " -1 & 1 & -1\\\\\n", " 1 & & & & -1\\\\\n", " -1 & & & 1 & -1\\\\\n", " 1\\\\\n", " -1 & & & & & 1\\\\\n", " 1 & & & & & & 1\\\\\n", " -1 & & & & & & & 1\n", " \\end{bmatrix} \\quad & \\quad\n", " N &= \\begin{bmatrix}\n", " 1\\\\\n", " \\\\\n", " & 1\\\\\n", " \\\\\n", " & & 1 & -1\\\\\n", " & & & -1\\\\\n", " & & & & -1\\\\\n", " & & & & -1\n", " \\end{bmatrix}\\\\\n", " x^*_\\mathcal{B} &= \\begin{bmatrix}\n", " 4\\\\ 6\\\\ 3\\\\ 4\\\\ 2\\\\ 0\\\\ 4\\\\ -4\n", " \\end{bmatrix}\n", " \\quad & \\quad\n", " z^*_\\mathcal{N} &= \\begin{bmatrix}\n", " 1\\\\ 1\\\\ -2 - n\\mu'\\\\ 2 + n\\mu'\\\\ 0\n", " \\end{bmatrix}\\\\\n", " \\bar{x}_\\mathcal{B} &= \\begin{bmatrix}\n", " 1\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2\n", " \\end{bmatrix}\n", " \\quad & \\quad\n", " \\bar{z}_\\mathcal{N} &= \\begin{bmatrix} 1\\\\ 1\\\\ -3\\\\ 4\\\\ 1 \\end{bmatrix}.\n", "\n", " The current range for :math:`\\mu` is\n", "\n", " .. math::\n", "\n", " \\frac{-2 - n\\mu'}{4} \\leq \\mu \\leq \\frac{-2 - n\\mu'}{3}.\n", "\n", " By inspection, the dictionary is optimal as long as\n", " :math:`\\mu' \\leq -\\frac{10}{n}`.\n", "\n", "Iteration 4\n", "^^^^^^^^^^^\n", "\n", "Step 1\n", " Computing\n", "\n", " .. math::\n", "\n", " \\mu^*\n", " &= \\max\\left\\{\n", " \\max_{j \\in \\mathcal{N}, \\bar{z}_j > 0} -\\frac{z^*_j}{\\bar{z}_j},\n", " \\max_{i \\in \\mathcal{B}, \\bar{x}_i > 0} -\\frac{x^*_i}{\\bar{x}_i},\n", " \\right\\}\\\\\n", " &= \\max\\left\\{\n", " \\max\\left\\{\n", " -\\frac{1}{1}, -\\frac{1}{1}, -\\frac{2 + n\\mu'}{4}, -\\frac{0}{1}\n", " \\right\\},\n", " \\max\\left\\{\n", " -\\frac{4}{1}, -\\frac{6}{2},\n", " -\\frac{4}{2}, -\\frac{0}{2}, -\\frac{-4}{2}\n", " \\right\\}\n", " \\right\\}\\\\\n", " &= -\\frac{2 + n\\mu'}{4}\n", "\n", " shows that :math:`x_3` is the entering variable for\n", " :math:`\\mu' < -\\frac{10}{n}`.\n", "\n", "Step 2\n", " The maximum is achieved by :math:`j = 3 \\in \\mathcal{N}`, so do one step of\n", " the primal simplex method.\n", "\n", " The primal step direction is\n", "\n", " .. math::\n", "\n", " \\Delta x_\\mathcal{B} =\n", " B^{-1} N e_j =\n", " \\begin{bmatrix}\n", " & & 1 & -1\\\\\n", " -1 & & 2 & -2\\\\\n", " -1 & & 1 & -1\\\\\n", " & -1 & 2 & -2\\\\\n", " & -1 & 1 & -1\\\\\n", " & & 1 & -2\\\\\n", " & & -1 & 1 & -1\\\\\n", " & & 1 & -1 & -1\n", " \\end{bmatrix} \\begin{bmatrix} 0\\\\ 0\\\\ 0\\\\ 1\\\\ 0 \\end{bmatrix} =\n", " \\begin{bmatrix} -1\\\\ -2\\\\ -1\\\\ -2\\\\ -1\\\\ -2\\\\ 1\\\\ -1 \\end{bmatrix}.\n", "\n", " The index of the leaving variable is\n", "\n", " .. math::\n", "\n", " \\argmin_{i \\in \\mathcal{B}}\n", " \\left\\{\n", " \\frac{x^*_i + \\mu^* \\bar{x}_i}{\\Delta x_i} \\colon \\Delta x_i > 0\n", " \\right\\}\n", " &= \\argmin_{i \\in \\mathcal{B}}\\left\\{\n", " \\frac{4 + \\mu^* (0)}{1}\n", " \\right\\}\\\\\n", " &= 11.\n", "\n", " The dual step direction is\n", "\n", " .. math::\n", "\n", " \\Delta z_\\mathcal{N} =\n", " -\\left( B^{-1} N \\right)^\\top e_i =\n", " -\\begin{bmatrix}\n", " & -1 & -1\\\\\n", " & & & -1 & -1\\\\\n", " 1 & 2 & 1 & 2 & 1 & 1 & -1 & 1\\\\\n", " -1 & -2 & -1 & -2 & -1 & -2 & 1 & -1\\\\\n", " & & & & & & -1 & -1\n", " \\end{bmatrix} \\begin{bmatrix}\n", " 0\\\\ 0\\\\ 0\\\\ 0\\\\ 0\\\\ 0\\\\ 1\\\\ 0\n", " \\end{bmatrix} =\n", " \\begin{bmatrix} 0\\\\ 0\\\\ 1\\\\ -1\\\\ 1 \\end{bmatrix}.\n", "\n", "Step 3\n", " The step length adjustments are given by\n", "\n", " .. math::\n", "\n", " t &= \\frac{x^*_i}{\\Delta x_i} = \\frac{4}{1} = 4\n", " \\quad & \\quad\n", " \\bar{t} = \\frac{\\bar{x}_i}{\\Delta x_i} = \\frac{0}{1} = 0\\\\\n", " s &= \\frac{z^*_j}{\\Delta z_j} = \\frac{2 + n\\mu'}{-1} = -2 - n\\mu'\n", " \\quad & \\quad\n", " \\bar{s} = \\frac{\\bar{z}_j}{\\Delta z_j} = \\frac{4}{-1} = -4.\n", "\n", "Step 4\n", " The current primal and dual solutions are updated to\n", "\n", " .. math::\n", "\n", " x^*_\\mathcal{B} &\\leftarrow x^*_\\mathcal{B} - t \\Delta x_\\mathcal{B} =\n", " \\begin{bmatrix} 4\\\\ 6\\\\ 3\\\\ 4\\\\ 2\\\\ 0\\\\ 4\\\\ -4 \\end{bmatrix} -\n", " (4) \\begin{bmatrix}\n", " -1\\\\ -2\\\\ -1\\\\ -2\\\\ -1\\\\ -2\\\\ 1\\\\ -1\n", " \\end{bmatrix} =\n", " \\begin{bmatrix} 8\\\\ 14\\\\ 7\\\\ 12\\\\ 6\\\\ 8\\\\ 0\\\\ 0 \\end{bmatrix}\n", " \\quad & \\quad\n", " x^*_j &\\leftarrow t = 4\\\\\n", " \\bar{x}_\\mathcal{B} &\\leftarrow\n", " \\bar{x}_\\mathcal{B} - \\bar{t} \\Delta x_\\mathcal{B} =\n", " \\begin{bmatrix} 1\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2 \\end{bmatrix} -\n", " (0) \\begin{bmatrix}\n", " -1\\\\ -2\\\\ -1\\\\ -2\\\\ -1\\\\ -2\\\\ 1\\\\ -1\n", " \\end{bmatrix} =\n", " \\begin{bmatrix} 1\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2 \\end{bmatrix}\n", " \\quad & \\quad\n", " \\bar{x}_j &\\leftarrow \\bar{t} = 0\n", "\n", " and\n", "\n", " .. math::\n", "\n", " z^*_\\mathcal{N} &\\leftarrow z^*_\\mathcal{N} - s \\Delta z_\\mathcal{N} =\n", " \\begin{bmatrix} 1\\\\ 1\\\\ -2 - n\\mu'\\\\ 2 + n\\mu'\\\\ 0 \\end{bmatrix} -\n", " (-2 - n\\mu') \\begin{bmatrix} 0\\\\ 0\\\\ 1\\\\ -1\\\\ 1 \\end{bmatrix} =\n", " \\begin{bmatrix} 1\\\\ 1\\\\ 0\\\\ 0\\\\ 2 + n\\mu' \\end{bmatrix}\n", " \\quad & \\quad\n", " z^*_i &\\leftarrow s = -2 - n\\mu'\\\\\n", " \\bar{z}_\\mathcal{N} &\\leftarrow\n", " \\bar{z}_\\mathcal{N} - \\bar{s} \\Delta z_\\mathcal{N} =\n", " \\begin{bmatrix} 1\\\\ 1\\\\ -3\\\\ 4\\\\ 1 \\end{bmatrix} -\n", " (-4) \\begin{bmatrix} 0\\\\ 0\\\\ 1\\\\ -1\\\\ 1 \\end{bmatrix} =\n", " \\begin{bmatrix} 1\\\\ 1\\\\ 1\\\\ 0\\\\ 5 \\end{bmatrix}\n", " \\quad & \\quad\n", " \\bar{z}_i &\\leftarrow \\bar{s} = -4.\n", "\n", "Step 5\n", " The current basis is updated to\n", "\n", " .. math::\n", "\n", " \\mathcal{B} &= \\{ 0, 6, 1, 8, 2, 10, 3, 12 \\} \\quad & \\quad\n", " \\mathcal{N} &= \\{ 5, 7, 9, 11, 4 \\}\\\\\n", " B &= \\begin{bmatrix}\n", " 1 & & -1\\\\\n", " -1 & 1 & -1\\\\\n", " 1 & & & & -1\\\\\n", " -1 & & & 1 & -1\\\\\n", " 1 & & & & & & -1\\\\\n", " -1 & & & & & 1 & -1\\\\\n", " 1\\\\\n", " -1 & & & & & & & 1\n", " \\end{bmatrix} \\quad & \\quad\n", " N &= \\begin{bmatrix}\n", " 1\\\\\n", " \\\\\n", " & 1\\\\\n", " \\\\\n", " & & 1\\\\\n", " \\\\\n", " & & & 1 & -1\\\\\n", " & & & & -1\n", " \\end{bmatrix}\\\\\n", " x^*_\\mathcal{B} &= \\begin{bmatrix}\n", " 8\\\\ 14\\\\ 7\\\\ 12\\\\ 6\\\\ 8\\\\ 4\\\\ 0\n", " \\end{bmatrix}\n", " \\quad & \\quad\n", " z^*_\\mathcal{N} &= \\begin{bmatrix}\n", " 1\\\\ 1\\\\ 0\\\\ -2 - n\\mu'\\\\ 2 + n\\mu'\n", " \\end{bmatrix}\\\\\n", " \\bar{x}_\\mathcal{B} &= \\begin{bmatrix}\n", " 1\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2\\\\ 0\\\\ 2\n", " \\end{bmatrix}\n", " \\quad & \\quad\n", " \\bar{z}_\\mathcal{N} &= \\begin{bmatrix} 1\\\\ 1\\\\ 1\\\\ -4\\\\ 5 \\end{bmatrix}.\n", "\n", " The current range for :math:`\\mu` is\n", "\n", " .. math::\n", "\n", " \\frac{-2 - n\\mu'}{5} \\leq \\mu \\leq \\frac{-2 - n\\mu'}{4}.\n", "\n", " By inspection, the dictionary is optimal as long as\n", " :math:`\\mu' \\leq -\\frac{2}{n}`.\n", "\n", "Iteration 5\n", "^^^^^^^^^^^\n", "\n", "Step 1\n", " Computing\n", "\n", " .. math::\n", "\n", " \\mu^*\n", " &= \\max\\left\\{\n", " \\max_{j \\in \\mathcal{N}, \\bar{z}_j > 0} -\\frac{z^*_j}{\\bar{z}_j},\n", " \\max_{i \\in \\mathcal{B}, \\bar{x}_i > 0} -\\frac{x^*_i}{\\bar{x}_i},\n", " \\right\\}\\\\\n", " &= \\max\\left\\{\n", " \\max\\left\\{\n", " -\\frac{1}{1}, -\\frac{1}{1}, -\\frac{0}{1}, -\\frac{2 + n\\mu'}{5},\n", " \\right\\},\n", " \\max\\left\\{\n", " -\\frac{8}{1}, -\\frac{14}{2},\n", " -\\frac{12}{2}, -\\frac{8}{2}, -\\frac{0}{2}\n", " \\right\\}\n", " \\right\\}\\\\\n", " &= -\\frac{2 + n\\mu'}{5}\n", "\n", " shows that :math:`x_4` is the entering variable for\n", " :math:`\\mu' < -\\frac{2}{n}`.\n", "\n", "Step 2\n", " The maximum is achieved by :math:`j = 4 \\in \\mathcal{N}`, so do one step of\n", " the primal simplex method.\n", "\n", " The primal step direction is\n", "\n", " .. math::\n", "\n", " \\Delta x_\\mathcal{B} =\n", " B^{-1} N e_j =\n", " \\begin{bmatrix}\n", " & & & 1 & -1\\\\\n", " -1 & & & 2 & -2\\\\\n", " -1 & & & 1 & -1\\\\\n", " & -1 & & 2 & -2\\\\\n", " & -1 & & 1 & -1\\\\\n", " & & -1 & 2 & -2\\\\\n", " & & -1 & 1 & -1\\\\\n", " & & & 1 & -2\n", " \\end{bmatrix} \\begin{bmatrix} 0\\\\ 0\\\\ 0\\\\ 0\\\\ 1 \\end{bmatrix} =\n", " \\begin{bmatrix} -1\\\\ -2\\\\ -1\\\\ -2\\\\ -1\\\\ -2\\\\ -1\\\\ -2 \\end{bmatrix}.\n", "\n", " The index of the leaving variable is\n", "\n", " .. math::\n", "\n", " \\argmin_{i \\in \\mathcal{B}}\n", " \\left\\{\n", " \\frac{x^*_i + \\mu^* \\bar{x}_i}{\\Delta x_i} \\colon \\Delta x_i > 0\n", " \\right\\} =\n", " \\emptyset.\n", "\n", " Therefore, the primal is unbounded (i.e. the dual is infeasible) when\n", " :math:`\\mu' < -\\frac{2}{n}`.\n", "\n", "(d)\n", "---\n", "\n", "Even though there is a pattern for the :math:`k`-th iteration, it is unclear\n", "whether the results of (c) is correct because :math:`\\mu` and :math:`\\mu'` have\n", "strict boundary conditions." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy\n", "\n", "n = 10\n", "b = numpy.random.randint(-100, 100, size=n)\n", "mu = 0\n", "print(numpy.sort(b))\n", "from pymprog import model\n", "\n", "lp = model('L1-Regression LP')\n", "\n", "#lp.verbose(True)\n", "x = lp.var('x', bounds=(None, None))\n", "t = lp.var('t', n, bounds=(0, None))\n", "\n", "lp.minimize(sum(t_j + mu * (x - b_j) for b_j, t_j in zip(b, t)))\n", "\n", "for i in range(n):\n", " -t[i] <= x - b[i] <= t[i]\n", "\n", "lp.solve()\n", "lp.sensitivity()\n", "\n", "lp.end()" ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Exercise 12.11\n", "==============\n", "\n", "Already proved in :ref:`another exercise `." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ ".. rubric:: References\n", "\n", ".. bibliography:: chapter-12.bib" ] } ], "metadata": { "anaconda-cloud": {}, "celltoolbar": "Raw Cell Format", "kernelspec": { "display_name": "Python [default]", "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.5.2" } }, "nbformat": 4, "nbformat_minor": 1 }