group_lasso_and_sparse_group_lasso.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "female-verse",
   "metadata": {},
   "source": [
    "# Group Lasso\n",
    "Reference: Yuan, Ming, and Yi Lin. \"Model selection and estimation in regression with grouped variables.\" Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68.1 (2006): 49-67."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "exclusive-christianity",
   "metadata": {},
   "source": [
    "The objective function is $$min_{\\beta} \\frac12 \\| y-\\Sigma_{j=1}^JX_j\\beta_j \\|_2^2 + \\lambda \\| \\beta_j \\|_{K_j},$$\n",
    "where $\\| \\eta \\|_{K} = (\\eta^T K \\eta)^{1/2}$ for a vector $\\eta \\in R^d$ and a SPD matrix $K \\in R^{d\\times d}$,\n",
    "$\\beta_j$ is the coefficient vector of size $p_j$ and $p_j$ is the number of entries in $j$-th group.\n",
    "\n",
    "Let $K_j=p_jI_{p_j}$ then $\\| \\beta_j \\|_{K_j} = \\Sigma_{j=1}^J \\sqrt{p_j} \\| \\beta_j \\|_2$,\n",
    "**assume each $X_j$ is orthonormal as $X_j^T X_j = I_{p_j}$** for $j=1,...,J$,\n",
    "according to KKT condition the solution of $\\beta$  is\n",
    "$$ -X_j^T (y-\\Sigma_{j=1}^J X_j \\beta_j) + \\lambda \\sqrt{p_j} \\frac{\\beta_j}{\\|\\beta_j\\|_2} = 0, \\text{ if } \\beta_j \\neq 0$$\n",
    "otherwise, $\\beta_j = 0$, zero lies in the subgradient of L2 norm of $\\beta_j$ as $\\| \\frac{X_j^T \\left(y - X\\beta\\right)} {\\lambda \\sqrt{p_j}} \\|_2 \\leq 1$, then\n",
    "$$ \\| X_j^T \\left(y - X\\beta_{-j}\\right)\\|_2 \\leq \\lambda \\sqrt{p_j}, \\text{ if } \\beta_j = 0 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "mental-treasure",
   "metadata": {},
   "source": [
    "Denote $\\beta_{-j}=( \\beta_1^T, ..., \\beta_{j-1}^T, \\mathbf{0}, \\beta_{j+1}^T, \\beta_{J}^T)^T,\n",
    "r_j = y - X\\beta_{-j},\n",
    "s_j = X_j^T r_j$,\n",
    "recall $X_j^T X_j = I_{p_j}$,\n",
    "when $\\beta_j \\neq 0$,\n",
    "$$-s_j + X_j^T X_j\\beta_j + \\lambda \\sqrt{p_j} \\frac{\\beta_j}{\\|\\beta_j\\|_2} = 0$$\n",
    "$$s_j = (1 + \\frac{ \\lambda \\sqrt{p_j} } { \\|\\beta_j\\|_2 }) \\beta_j$$\n",
    "$s_j$ and $\\beta_j$ are parallel, thus $\\frac{\\beta_j}{\\|\\beta_j\\|_2} = \\frac{s_j}{\\|s_j\\|_2}$ and \n",
    "$\\beta_j = (1 - \\frac{ \\lambda \\sqrt{p_j} } { \\|s_j\\|_2 }) s_j$.\n",
    "\n",
    "Finally, combine solution of $\\beta_j$ in two cases,\n",
    "$$\\beta_j = (1 - \\frac{ \\lambda \\sqrt{p_j} } { \\|s_j\\|_2 })_+ s_j$$\n",
    "\n",
    "The algorithm is updated iteratively for $j=1,...,J$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "respective-wholesale",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "reliable-second",
   "metadata": {},
   "source": [
    "Simulation setting:\n",
    "$n=100, p=20$, $X$ is centered and scaled to variance 1 in each column, and X in each group holds orthonormality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "minute-screening",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "n = 100\n",
    "p = 20\n",
    "set.seed(0)\n",
    "# beta = rnorm(p, mean = 0, sd = 1)\n",
    "# print(mean(beta))\n",
    "# print(sd(beta))\n",
    "\n",
    "X = matrix( rnorm(n*p, mean = 0, sd = 1), n, p )\n",
    "# X = X - colMeans(X) # centering\n",
    "# Xsd <- sqrt(colMeans(X^2))\n",
    "# X = t( t(X)/Xsd ) /sqrt(p)\n",
    "X = scale(X, center=T, scale=T)\n",
    "\n",
    "## Note that if orthonormalize the whole matrix X, the update only happens in the first iteration, r_j would be updated since beta changes\n",
    "## but s_j=X_j^T(y-X\\beta_{-j}) remains the same since X_1^T X = [I, 0 , 0, 0] as an example then \\beta_j remains unchanged as well\n",
    "# X.svd = svd(X) \n",
    "# X = X.svd$u\n",
    "# # t(X) %*% X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "understood-solution",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "beta_true: 0.15 -0.33 0.25 -0.25 0.05 0 0 0 0.5 0.2 -0.25 0.12 -0.125 0 0 0 0 0 0 0 \n",
      "group_id: 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 \n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<ol>\n",
       "\t<li><style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>1</li><li>2</li><li>3</li><li>4</li><li>5</li></ol>\n",
       "</li>\n",
       "\t<li><style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>6</li><li>7</li><li>8</li><li>9</li><li>10</li></ol>\n",
       "</li>\n",
       "\t<li><style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>11</li><li>12</li><li>13</li><li>14</li><li>15</li></ol>\n",
       "</li>\n",
       "\t<li><style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>16</li><li>17</li><li>18</li><li>19</li><li>20</li></ol>\n",
       "</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate}\n",
       "\\item \\begin{enumerate*}\n",
       "\\item 1\n",
       "\\item 2\n",
       "\\item 3\n",
       "\\item 4\n",
       "\\item 5\n",
       "\\end{enumerate*}\n",
       "\n",
       "\\item \\begin{enumerate*}\n",
       "\\item 6\n",
       "\\item 7\n",
       "\\item 8\n",
       "\\item 9\n",
       "\\item 10\n",
       "\\end{enumerate*}\n",
       "\n",
       "\\item \\begin{enumerate*}\n",
       "\\item 11\n",
       "\\item 12\n",
       "\\item 13\n",
       "\\item 14\n",
       "\\item 15\n",
       "\\end{enumerate*}\n",
       "\n",
       "\\item \\begin{enumerate*}\n",
       "\\item 16\n",
       "\\item 17\n",
       "\\item 18\n",
       "\\item 19\n",
       "\\item 20\n",
       "\\end{enumerate*}\n",
       "\n",
       "\\end{enumerate}\n"
      ],
      "text/markdown": [
       "1. 1. 1\n",
       "2. 2\n",
       "3. 3\n",
       "4. 4\n",
       "5. 5\n",
       "\n",
       "\n",
       "\n",
       "2. 1. 6\n",
       "2. 7\n",
       "3. 8\n",
       "4. 9\n",
       "5. 10\n",
       "\n",
       "\n",
       "\n",
       "3. 1. 11\n",
       "2. 12\n",
       "3. 13\n",
       "4. 14\n",
       "5. 15\n",
       "\n",
       "\n",
       "\n",
       "4. 1. 16\n",
       "2. 17\n",
       "3. 18\n",
       "4. 19\n",
       "5. 20\n",
       "\n",
       "\n",
       "\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[[1]]\n",
       "[1] 1 2 3 4 5\n",
       "\n",
       "[[2]]\n",
       "[1]  6  7  8  9 10\n",
       "\n",
       "[[3]]\n",
       "[1] 11 12 13 14 15\n",
       "\n",
       "[[4]]\n",
       "[1] 16 17 18 19 20\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "group_number = 4\n",
    "card_per_group = p / group_number\n",
    "beta_true = c(0.15,-0.33,0.25,-0.25,0.05, 0,0,0,0.5,0.2, -0.25, 0.12,-0.125,0,0, 0,0,0,0,0)\n",
    "cat('beta_true:', beta_true, '\\n')\n",
    "group_id = rep(1:group_number, each=card_per_group)\n",
    "cat('group_id:', group_id, '\\n')\n",
    "J = group_number\n",
    "group_list = list()\n",
    "for(j in 1:J){\n",
    "  group_list[[j]] = c(1:p)[group_id==j]\n",
    "}\n",
    "group_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fundamental-bronze",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for(j in 1:J){\n",
    "  X_j = X[, group_list[[j]] ]\n",
    "  X_j.svd = svd(X_j)\n",
    "  X[, group_list[[j]] ] = X_j.svd$u\n",
    "}\n",
    "# round(t(X) %*% X, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "planned-baghdad",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "e = rnorm(n, 0, 0.1)\n",
    "y = X %*% beta_true + e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "manual-wesley",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter 1 16.69633 \n",
      "iter 2 6.479991 \n",
      "iter 3 2.08435 \n",
      "iter 4 0.8267275 \n",
      "iter 5 0.2383172 \n",
      "iter 6 0.04213397 \n",
      "iter 7 0.0140898 \n",
      "iter 8 0.003104337 \n",
      "iter 9 0.0006537933 \n",
      "iter 10 0.0001368071 \n",
      "iter 11 2.868036e-05 \n",
      "iter 12 6.026674e-06 \n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>0.145</li><li>-0.293</li><li>0.18</li><li>-0.089</li><li>0.115</li><li>0.073</li><li>0.095</li><li>0.133</li><li>0.356</li><li>0.246</li><li>-0.216</li><li>0.032</li><li>-0.148</li><li>-0.015</li><li>0.068</li><li>0</li><li>0</li><li>0</li><li>0</li><li>0</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 0.145\n",
       "\\item -0.293\n",
       "\\item 0.18\n",
       "\\item -0.089\n",
       "\\item 0.115\n",
       "\\item 0.073\n",
       "\\item 0.095\n",
       "\\item 0.133\n",
       "\\item 0.356\n",
       "\\item 0.246\n",
       "\\item -0.216\n",
       "\\item 0.032\n",
       "\\item -0.148\n",
       "\\item -0.015\n",
       "\\item 0.068\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 0.145\n",
       "2. -0.293\n",
       "3. 0.18\n",
       "4. -0.089\n",
       "5. 0.115\n",
       "6. 0.073\n",
       "7. 0.095\n",
       "8. 0.133\n",
       "9. 0.356\n",
       "10. 0.246\n",
       "11. -0.216\n",
       "12. 0.032\n",
       "13. -0.148\n",
       "14. -0.015\n",
       "15. 0.068\n",
       "16. 0\n",
       "17. 0\n",
       "18. 0\n",
       "19. 0\n",
       "20. 0\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       " [1]  0.145 -0.293  0.180 -0.089  0.115  0.073  0.095  0.133  0.356  0.246\n",
       "[11] -0.216  0.032 -0.148 -0.015  0.068  0.000  0.000  0.000  0.000  0.000"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# hyperparameters\n",
    "lambda = 0.1\n",
    "max_iter = 100\n",
    "threshold = 1e-5\n",
    "\n",
    "beta_hat = rep(10, p)\n",
    "\n",
    "iter = 0\n",
    "\n",
    "while(iter < max_iter){\n",
    "  beta_old = beta_hat\n",
    "  for(j in 1:J){ # for each group\n",
    "    r_j = y - X[, -group_list[[j]] ] %*% beta_hat[ -group_list[[j]] ] # r_j = y - \\Sigma_{k \\neq j} X \\beta_{-j}\n",
    "#     cat(iter, j, r_j, '\\n')\n",
    "#     cat('L2 norm of r_j', sqrt(sum(r_j^2)), '\\n')\n",
    "    X_j = X[, group_list[[j]] ]\n",
    "#     print(round(t(X_j) %*% X,2  ))\n",
    "    \n",
    "    s_j = t(X_j) %*% r_j # s_j = X_j^T ( y - \\Sigma_{k \\neq j} X \\beta_{-j} )\n",
    "#     print(s_j)\n",
    "    p_j = card_per_group\n",
    "    s_2_norm = sqrt( sum(s_j^2) )\n",
    "#     cat('L2 norm of s_j',s_2_norm, '\\n')\n",
    "\n",
    "    if(s_2_norm <= lambda * sqrt(p_j)){\n",
    "      beta_j = 0\n",
    "    }\n",
    "    else{\n",
    "      beta_j = (1 - lambda*sqrt(p_j)/s_2_norm) * s_j\n",
    "    }\n",
    "    beta_hat[ group_list[[j]] ] = beta_j\n",
    "#     print(beta_j)\n",
    "#     cat('L2 norm of beta_j', sqrt(sum(beta_j^2)), '\\n')\n",
    "    \n",
    "  }\n",
    "  iter = iter + 1\n",
    "  change = max( abs( beta_hat - beta_old ) )\n",
    "  cat('iter', iter, change, '\\n')\n",
    "  if(change < threshold){\n",
    "    break\n",
    "  }\n",
    "}\n",
    "round(beta_hat,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fresh-piano",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "theoretical-breakfast",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "adverse-hypothetical",
   "metadata": {},
   "source": [
    "# Sparse group lasso\n",
    "Reference: Friedman, Jerome, Trevor Hastie, and Robert Tibshirani. \"A note on the group lasso and a sparse group lasso.\" arXiv preprint arXiv:1001.0736 (2010).\n",
    "\n",
    "Intuition: sparse group lasso yield sparsity at both the group and individual feature levels."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cognitive-kitty",
   "metadata": {},
   "source": [
    "The objective function is\n",
    "$$min_{\\beta} \\frac12 \\| y-\\Sigma_{l=1}^L X^{(l)}\\beta^{(l)} \\|_2^2 + \\lambda_1 \\Sigma_{l=1}^L \\|\\beta^{(l)}\\|_2 + \\lambda_2 \\|\\beta^{(l)}\\|_1,$$\n",
    "\n",
    "the subgradient equations are\n",
    "$$ -X^{(l)T} (y-\\Sigma_{l=1}^L X^{(l)}\\beta^{(l)}) + \\lambda_1 s_l + \\lambda_2 t_l = 0. $$\n",
    "This is separable for groups so that block coordinate descent can be used, each time update one group only, with other group coefficients fixed.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "improved-making",
   "metadata": {},
   "source": [
    "For each group $l$ whose element $j=1,...,p_j$,\n",
    "$X^{(l)}=(Z_1,...,Z_{p_j}), \\beta^{(l)} = \\theta = (\\theta_1,...,\\theta_{p_j})^T$, let $r = y - \\Sigma_{i\\neq l} X^{(i)}\\beta^{(i)}$,\n",
    "\n",
    "the subgradient equations are\n",
    "$$ -Z_j^T(r - \\Sigma_k Z_k\\theta_k) + \\lambda_1 s_j + \\lambda_2 t_j = 0 $$\n",
    "\n",
    "Let $a=X^{(l)T} r$,\n",
    "\n",
    "$$\\theta=0 \\iff a_j=\\lambda_1 s_j + \\lambda_2 t_j \\text{ with } \\|s\\|\\leq1, t_j\\in[-1,1]$$\n",
    "\n",
    "$$\\iff min_t J(t) = \\Sigma_{j=1}^{p_j}s_j^2 = \\frac{1}{\\lambda_1^2} \\Sigma_{j=1}^{p_j} (a_j-\\lambda_2 t_j)^2, \\text{ s.t. }t_j\\in[-1,1],$$\n",
    "whose solution is $\\hat t_j=\\frac{a_j}{\\lambda_2}$ if $|\\frac{a_j}{\\lambda_2}|\\leq1$, otherwise $\\hat t_j=sign(\\frac{a_j}{\\lambda_2})$.\n",
    "\n",
    "If $\\hat t_j$ is a feasible solution with $J(\\hat t_j)\\leq1$, then $\\theta=0$.\n",
    "\n",
    "If $J(\\hat t_j)>1$, optimize each nonzero entries in $\\theta$ though\n",
    "$$ min_{\\theta_j} \\frac12 \\| r - \\Sigma_{j=1}^{p_j} Z_j \\theta_j \\|_2^2 + \\lambda_1 \\| \\theta \\|_2 + \\lambda_2 \\Sigma_{j=1}^{p_j} |\\theta_j|$$\n",
    "by one-dimensional search with **optimize** function in the R package"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "steady-accident",
   "metadata": {},
   "source": [
    "Algorithm: \n",
    "1. outer loop: block coordinate decent over group $l=1,...,L$,\n",
    "2. inner loop: coordinate decent over (nonzero) entries $\\theta_j$ for $j=1,...,p_j$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "effective-boxing",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "internal-exhibition",
   "metadata": {},
   "source": [
    "Simulation setting:\n",
    "$n=100, p=20$, $X$ is centered and scaled to variance $\\frac1p$ in each column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "diverse-shannon",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100\n",
    "p = 20\n",
    "set.seed(0)\n",
    "\n",
    "X = matrix( rnorm(n*p, mean = 0, sd = 1), n, p )\n",
    "X = X - colMeans(X) # centering\n",
    "Xsd <- sqrt(colMeans(X^2))\n",
    "X = t( t(X)/Xsd ) /sqrt(p)\n",
    "\n",
    "e = rnorm(n, 0, 0.1)\n",
    "y = X %*% beta_true + e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "scientific-blood",
   "metadata": {},
   "outputs": [],
   "source": [
    "group_number = 4\n",
    "card_per_group = p / group_number\n",
    "beta_true = c(0.15,-0.33,0.25,-0.25,0.05, 0,0,0,0.5,0.2, -0.25, 0.12,-0.125,0,0, 0,0,0,0,0)\n",
    "# cat('beta:', beta, '\\n')\n",
    "group_id = rep(1:group_number, each=card_per_group)\n",
    "# cat('group_id:', group_id, '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "earlier-compact",
   "metadata": {},
   "outputs": [],
   "source": [
    "L = group_number\n",
    "group_list = list()\n",
    "for(l in 1:L){\n",
    "  group_list[[l]] = c(1:p)[group_id==l]\n",
    "}\n",
    "# group_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "tough-check",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1]  2.719209e-02 -3.893512e-01  3.046806e-01 -1.735826e-01  1.346569e-06\n",
      "[1] -0.01562881 -0.12591371  0.00000000  0.45096282  0.08911321\n",
      "[1] -0.18672487  0.06802142 -0.08516696  0.04588338  0.00000000\n",
      "[1]  0.00000e+00  0.00000e+00  0.00000e+00 -8.22812e-06  0.00000e+00\n",
      "iter 1 0.4509628 \n",
      "[1]  0.03061995 -0.35211492  0.25070759 -0.18414058  0.00000000\n",
      "[1] -5.228342e-05 -7.127099e-02  0.000000e+00  4.413935e-01  6.923042e-02\n",
      "[1] -0.19879879  0.07278474 -0.09648731  0.05922917  0.00000000\n",
      "[1]  0.000000e+00  0.000000e+00  0.000000e+00 -3.698008e-05 -7.812010e-06\n",
      "iter 2 0.05464272 \n",
      "[1]  0.02954488 -0.34953557  0.24797525 -0.18503279  0.00000000\n",
      "[1] -6.886815e-06 -6.584137e-02  0.000000e+00  4.393237e-01  6.845442e-02\n",
      "[1] -0.19987972  0.07360516 -0.09757684  0.06039204  0.00000000\n",
      "[1]  0.000000e+00  0.000000e+00  0.000000e+00  2.452883e-07 -3.605729e-05\n",
      "iter 3 0.005429615 \n",
      "[1]  0.02874492 -0.34914443  0.24763086 -0.18495681  0.00000000\n",
      "[1] -1.001736e-05 -6.522575e-02  0.000000e+00  4.390715e-01  6.860152e-02\n",
      "[1] -0.20009598  0.07378574 -0.09766248  0.06058278  0.00000000\n",
      "[1]  0.000000e+00  0.000000e+00  0.000000e+00  4.362841e-08 -1.289102e-05\n",
      "iter 4 0.0007999633 \n",
      "[1]  0.02856973 -0.34910394  0.24755630 -0.18492607  0.00000000\n",
      "[1] -1.350522e-05 -6.512050e-02  0.000000e+00  4.390371e-01  6.863401e-02\n",
      "[1] -0.20014428  0.07381797 -0.09767460  0.06061230  0.00000000\n",
      "[1]  0.000000e+00  0.000000e+00  0.000000e+00 -8.090219e-07  3.439071e-06\n",
      "iter 5 0.0001751909 \n",
      "[1]  0.02853965 -0.34909695  0.24754246 -0.18492012  0.00000000\n",
      "[1] -1.406752e-05 -6.510099e-02  0.000000e+00  4.390319e-01  6.863832e-02\n",
      "[1] -0.20015480  0.07382259 -0.09767829  0.06061572  0.00000000\n",
      "[1]  0.000000e+00  0.000000e+00  0.000000e+00 -1.021963e-06 -3.889179e-07\n",
      "iter 6 3.007139e-05 \n",
      "[1]  0.02853288 -0.34909524  0.24753940 -0.18491877  0.00000000\n",
      "[1] -1.404733e-05 -6.509779e-02  0.000000e+00  4.390303e-01  6.863836e-02\n",
      "[1] -0.20015627  0.07382388 -0.09767799  0.06061758  0.00000000\n",
      "[1]  0.000000e+00  0.000000e+00  0.000000e+00 -1.074452e-06 -1.613236e-06\n",
      "iter 7 6.779386e-06 \n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>0.029</li><li>-0.349</li><li>0.248</li><li>-0.185</li><li>0</li><li>0</li><li>-0.065</li><li>0</li><li>0.439</li><li>0.069</li><li>-0.2</li><li>0.074</li><li>-0.098</li><li>0.061</li><li>0</li><li>0</li><li>0</li><li>0</li><li>0</li><li>0</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 0.029\n",
       "\\item -0.349\n",
       "\\item 0.248\n",
       "\\item -0.185\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item -0.065\n",
       "\\item 0\n",
       "\\item 0.439\n",
       "\\item 0.069\n",
       "\\item -0.2\n",
       "\\item 0.074\n",
       "\\item -0.098\n",
       "\\item 0.061\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 0.029\n",
       "2. -0.349\n",
       "3. 0.248\n",
       "4. -0.185\n",
       "5. 0\n",
       "6. 0\n",
       "7. -0.065\n",
       "8. 0\n",
       "9. 0.439\n",
       "10. 0.069\n",
       "11. -0.2\n",
       "12. 0.074\n",
       "13. -0.098\n",
       "14. 0.061\n",
       "15. 0\n",
       "16. 0\n",
       "17. 0\n",
       "18. 0\n",
       "19. 0\n",
       "20. 0\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       " [1]  0.029 -0.349  0.248 -0.185  0.000  0.000 -0.065  0.000  0.439  0.069\n",
       "[11] -0.200  0.074 -0.098  0.061  0.000  0.000  0.000  0.000  0.000  0.000"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# hyperparameters\n",
    "lambda.1 = 0.1\n",
    "lambda.2 = 0.2\n",
    "max_iter = 100\n",
    "threshold = 1e-5\n",
    "\n",
    "beta_hat = rep(0, p)\n",
    "\n",
    "iter = 0\n",
    "\n",
    "while(iter < max_iter){\n",
    "  beta_old = beta_hat\n",
    "  for(l in 1:L){ # for each group\n",
    "    r_l = y - X[, -group_list[[l]] ] %*% beta_hat[ -group_list[[l]] ] # r_l = y - \\Sigma_{k \\neq j} X \\beta_{-j}\n",
    "#     cat(iter, l, r_l, '\\n')\n",
    "    X_l = X[, group_list[[l]] ]\n",
    "    a = t(X_l) %*% r_l\n",
    "\n",
    "#     a_l2 = a / lambda.2\n",
    "# #     print(a_l2)\n",
    "#     t = a_l2\n",
    "#     t[t>1] = 1\n",
    "#     t[t<(-1)] = -1\n",
    "# #     print(t)\n",
    "#     a_l2t_2 = (a - lambda.2*t)^2 # (a_j - \\lambda_2 t_j)^2\n",
    "# #     print(a_l2t_2)\n",
    "#     J = sum(a_l2t_2) / (lambda.1**2)\n",
    "#     print(J)\n",
    "#     if(J <= 1){\n",
    "#       theta = 0\n",
    "#     }\n",
    "\n",
    "    p_l = card_per_group\n",
    "    theta = rep(0, p_l)\n",
    "    for(j in 1:p_l){\n",
    "      if(abs(a[j]) <= lambda.2){\n",
    "        theta[j] = 0\n",
    "      }\n",
    "      else{\n",
    "        Z_j = X_l[, -j] # Z_{-j}\n",
    "        theta_j = theta[-j] # \\theta_{-j}\n",
    "        f <- function(x) { 1/2* sum(( r_l - Z_j%*%theta_j - X_l[,j]*x )^2) +\n",
    "                          lambda.1 * sqrt(sum(theta_j**2) + x**2) +\n",
    "                          lambda.2 * (sum(abs(theta_j)) + abs(x))  }\n",
    "        theta[j] = optimize(f, c(-10,10))$minimum\n",
    "      }\n",
    "    }\n",
    "    print(theta)\n",
    "    beta_hat[ group_list[[l]] ] = theta\n",
    "  \n",
    "  }\n",
    "\n",
    "\n",
    "\n",
    "  iter = iter + 1\n",
    "  change = max( abs( beta_hat - beta_old ) )\n",
    "  cat('iter', iter, change, '\\n')\n",
    "  if(change < threshold){\n",
    "    break\n",
    "  }\n",
    "}\n",
    "round(beta_hat,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "executive-milton",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "scenic-output",
   "metadata": {},
   "outputs": [],
   "source": [
    "# library(SGL)\n",
    "# SGL(data = list(x=X, y=y), index = group_id)\n",
    "# SGL(data = list(x=X, y=y), index = group_id, lambdas = (lambda.1+lambda.2), alpha=lambda.2/(lambda.1+lambda.2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aerial-least",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "stable-insulation",
   "metadata": {},
   "source": [
    "Conclusion:\n",
    "1. Group lasso put assumption on groupwise sparsity and successfully reveal all zeroes in the last group, but for nonzero groups like group 2 and group 3 in the simulation setting, elements are shrinkaged to be closer, are all nonzero in these groups.\n",
    "2. Sparse group lasso handles both groupwise sparsity and the sparsity within each group, which performs better especially for the case when both zero and nonzero elements exist within group, just like group 2 and group 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "vital-mountain",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}