{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Method of multipliers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The method of multipliers is an algorithm for solving convex optimization problems. \n",
    "Suppose we have a problem of the form\n",
    "\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize} & f(x)\\\\\n",
    "\\mbox{subject to} & Ax = b,\n",
    "\\end{array}\n",
    "\n",
    "where $f$ is convex, $x \\in \\mathbf{R}^n$ is the optimization variable, and $A \\in \\mathbf{R}^{m \\times n}$ and $b \\in \\mathbf{R}^m$ are problem data.\n",
    "\n",
    "To apply the method of multipliers, we first form the augmented Lagrangian \n",
    "\n",
    "$$L_{\\rho}(x,y) = f(x) + y^T(Ax - b) + (\\rho/2)\\|Ax-b\\|^2_2.$$\n",
    "\n",
    "The dual function associated with the augmented Lagrangian is $g_{\\rho}(y) = \\inf_x L_{\\rho}(x,y)$. \n",
    "The dual function $g_{\\rho}(y)$ is concave and its maximal value is the same as the optimal value of the original problem.\n",
    "\n",
    "We maximize the dual function using gradient ascent. Each step of gradient ascent reduces to the $x$ and $y$ updates\n",
    "\n",
    "\\begin{array}{lll}\n",
    "x^{k+1} & := & \\mathop{\\rm argmin}_{x}\\left(f(x) + (y^k)^T(Ax - b) + (\\rho/2)\\left\\|Ax-b\\right\\|^2_2 \\right) \\\\\n",
    "y^{k+1} & := & y^{k} + \\rho(Ax^{k+1}-b)\n",
    "\\end{array}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following CVXPY script implements the method of multipliers and uses it to solve an optimization problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal value from CVXPY: 5.5905035557463005\n",
      "Optimal value from method of multipliers: 5.572761551213633\n"
     ]
    }
   ],
   "source": [
    "import cvxpy as cp\n",
    "import numpy as np\n",
    "np.random.seed(1)\n",
    "\n",
    "# Initialize data.\n",
    "MAX_ITERS = 10\n",
    "rho = 1.0\n",
    "n = 20\n",
    "m = 10\n",
    "A = np.random.randn(m,n)\n",
    "b = np.random.randn(m)\n",
    "\n",
    "# Initialize problem.\n",
    "x = cp.Variable(shape=n)\n",
    "f = cp.norm(x, 1)\n",
    "\n",
    "# Solve with CVXPY.\n",
    "cp.Problem(cp.Minimize(f), [A*x == b]).solve()\n",
    "print(\"Optimal value from CVXPY: {}\".format(f.value))\n",
    "\n",
    "# Solve with method of multipliers.\n",
    "resid = A*x - b\n",
    "y = cp.Parameter(shape=(m)); y.value = np.zeros(m)\n",
    "aug_lagr = f + y.T*resid + (rho/2)*cp.sum_squares(resid)\n",
    "for t in range(MAX_ITERS):\n",
    "    cp.Problem(cp.Minimize(aug_lagr)).solve()\n",
    "    y.value += rho*resid.value\n",
    "    \n",
    "print(\"Optimal value from method of multipliers: {}\".format(f.value))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
