tutorial_wLasso.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "This notebook is a tutorial on using VF-iDCA to solve the hyperparameter selection problem. And we take the following w-Lasso problem for example,\n",
    "$$\\begin{aligned}\n",
    "\\min_{\\lambda}\\ & \n",
    "\\|X_{\\text{val}}\\hat{w} - y_{\\text{val}}\\|^2 \\\\ \n",
    "\\text{where } & \n",
    "\\hat{w} \\in \\mathop{\\arg\\min}_{w} \n",
    "\\|X_{\\text{tr}}w - y_{\\text{tr}}\\|^2 + \\sum_{i=1}^p \\lambda_i |w_i|\n",
    "\\end{aligned}$$\n",
    "With hyperparameter decoupling, we are solving the following problem\n",
    "$$\\begin{aligned}\n",
    "& \\min_{\\lambda}\\ \n",
    "\\|X_{\\text{val}}\\hat{w} - y_{\\text{val}}\\|^2 \\\\ \n",
    "& \\begin{aligned}\n",
    "\\text{where } \\hat{w} \\in \\mathop{\\arg\\min}_{w}\\ &\n",
    "\\|X_{\\text{tr}}w - y_{\\text{tr}}\\|^2 \\\\\n",
    "\\text{s.t. } & |w_i| \\le r_i, \\quad i = 1,\\dots, p\n",
    "\\end{aligned}\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Structure\n",
    "Here is the structure of this tutorial:\n",
    "- How to construct and solve a subproblem with `cvxpy`\n",
    "- How to use VF-iDCA\n",
    "- Some advice for acceleration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving w-Lasso problem with `cvxpy`\n",
    "We firstly generate some synetic data for illustration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate Data\n",
    "import numpy as np\n",
    "\n",
    "np.random.seed(42) # for reproduction\n",
    "\n",
    "n, p = 600, 100\n",
    "beta_nonzero = np.array([1,2,3,4,5])\n",
    "beta_real = np.concatenate([beta_nonzero, np.zeros(p - len(beta_nonzero))])\n",
    "X = np.random.randn(n, p)\n",
    "y_true = X @ beta_real\n",
    "\n",
    "# add noise\n",
    "snr = 2\n",
    "epsilon = np.random.randn(n)\n",
    "SNR_factor = snr / np.linalg.norm(y_true) * np.linalg.norm(epsilon)\n",
    "y = y_true + 1.0 / SNR_factor * epsilon\n",
    "\n",
    "# split data\n",
    "nTr, nVal = 100, 100\n",
    "XTr, XVal, XTest = X[0:nTr], X[nTr:nTr + nVal], X[nTr + nVal:]\n",
    "yTr, yVal, yTest = y[0:nTr], y[nTr:nTr + nVal], y[nTr + nVal:]\n",
    "\n",
    "## define err function\n",
    "def ls_err(X, y, w):\n",
    "    return np.sum(np.square(X @ w - y))/len(y)\n",
    "\n",
    "def err_report(w):\n",
    "    print(\"%20s %8.3f\" % (\"train err:\", ls_err(XTr, yTr, w)))\n",
    "    print(\"%20s %8.3f\" % (\"validate err:\", ls_err(XVal, yVal, w)))\n",
    "    print(\"%20s %8.3f\" % (\"test err:\", ls_err(XTest, yTest, w)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The original wlasso problem could be solving as "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "optimal objective is 24.922\n",
      "first 3 values of optimal w are [6.09002687e-08 1.87646683e+00 2.46803579e+00]\n",
      "          train err:   10.772\n",
      "       validate err:   12.111\n",
      "           test err:   15.348\n"
     ]
    }
   ],
   "source": [
    "import cvxpy as cp \n",
    "\n",
    "## define the wlasso problem\n",
    "w = cp.Variable(p) # define the variables\n",
    "lam = cp.Parameter(p, nonneg=True) # define the regularization parameters, nonneg=True means its value will not be negtive\n",
    "wlasso_obj = cp.sum_squares(XTr @ w - yTr) / len(yTr) + lam @ cp.abs(w) # define the objective of w-lasso problem, to check the supported function, one should refer to https://www.cvxpy.org/tutorial/functions/index.html#functions\n",
    "wlasso = cp.Problem(cp.Minimize(wlasso_obj)) # define the problem\n",
    "\n",
    "## solve a wlasso problem\n",
    "lam.value = np.ones(p) # clarify the value of parameters, here we use [1,1,..,1]\n",
    "wlasso.solve(solver = cp.ECOS) # use ECOS solver to solve the problem\n",
    "print(\"optimal objective is %.3f\" % (wlasso.value))\n",
    "print(\"first 3 values of optimal w are\", w.value[:3])\n",
    "err_report(w.value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The decoupled wlasso problem could be solving as "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "optimal objective is 8.203\n",
      "first 3 values of optimal w are [0.70813074 0.99999992 0.99999999]\n",
      "first value of the corresponding lam are 6.019831201251446e-08\n",
      "          train err:    8.203\n",
      "       validate err:  107.119\n",
      "           test err:  135.303\n"
     ]
    }
   ],
   "source": [
    "## define the wlasso problem\n",
    "w = cp.Variable(p) # define the variables\n",
    "r = cp.Parameter(p, nonneg=True) # define the regularization parameters\n",
    "wlasso_obj = cp.sum_squares(XTr @ w - yTr) / len(yTr)\n",
    "cons = [w[i] <= r[i] for i in range(p)]\n",
    "wlasso = cp.Problem(cp.Minimize(wlasso_obj), cons) # define the problem\n",
    "\n",
    "## solve a wlasso problem\n",
    "r.value = np.ones(p) # clarify the value of parameters, here we use [1,1,..,1]\n",
    "wlasso.solve(solver = cp.ECOS) # use ECOS solver to solve the problem\n",
    "print(\"optimal objective is %.3f\" % (wlasso.value))\n",
    "print(\"first 3 values of optimal w are\", w.value[:3])\n",
    "print(\"first value of the corresponding lam are\", cons[0].dual_value)\n",
    "err_report(w.value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use VF-iDCA to select hyperparamter\n",
    "Before we apply the VF-iDCA method, we would like to clarify how to solve the following approximated problem\n",
    "$$\\begin{aligned}\n",
    "\\min_{w, r} \\phi_k(w, r) = & \\|X_{\\text{Val}}w - y\\|^2 + \\frac{\\rho}2 (\\|w - w^k\\|^2 + \\|r - r^k\\|^2)\\\\  & + \\alpha_k \\max\\{0,V_k(w, r), \\max_{i=1,\\dots,p}|w_i| - r_i\\}\n",
    "\\end{aligned}$$\n",
    "where\n",
    "$$\\begin{aligned}\n",
    "V_k(w, r) = & \\|X_{\\text{Tr}}w - y\\|^2 - \\|X_{\\text{Tr}}\\hat{w}^k - y\\|^2\n",
    "+ \\langle \\gamma^k, r - r^k \\rangle \\\\\n",
    "= & \\|X_{\\text{Tr}}w - y\\|^2 - l(\\hat{w}^k)\n",
    "+ \\langle \\gamma^k, r - r^k \\rangle\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "d:\\_Tools\\anaconda3\\envs\\cvxpy\\lib\\site-packages\\cvxpy\\reductions\\solvers\\solving_chain.py:167: UserWarning: You are solving a parameterized problem that is not DPP. Because the problem is not DPP, subsequent solves will not be faster than the first one. For more information, see the documentation on Discplined Parametrized Programming, at\n",
      "\thttps://www.cvxpy.org/tutorial/advanced/index.html#disciplined-parametrized-programming\n",
      "  warnings.warn(dpp_error_msg)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "optimal objective is 26.108\n",
      "first 3 values of optimal w are [0.23810928 1.24042975 2.41663681]\n",
      "first 3 values of optimal r are [0.99999834 0.99999825 0.99999765]\n",
      "          train err:   10.654\n",
      "       validate err:    6.942\n",
      "           test err:   23.501\n"
     ]
    }
   ],
   "source": [
    "## define the approximated problem\n",
    "rho = 1. \n",
    "wU = cp.Variable(p) # define the variables\n",
    "rU = cp.Variable(p) # define the regularization parameters\n",
    "fL = cp.Parameter()\n",
    "wk = cp.Parameter(p)\n",
    "rk = cp.Parameter(p, nonneg=True)\n",
    "gamma = cp.Parameter(p)\n",
    "alpha = cp.Parameter(nonneg=True)\n",
    "\n",
    "ls_val = cp.sum_squares(XVal @ wU - yVal) / len(yVal)\n",
    "prox = cp.sum_squares(wU - wk) + cp.sum_squares(rU - rk)\n",
    "Vk = cp.sum_squares(XTr @ wU - yTr) / len(yTr) - fL + gamma @ (r - rk)\n",
    "violation = cp.maximum(*([0, Vk] + [cp.abs(wU[i]) - rU[i] for i in range(p)]))\n",
    "\n",
    "phi_k = ls_val + rho/2 * prox + alpha * violation\n",
    "bi_cons = [rU >= 0]\n",
    "\n",
    "dc_app = cp.Problem(cp.Minimize(phi_k), bi_cons)\n",
    "\n",
    "fL.value = wlasso.value\n",
    "wk.value = np.zeros(p)\n",
    "rk.value = np.ones(p)\n",
    "gamma.value = np.array([float(cons[i].dual_value) for i in range(p)])\n",
    "alpha.value = 1.\n",
    "\n",
    "dc_app.solve(solver = cp.ECOS)\n",
    "print(\"optimal objective is %.3f\" % (dc_app.value))\n",
    "print(\"first 3 values of optimal w are\", wU.value[:3])\n",
    "print(\"first 3 values of optimal r are\", rU.value[:3])\n",
    "err_report(wU.value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "820.303557332848"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fL.value "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To use VF-iDCA more easily, we need define both lower-level problem and approximated problem in the following way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8.20303667928094\n",
      "24.644251976385892\n"
     ]
    }
   ],
   "source": [
    "class DC_lower():\n",
    "    def __init__(self):\n",
    "        self.wL = cp.Variable(p)\n",
    "        self.rL = cp.Parameter(p, nonneg=True)\n",
    "        LSTr = cp.sum_squares(XTr @ self.wL - yTr) / len(yTr)\n",
    "        self.cons = [self.wL[i] <= self.rL[i] for i in range(p)]\n",
    "        self.dc_lower = cp.Problem(cp.Minimize(LSTr), self.cons)\n",
    "\n",
    "    def solve(self, r):\n",
    "        self.rL.value = r\n",
    "        result = self.dc_lower.solve(solver = cp.ECOS)\n",
    "        return result\n",
    "\n",
    "    def dual_value(self):\n",
    "        return np.array([float(self.cons[i].dual_value) for i in range(p)])\n",
    "\n",
    "class DC_approximated():\n",
    "    def __init__(self):\n",
    "        self.rho = 1.\n",
    "        self.delta = 5.\n",
    "        self.c_alpha = 1. \n",
    "\n",
    "        self.wU = cp.Variable(p)\n",
    "        self.rU = cp.Variable(p)\n",
    "        self.fL = cp.Parameter()\n",
    "        self.wk = cp.Parameter(p)\n",
    "        self.rk = cp.Parameter(p, nonneg=True)\n",
    "        self.gamma = cp.Parameter(p)\n",
    "        self.alpha = cp.Parameter(nonneg=True)\n",
    "\n",
    "        self.alpha.value = 1.\n",
    "\n",
    "        LSVal = cp.sum_squares(XVal @ self.wU - yVal) / len(yVal)\n",
    "        prox = cp.sum_squares(self.wU - self.wk) + cp.sum_squares(self.rU - self.rk)\n",
    "        Vk = cp.sum_squares(XTr @ self.wU - yTr) / len(yTr) - self.fL + self.gamma @ (self.rU - self.rk)\n",
    "        self.violation = cp.maximum(*([0, Vk] + [cp.abs(self.wU[i]) - self.rU[i] for i in range(p)]))\n",
    "\n",
    "        phi_k = LSVal + self.rho/2 * prox + self.alpha * self.violation\n",
    "        bi_cons = [self.rU >= 0]\n",
    "\n",
    "        self.dc_app = cp.Problem(cp.Minimize(phi_k), bi_cons)\n",
    "    \n",
    "    def solve(self):\n",
    "        self.dc_app.solve(solver = cp.ECOS)\n",
    "        return self.dc_app.value, self.wU.value, np.maximum(0, self.rU.value)\n",
    "    \n",
    "    def clare_variable(self, w, r):\n",
    "        self.wk.value = w \n",
    "        self.rk.value = r\n",
    "    \n",
    "    def clare_V(self, fL, gamma):\n",
    "        self.fL.value = fL\n",
    "        self.gamma.value = gamma\n",
    "\n",
    "# Test the function\n",
    "dclower = DC_lower()\n",
    "dcapp = DC_approximated()\n",
    "print(dclower.solve(np.ones(p)))\n",
    "fL = dclower.solve(np.ones(p))\n",
    "gamma = dclower.dual_value()\n",
    "dcapp.clare_variable(np.zeros(p), np.ones(p))\n",
    "dcapp.clare_V(fL, gamma)\n",
    "print(dcapp.solve()[0])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def iteration_err(w, r, wp, rp):\n",
    "    return np.sqrt(\n",
    "        np.sum(np.square(w - wp)) + np.sum(np.square(r - rp))\n",
    "    ) / np.sqrt(\n",
    "        1 + np.sum(np.square(w)) + np.sum(np.square(r))\n",
    "    )\n",
    "\n",
    "def VF_iDCA(dclower, dcapp, DC_Setting = dict()):\n",
    "    MAX_ITERATION = DC_Setting[\"MAX_ITERATION\"] if \"MAX_ITERATION\" in DC_Setting.keys() else 10\n",
    "    TOL = DC_Setting[\"TOL\"] if \"TOL\" in DC_Setting.keys() else 1e-2\n",
    "    dcapp.alpha.value = DC_Setting[\"alpha\"] if \"alpha\" in DC_Setting.keys() else 1.\n",
    "    dcapp.c_alpha = DC_Setting[\"c_alpha\"] if \"c_alpha\" in DC_Setting.keys() else 1.\n",
    "    dcapp.delta = DC_Setting[\"delta\"] if \"delta\" in DC_Setting.keys() else 5. \n",
    "    r = DC_Setting[\"initial_r\"] if \"initial_r\" in DC_Setting.keys() else np.ones(dcapp.rU.shape)\n",
    "    w = DC_Setting[\"initial_w\"] if \"initial_w\" in DC_Setting.keys() else np.zeros(dcapp.wU.shape)\n",
    "\n",
    "    for _ in range(MAX_ITERATION):\n",
    "        fL = dclower.solve(r)\n",
    "        gamma = dclower.dual_value()\n",
    "        dcapp.clare_variable(w, r)\n",
    "        dcapp.clare_V(fL, gamma)\n",
    "        _, wp, rp = dcapp.solve()\n",
    "\n",
    "        err = iteration_err(w, r, wp, rp) \n",
    "        penalty = dcapp.violation.value\n",
    "\n",
    "        if err < TOL and penalty < TOL:\n",
    "            break\n",
    "\n",
    "        if err * dcapp.alpha.value <= dcapp.c_alpha * min( 1., dcapp.alpha.value * dcapp.violation.value ):\n",
    "            dcapp.alpha.value = dcapp.alpha.value + dcapp.delta\n",
    "\n",
    "        w, r = wp, rp\n",
    "\n",
    "    return wp, rp  \n",
    "        \n",
    "dclower = DC_lower()\n",
    "dcapp = DC_approximated()\n",
    "wp, rp = VF_iDCA(dclower, dcapp, dict())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          train err:    5.856\n",
      "       validate err:    6.405\n",
      "           test err:   22.499\n"
     ]
    }
   ],
   "source": [
    "err_report(wp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to accelerate VF-iDCA\n",
    "However, such a programming works slowly since the parameterized problem is not DPP. \n",
    "\n",
    "Because the problem is not DPP, subsequent solves will not be faster than the first one. \n",
    "\n",
    "Documentation about DPP is at \n",
    "https://www.cvxpy.org/tutorial/advanced/index.html#disciplined-parametrized-programming\n",
    "\n",
    "The main problem is that DPP don't allow the multiplication between parameters (even when the signs are given), so you need to take good care on `alpha*violation`.\n",
    "\n",
    "For more details on acceleration, please refer to the documentation and the example codes of other problems.\n",
    "\n",
    "Any question about programming, please contact 11930905@mail.sustech.edu.cn."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.6 ('cvxpy')",
   "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.6"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "91f36096c2d25ba640567d5993293a7d6f8015c5f45ad7887fa4b03526a8efc2"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}