Comparison to regular sparse group LASSO

Compare directly to SGL.ipynb

Compare to sglasso.ipynb

Comparison to regular sparse group LASSO.ipynb

Cox stacked block - Ruilin penalty.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "import numpy as np, time\n",
    "from regreg.smooth.cox import cox_loglike\n",
    "import regreg.api as rr\n",
    "import regreg.affine as ra\n",
    "%load_ext rpy2.ipython\n",
    "toc = time.time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: Loading required package: rms\n",
      "\n",
      "R[write to console]: Loading required package: Hmisc\n",
      "\n",
      "R[write to console]: Loading required package: lattice\n",
      "\n",
      "R[write to console]: Loading required package: survival\n",
      "\n",
      "R[write to console]: Loading required package: Formula\n",
      "\n",
      "R[write to console]: Loading required package: ggplot2\n",
      "\n",
      "R[write to console]: \n",
      "Attaching package: ‘Hmisc’\n",
      "\n",
      "\n",
      "R[write to console]: The following objects are masked from ‘package:base’:\n",
      "\n",
      "    format.pval, units\n",
      "\n",
      "\n",
      "R[write to console]: Loading required package: SparseM\n",
      "\n",
      "R[write to console]: \n",
      "Attaching package: ‘SparseM’\n",
      "\n",
      "\n",
      "R[write to console]: The following object is masked from ‘package:base’:\n",
      "\n",
      "    backsolve\n",
      "\n",
      "\n",
      "R[write to console]: Loading required package: mgcv\n",
      "\n",
      "R[write to console]: Loading required package: nlme\n",
      "\n",
      "R[write to console]: This is mgcv 1.8-28. For overview type 'help(\"mgcv-package\")'.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%R library(coxed)#install.packages('coxed', repos='http://cloud.r-project.org')\n",
    "toc = time.time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "K = 20\n",
    "p = 5000\n",
    "X_list = list()\n",
    "censor_list = list()\n",
    "y_list = list()\n",
    "s = 100\n",
    "true_beta = matrix(rnorm(K*p),K,p)\n",
    "\n",
    "for (i in 1:K){\n",
    "  N = (1000+10*i)\n",
    "  X = matrix(as.numeric(rbinom(N*p, 1, 0.5)), N, p)\n",
    "  simdata = sim.survdata(T=120, num.data.frames=1, X = X, beta=true_beta[i,])\n",
    "  df = simdata$data\n",
    "  y =  df$y\n",
    "  censor_list[[i]] = as.numeric(df$failed)\n",
    "  X_list[[i]] = X\n",
    "  y_list[[i]] = df$y\n",
    "}\n",
    "save(X_list, y_list, censor_list, file='instance.RData')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data generation time: 48.116044998168945\n"
     ]
    }
   ],
   "source": [
    "tic = time.time()\n",
    "print('Data generation time:', tic-toc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading time: 4.379420042037964\n"
     ]
    }
   ],
   "source": [
    "toc = time.time()\n",
    "%R load('instance.RData')\n",
    "def load_data(idx):\n",
    "    %R -i idx -o X X = X_list[[idx]]\n",
    "    %R -o Y Y = y_list[[idx]]\n",
    "    %R -o C C = censor_list[[idx]]\n",
    "    return X, Y, C\n",
    "datasets = [load_data(idx) for idx in range(1, 21)]\n",
    "tic = time.time()\n",
    "print('Loading time:', tic-toc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "losses = [rr.cox_loglike(Y.shape[0], Y.reshape(-1), C.reshape(-1), coef=1./Y.shape[0]) for _, Y, C in datasets]\n",
    "Xblock = rr.block_columns([X for X, _, _ in datasets])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((5000, 20), (22100,))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xblock.input_shape, Xblock.output_shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(22100,)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "class cox_stacked(rr.smooth_atom):\n",
    "\n",
    "    def __init__(self,\n",
    "                 losses,\n",
    "                 X,\n",
    "                 quadratic=None, \n",
    "                 initial=None,\n",
    "                 offset=None):\n",
    "        \n",
    "        self.losses = losses\n",
    "        self.ndisease = len(losses)\n",
    "        self.nfeature = X.shape[0]\n",
    "\n",
    "        self.X, self.X_T = X, X.T\n",
    "        \n",
    "        rr.smooth_atom.__init__(self,\n",
    "                                self.X.output_shape,\n",
    "                                offset=offset,\n",
    "                                quadratic=quadratic,\n",
    "                                initial=initial)\n",
    "        self._gradient = np.zeros(X.output_shape)\n",
    "\n",
    "    def smooth_objective(self, arg, mode='both', check_feasibility=False):\n",
    "\n",
    "        arg = self.apply_offset(arg) # (nfeature, ndisease)\n",
    "        linpred = self.X.dot(arg) # (ndisease, ncase)\n",
    "        if mode == 'grad':\n",
    "            for d, slice in enumerate(self.X._slices):\n",
    "                self._gradient[slice] = self.losses[d].smooth_objective(linpred[slice], 'grad')\n",
    "            return self.scale(self.X_T.dot(self._gradient))\n",
    "        elif mode == 'func':\n",
    "            value = 0\n",
    "            for d, slice in enumerate(self.X._slices):\n",
    "                value += self.losses[d].smooth_objective(linpred[slice], 'func')\n",
    "            return self.scale(value)\n",
    "        elif mode == 'both':\n",
    "            value = 0\n",
    "            for d, slice in enumerate(self.X._slices):\n",
    "                f, g = self.losses[d].smooth_objective(linpred[slice], 'both')\n",
    "                self._gradient[slice] = g\n",
    "                value += f\n",
    "            return self.scale(value), self.scale(self.X_T.dot(self._gradient))\n",
    "        else:\n",
    "            raise ValueError(\"mode incorrectly specified\")\n",
    "\n",
    "loss = cox_stacked(losses, Xblock)\n",
    "loss.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check the loss can be computed\n",
    "\n",
    "- We'll use `G` to compute $\\lambda_{\\max}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5000, 20)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F, G = loss.smooth_objective(np.zeros(Xblock.input_shape), 'both')\n",
    "G.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.020851269364356995"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nfeature = Xblock.input_shape[0]\n",
    "penalty = rr.sparse_group_block(Xblock.input_shape, l1_weight=1, \n",
    "                                l2_weight=1, lagrange=1.)\n",
    "dual = penalty.conjugate\n",
    "lambda_max = dual.seminorm(G, lagrange=1)\n",
    "penalty.lagrange = lambda_max\n",
    "lambda_max"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.99999"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "problem = rr.simple_problem(loss, penalty)\n",
    "soln = problem.solve(tol=1.e-9)\n",
    "np.mean(soln == 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First 10 values on logscale of length 100 down to 0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.02085127, 0.01990355, 0.0189989 , 0.01813537, 0.01731109,\n",
       "       0.01652427, 0.01577322, 0.0150563 , 0.01437197, 0.01371874])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lagrange_vals = np.exp(np.linspace(0, np.log(0.01), 100))[:10] * lambda_max\n",
    "lagrange_vals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Timing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 1 0.020851269364356995 0.020851269364356995\n",
      "2 4 0.019901469349861145 0.01990354699118715\n",
      "4 6 0.018999144434928894 0.018998900062533984\n",
      "6 10 0.01813565194606781 0.018135370732964236\n",
      "10 17 0.017311230301856995 0.017311090144141207\n",
      "18 38 0.016524866223335266 0.016524274380223875\n",
      "41 93 0.015773668885231018 0.015773220607099382\n",
      "72 190 0.015056952834129333 0.015056303387093354\n",
      "109 309 0.01437290757894516 0.014371971161182306\n",
      "194 588 0.013719134032726288 0.013718742891094956\n",
      "time: 41.756864\n"
     ]
    }
   ],
   "source": [
    "from time import time\n",
    "toc = time()\n",
    "solns = []\n",
    "problem.coefs[:] = 0\n",
    "for lagrange in lagrange_vals:\n",
    "    penalty.lagrange = lagrange\n",
    "    soln = problem.solve(tol=1.e-12)\n",
    "    solns.append(soln.copy())\n",
    "    print(np.sum(np.sum(soln**2, 1) > 0), np.sum(soln != 0), dual.seminorm(loss.smooth_objective(soln, 'grad'), lagrange=1.), lagrange)\n",
    "tic = time()\n",
    "print('time: %f' % (tic-toc))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "all,-slideshow",
   "formats": "ipynb,Rmd"
  },
  "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.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

Cox stacked block.ipynb

Cox stacked.ipynb

Generating an instance.ipynb