<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/classification-from-scratch-penalized-lasso-logistic-5-8/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://freakonometrics.hypotheses.org/52894"> R-english – Freakonometrics</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Fifth post of our series on <a rel="nofollow" href="https://freakonometrics.hypotheses.org/52731" target="_blank">classification from scratch</a>, following the previous post on penalization using the <img alt="\ell_2" title="\ell_2" src="https://s0.wp.com/latex.php?latex=%5Cell_2&bg=ffffff&fg=000&s=0" class="latex"> norm (so-called <a rel="nofollow" href="https://freakonometrics.hypotheses.org/52773" target="_blank">Ridge regression</a>), this time, we will discuss penalization based on the <img alt="\ell_1" title="\ell_1" src="https://s0.wp.com/latex.php?latex=%5Cell_1&bg=ffffff&fg=000&s=0" class="latex"> norm (the so-called Lasso regression).</p>
<p>First of all, one should admit that if the name stands for <a rel="nofollow" href="https://en.wikipedia.org/wiki/Lasso_(statistics)" target="_blank">least absolute shrinkage and selection operator</a>, that’s actually a very cool name… Funny story, a few years before, Leo Breiman introduce a concept of <a rel="nofollow" href="https://t.co/THHHYIwgCL" target="_blank">garrote technique</a>… “The garrote eliminates some variables, shrinks others, and is relatively stable”.</p>
<p><img alt="" src="https://i2.wp.com/f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2018/06/DejmDd2XUAAR5_A.jpeg?w=456&ssl=1" data-recalc-dims="1"></p>
<p>I guess that somehow, the lasso is the extension of the garotte technique</p>
<h2>Normalization of the covariates</h2>
<p>As <a rel="nofollow" href="https://freakonometrics.hypotheses.org/52773" target="_blank">previously</a>, the first step will be to consider linear transformations of all covariates <img alt="x_j" title="x_j" src="https://s0.wp.com/latex.php?latex=x_j&bg=ffffff&fg=000&s=0" class="latex"> to get centered and scaled variables (with unit variance)</p>
<div><table><tr><td><pre>1
2
3
4
</pre></td><td><pre>y <span>=</span> myocarde$PRONO
X <span>=</span> myocarde<span>[</span>,<span>1</span><span>:</span><span>7</span><span>]</span>
<span>for</span><span>(</span>j <span>in</span> <span>1</span><span>:</span><span>7</span><span>)</span> X<span>[</span>,j<span>]</span> <span>=</span> <span>(</span>X<span>[</span>,j<span>]</span><span>-</span><span>mean</span><span>(</span>X<span>[</span>,j<span>]</span><span>)</span><span>)</span><span>/</span><span>sd</span><span>(</span>X<span>[</span>,j<span>]</span><span>)</span>
X <span>=</span> <span>as.<span>matrix</span></span><span>(</span>X<span>)</span></pre></td></tr></table></div>
<h2>Ridge Regression (from scratch)</h2>
<p>The heuristics about Lasso regression is the following graph. In the background, we can visualize the (two-dimensional) log-likelihood of the logistic regression, and the blue square is the constraint we have, if we rewite the optimization problem as a contrained optimization problem,</p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
8
9
10
</pre></td><td><pre>LogLik <span>=</span> <span>function</span><span>(</span>bbeta<span>)</span><span>{</span>
  b0<span>=</span>bbeta<span>[</span><span>1</span><span>]</span>
  <span>beta</span><span>=</span>bbeta<span>[</span><span>-</span><span>1</span><span>]</span>
  <span>sum</span><span>(</span><span>-</span>y<span>*</span><span>log</span><span>(</span><span>1</span> <span>+</span> <span>exp</span><span>(</span><span>-</span><span>(</span>b0<span>+</span>X<span>%*%</span><span>beta</span><span>)</span><span>)</span><span>)</span> <span>-</span> 
  <span>(</span><span>1</span><span>-</span>y<span>)</span><span>*</span><span>log</span><span>(</span><span>1</span> <span>+</span> <span>exp</span><span>(</span>b0<span>+</span>X<span>%*%</span><span>beta</span><span>)</span><span>)</span><span>)</span><span>}</span>
u <span>=</span> <span>seq</span><span>(</span><span>-</span><span>4</span>,<span>4</span>,<span>length</span><span>=</span><span>251</span><span>)</span>
v <span>=</span> <span>outer</span><span>(</span>u,u,<span>function</span><span>(</span>x,y<span>)</span> LogLik<span>(</span><span>c</span><span>(</span><span>1</span>,x,y<span>)</span><span>)</span><span>)</span>
<span>image</span><span>(</span>u,u,v,<span>col</span><span>=</span><span>rev</span><span>(</span><span>heat.<span>colors</span></span><span>(</span><span>25</span><span>)</span><span>)</span><span>)</span>
<span>contour</span><span>(</span>u,u,v,add<span>=</span>TRUE<span>)</span>
<span>polygon</span><span>(</span><span>c</span><span>(</span><span>-</span><span>1</span>,<span>0</span>,<span>1</span>,<span>0</span><span>)</span>,<span>c</span><span>(</span><span>0</span>,<span>1</span>,<span>0</span>,<span>-</span><span>1</span><span>)</span>,border<span>=</span><span>&quot;blue&quot;</span><span>)</span></pre></td></tr></table></div>
<p><img alt="" src="https://i1.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-01-a%CC%80-15.30.15.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>The nice thing here is that is works as a variable selection tool, since some components can be null here. That’s the idea behind the following (popular) graph</p>
<p><img alt="" src="https://i2.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-04-a%CC%80-14.53.39.png?w=456&ssl=1" data-recalc-dims="1"><br>
(with lasso on the left, and ridge on the right).</p>
<p>Heuristically, the maths explanation is the following. Consider a simple regression <img alt="y_i=x_i\beta+\varepsilon" title="y_i=x_i\beta+\varepsilon" src="https://s0.wp.com/latex.php?latex=y_i%3Dx_i%5Cbeta%2B%5Cvarepsilon&bg=ffffff&fg=000&s=0" class="latex">, with <img alt="\ell_1" title="\ell_1" src="https://s0.wp.com/latex.php?latex=%5Cell_1&bg=ffffff&fg=000&s=0" class="latex">-penality and a <img alt="\ell_2" title="\ell_2" src="https://s0.wp.com/latex.php?latex=%5Cell_2&bg=ffffff&fg=000&s=0" class="latex">-loss fuction. The optimization problem becomes<img alt="\min\big\{\mathbf{y}^T\mathbf{y}-2\mathbf{y}^T\mathbf{x}\beta+\beta\mathbf{x}^T\mathbf{x}\beta+2\lambda{\color{red}{|}}\beta{\color{red}{|}}\big\}" title="\min\big\{\mathbf{y}^T\mathbf{y}-2\mathbf{y}^T\mathbf{x}\beta+\beta\mathbf{x}^T\mathbf{x}\beta+2\lambda{\color{red}{|}}\beta{\color{red}{|}}\big\}" src="https://s0.wp.com/latex.php?latex=%5Cmin%5Cbig%5C%7B%5Cmathbf%7By%7D%5ET%5Cmathbf%7By%7D-2%5Cmathbf%7By%7D%5ET%5Cmathbf%7Bx%7D%5Cbeta%2B%5Cbeta%5Cmathbf%7Bx%7D%5ET%5Cmathbf%7Bx%7D%5Cbeta%2B2%5Clambda%7B%5Ccolor%7Bred%7D%7B%7C%7D%7D%5Cbeta%7B%5Ccolor%7Bred%7D%7B%7C%7D%7D%5Cbig%5C%7D&bg=ffffff&fg=000&s=0" class="latex">The first order condition can be written<img alt="-2\mathbf{y}^T\mathbf{x}+2\mathbf{x}^T\mathbf{x}\widehat{\beta}{\color{red}{\pm} }2\lambda=0" title="-2\mathbf{y}^T\mathbf{x}+2\mathbf{x}^T\mathbf{x}\widehat{\beta}{\color{red}{\pm} }2\lambda=0" src="https://s0.wp.com/latex.php?latex=-2%5Cmathbf%7By%7D%5ET%5Cmathbf%7Bx%7D%2B2%5Cmathbf%7Bx%7D%5ET%5Cmathbf%7Bx%7D%5Cwidehat%7B%5Cbeta%7D%7B%5Ccolor%7Bred%7D%7B%5Cpm%7D+%7D2%5Clambda%3D0&bg=ffffff&fg=000&s=0" class="latex">(the sign in <img alt="{\color{red}{\pm}}" title="{\color{red}{\pm}}" src="https://s0.wp.com/latex.php?latex=%7B%5Ccolor%7Bred%7D%7B%5Cpm%7D%7D&bg=ffffff&fg=000&s=0" class="latex"> being the sign of <img alt="\widehat{\beta}" title="\widehat{\beta}" src="https://s0.wp.com/latex.php?latex=%5Cwidehat%7B%5Cbeta%7D&bg=ffffff&fg=000&s=0" class="latex">).<br>
Assume that <img alt="\mathbf{y}^T\mathbf{x}&gt;0" title="\mathbf{y}^T\mathbf{x}&gt;0" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7By%7D%5ET%5Cmathbf%7Bx%7D%3E0&bg=ffffff&fg=000&s=0" class="latex">, then solution is<br>
<img alt="\widehat{\beta}_{\lambda}^{lasso}=\max\left\lbrace\frac{\mathbf{y}^T\mathbf{x}-\lambda}{\mathbf{x}^T\mathbf{x}},0\right\rbrace" title="\widehat{\beta}_{\lambda}^{lasso}=\max\left\lbrace\frac{\mathbf{y}^T\mathbf{x}-\lambda}{\mathbf{x}^T\mathbf{x}},0\right\rbrace" src="https://s0.wp.com/latex.php?latex=%5Cwidehat%7B%5Cbeta%7D_%7B%5Clambda%7D%5E%7Blasso%7D%3D%5Cmax%5Cleft%5Clbrace%5Cfrac%7B%5Cmathbf%7By%7D%5ET%5Cmathbf%7Bx%7D-%5Clambda%7D%7B%5Cmathbf%7Bx%7D%5ET%5Cmathbf%7Bx%7D%7D%2C0%5Cright%5Crbrace&bg=ffffff&fg=000&s=0" class="latex">(we get a corner solution when <img alt="\lambda" title="\lambda" src="https://s0.wp.com/latex.php?latex=%5Clambda&bg=ffffff&fg=000&s=0" class="latex"> is large).</p>
<h2>Optimization routine</h2>
<p>As in our previous post, let us start with standard (R) optimization routines, such as <a rel="nofollow" href="https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm" target="_blank">BFGS</a></p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
</pre></td><td><pre>PennegLogLik <span>=</span> <span>function</span><span>(</span>bbeta,lambda<span>=</span><span>0</span><span>)</span><span>{</span>
  b0<span>=</span>bbeta<span>[</span><span>1</span><span>]</span>
  <span>beta</span><span>=</span>bbeta<span>[</span><span>-</span><span>1</span><span>]</span>
 <span>-</span><span>sum</span><span>(</span><span>-</span>y<span>*</span><span>log</span><span>(</span><span>1</span> <span>+</span> <span>exp</span><span>(</span><span>-</span><span>(</span>b0<span>+</span>X<span>%*%</span><span>beta</span><span>)</span><span>)</span><span>)</span> <span>-</span> 
<span>(</span><span>1</span><span>-</span>y<span>)</span><span>*</span><span>log</span><span>(</span><span>1</span> <span>+</span> <span>exp</span><span>(</span>b0<span>+</span>X<span>%*%</span><span>beta</span><span>)</span><span>)</span><span>)</span><span>+</span>lambda<span>*</span><span>sum</span><span>(</span><span>abs</span><span>(</span><span>beta</span><span>)</span><span>)</span>
<span>}</span>
opt_lasso <span>=</span> <span>function</span><span>(</span>lambda<span>)</span><span>{</span>
beta_init <span>=</span> <span>lm</span><span>(</span>PRONO~.,<span>data</span><span>=</span>myocarde<span>)</span>$coefficients
logistic_opt <span>=</span> <span>optim</span><span>(</span><span>par</span> <span>=</span> beta_init<span>*</span><span>0</span>, <span>function</span><span>(</span>x<span>)</span> PennegLogLik<span>(</span>x,lambda<span>)</span>, 
hessian<span>=</span>TRUE, method <span>=</span> <span>&quot;BFGS&quot;</span>, control<span>=</span><span>list</span><span>(</span>abstol<span>=</span>1e<span>-</span>9<span>)</span><span>)</span>
logistic_opt$par<span>[</span><span>-</span><span>1</span><span>]</span>
<span>}</span>
v_lambda<span>=</span><span>c</span><span>(</span><span>exp</span><span>(</span><span>seq</span><span>(</span><span>-</span><span>4</span>,<span>2</span>,<span>length</span><span>=</span><span>61</span><span>)</span><span>)</span><span>)</span>
est_lasso<span>=</span><span>Vectorize</span><span>(</span>opt_lasso<span>)</span><span>(</span>v_lambda<span>)</span>
<span>library</span><span>(</span><span>&quot;RColorBrewer&quot;</span><span>)</span>
colrs<span>=</span>brewer.<span>pal</span><span>(</span><span>7</span>,<span>&quot;Set1&quot;</span><span>)</span>
<span>plot</span><span>(</span>v_lambda,est_lasso<span>[</span><span>1</span>,<span>]</span>,<span>col</span><span>=</span>colrs<span>[</span><span>1</span><span>]</span>,type<span>=</span><span>&quot;l&quot;</span><span>)</span>
<span>for</span><span>(</span>i <span>in</span> <span>2</span><span>:</span><span>7</span><span>)</span> <span>lines</span><span>(</span>v_lambda,est_lasso<span>[</span>i,<span>]</span>,<span>col</span><span>=</span>colrs<span>[</span>i<span>]</span>,lwd<span>=</span><span>2</span><span>)</span></pre></td></tr></table></div>
<p><img alt="" src="https://i2.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-02-a%CC%80-22.38.08.png?w=456&ssl=1" data-recalc-dims="1"><br>
But it is very heratic… or non stable.</p>
<h2>Using glmnet</h2>
<p>Just to compare, with R routines dedicated to lasso, we get the following</p>
<div><table><tr><td><pre>1
2
3
</pre></td><td><pre><span>library</span><span>(</span>glmnet<span>)</span>
glm_lasso <span>=</span> glmnet<span>(</span>X, y, alpha<span>=</span><span>1</span><span>)</span>
<span>plot</span><span>(</span>glm_lasso,xvar<span>=</span><span>&quot;lambda&quot;</span>,<span>col</span><span>=</span>colrs,lwd<span>=</span><span>2</span><span>)</span></pre></td></tr></table></div>
<p><img alt="" src="https://i0.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-02-a%CC%80-22.37.35.png?w=456&ssl=1" data-recalc-dims="1"></p>
<div><table><tr><td><pre>1
</pre></td><td><pre><span>plot</span><span>(</span>glm_lasso,<span>col</span><span>=</span>colrs,lwd<span>=</span><span>2</span><span>)</span></pre></td></tr></table></div>
<p><img alt="" src="https://i2.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-02-a%CC%80-22.37.49.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>If we look carefully what’s in the ouput, we can see that there is variable selection, in the sense that some <img alt="\widehat{\beta}_{j,\lambda}=0" title="\widehat{\beta}_{j,\lambda}=0" src="https://s0.wp.com/latex.php?latex=%5Cwidehat%7B%5Cbeta%7D_%7Bj%2C%5Clambda%7D%3D0&bg=ffffff&fg=000&s=0" class="latex">, in the sense “really null”</p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
8
9
10
</pre></td><td><pre>glmnet<span>(</span>X, y, alpha<span>=</span><span>1</span>,lambda<span>=</span><span>exp</span><span>(</span><span>-</span><span>4</span><span>)</span><span>)</span>$beta
7x1 sparse Matrix of <span>class</span> <span>&quot;dgCMatrix&quot;</span>
               s0
FRCAR  .         
<span>INCAR</span>  <span>0.11005070</span>
INSYS  <span>0.03231929</span>
PRDIA  .         
<span>PAPUL</span>  .         
<span>PVENT</span> <span>-</span><span>0.03138089</span>
REPUL <span>-</span><span>0.20962611</span></pre></td></tr></table></div>
<p>Of course, with out optimization routine, we cannot expect to have null values</p>
<div><table><tr><td><pre>1
2
3
4
5
</pre></td><td><pre>opt_lasso<span>(</span>.2<span>)</span>
         FRCAR         INCAR         INSYS         PRDIA
  <span>0.4810999782</span>  <span>0.0002813658</span>  <span>1.9117847987</span> <span>-</span><span>0.3873926427</span>
          PAPUL         PVENT        REPUL 
 <span>-</span><span>0.0863050787</span> <span>-</span><span>0.4144139379</span> <span>-</span><span>1.3849264055</span></pre></td></tr></table></div>
<p>So clearly, it will be necessary to spend more time today, to understand how it works…</p>
<h2>Orthogonal covariates</h2>
<p>Before getting into the maths, observe that when covariates are orthogonal, there is some very clear “variable” selection process,</p>
<div><table><tr><td><pre>1
2
3
4
5
6
</pre></td><td><pre><span>library</span><span>(</span>factoextra<span>)</span>
pca <span>=</span> <span>princomp</span><span>(</span>X<span>)</span>
pca_X <span>=</span> get_pca_ind<span>(</span>pca<span>)</span>$coord
glm_lasso <span>=</span> glmnet<span>(</span>pca_X, y, alpha<span>=</span><span>1</span><span>)</span>
<span>plot</span><span>(</span>glm_lasso,xvar<span>=</span><span>&quot;lambda&quot;</span>,<span>col</span><span>=</span>colrs<span>)</span>
<span>plot</span><span>(</span>glm_lasso,<span>col</span><span>=</span>colrs<span>)</span></pre></td></tr></table></div>
<p><img alt="" src="https://i1.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-02-a%CC%80-23.25.08.png?w=456&ssl=1" data-recalc-dims="1"></p>
<h2>Interior Point approach</h2>
<p>The penalty is now expressed using the <img alt="\ell_1" title="\ell_1" src="https://s0.wp.com/latex.php?latex=%5Cell_1&bg=ffffff&fg=000&s=0" class="latex"> so intuitively, it should be possible to consider algorithms related to linear programming. That was actually suggested in <a rel="nofollow" href="http://jmlr.csail.mit.edu/papers/volume8/koh07a/koh07a.pdf" target="_blank">Koh, Kim &amp; Boyd (2007)</a>, with some implementation in matlab, see <a rel="nofollow" href="http://web.stanford.edu/~boyd/l1_logreg/" target="_blank">http://web.stanford.edu/~boyd/l1_logreg/</a>. If I can find some time, later one, maybe I will try to recode it. But actually, it is not the technique used in most R functions.</p>
<p>Now, o be honest, we face a double challenge today: the first one is to understand how lasso works for the “standard” (least square) problem, the second one is to see how to adapt it to the logistic case.</p>
<h2>Standard lasso (with weights)</h2>
<p>If we get back to the original Lasso approach, the goal was to solve<img alt="\min\left\lbrace\frac{1}{2n}\sum_{i=1}^n [y_i-(\beta_0+\mathbf{x}_i^T\mathbf{\beta})]^2+\lambda \sum_j |\beta_j|\right\rbrace" title="\min\left\lbrace\frac{1}{2n}\sum_{i=1}^n [y_i-(\beta_0+\mathbf{x}_i^T\mathbf{\beta})]^2+\lambda \sum_j |\beta_j|\right\rbrace" src="https://s0.wp.com/latex.php?latex=%5Cmin%5Cleft%5Clbrace%5Cfrac%7B1%7D%7B2n%7D%5Csum_%7Bi%3D1%7D%5En+%5By_i-%28%5Cbeta_0%2B%5Cmathbf%7Bx%7D_i%5ET%5Cmathbf%7B%5Cbeta%7D%29%5D%5E2%2B%5Clambda+%5Csum_j+%7C%5Cbeta_j%7C%5Cright%5Crbrace&bg=ffffff&fg=000&s=0" class="latex">(with standard notions, as in <a rel="nofollow" href="https://en.wikipedia.org/wiki/Lasso_(statistics)" target="_blank">wikipedia</a> or <a> Jocelyn Chi’s post</a> – most of the code in this section is inspired by Jocelyn’s great post).</p>
<p>Observe that the intercept is not subject to the penalty. The first order condition is then<img alt="\frac{\partial}{\partial\beta_0}\|\mathbf{y}-\mathbf{X}\mathbf{\beta}-\beta_0\mathbf{1}\|^2=(\mathbf{X}\mathbf{\beta}-\mathbf{y})^T\mathbf{1}+\beta_0\|\mathbf{1}\|^2=0" title="\frac{\partial}{\partial\beta_0}\|\mathbf{y}-\mathbf{X}\mathbf{\beta}-\beta_0\mathbf{1}\|^2=(\mathbf{X}\mathbf{\beta}-\mathbf{y})^T\mathbf{1}+\beta_0\|\mathbf{1}\|^2=0" src="https://s0.wp.com/latex.php?latex=%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%5Cbeta_0%7D%5C%7C%5Cmathbf%7By%7D-%5Cmathbf%7BX%7D%5Cmathbf%7B%5Cbeta%7D-%5Cbeta_0%5Cmathbf%7B1%7D%5C%7C%5E2%3D%28%5Cmathbf%7BX%7D%5Cmathbf%7B%5Cbeta%7D-%5Cmathbf%7By%7D%29%5ET%5Cmathbf%7B1%7D%2B%5Cbeta_0%5C%7C%5Cmathbf%7B1%7D%5C%7C%5E2%3D0&bg=ffffff&fg=000&s=0" class="latex">i.e.<img alt="\beta_0=\frac{1}{n^2}(\mathbf{X}\mathbf{\beta}-\mathbf{y})^T\mathbf{1}" title="\beta_0=\frac{1}{n^2}(\mathbf{X}\mathbf{\beta}-\mathbf{y})^T\mathbf{1}" src="https://s0.wp.com/latex.php?latex=%5Cbeta_0%3D%5Cfrac%7B1%7D%7Bn%5E2%7D%28%5Cmathbf%7BX%7D%5Cmathbf%7B%5Cbeta%7D-%5Cmathbf%7By%7D%29%5ET%5Cmathbf%7B1%7D&bg=ffffff&fg=000&s=0" class="latex">Assume now that <a rel="nofollow" href="https://en.wikipedia.org/wiki/Karush%E2%80%93Kuhn%E2%80%93Tucker_conditions" target="_blank">KKT conditions</a> are satisfied, since we cannot differentiate (to find points where the gradient is <img alt="\mathbf{0}" title="\mathbf{0}" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7B0%7D&bg=ffffff&fg=000&s=0" class="latex">), we can check if <img alt="\mathbf{0}" title="\mathbf{0}" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7B0%7D&bg=ffffff&fg=000&s=0" class="latex"> contains the subdifferential at the minimum.<br>
<img src="https://i1.wp.com/f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2018/06/subgradient.png?resize=331%2C155&ssl=1" alt="" width="331" data-recalc-dims="1" height="155"></p>
<p>Namely<img alt="\mathbf{0}\in\partial \left(\frac{1}{2}\|\mathbf{y}-\mathbf{X}\mathbf{\beta}\|^2+\lambda\|\mathbf{\beta}\|_{\ell_1}\right)=\frac{1}{2}\nabla\|\mathbf{y}-\mathbf{X}\mathbf{\beta}\|^2+\partial(\lambda\|\mathbf{\beta}\|_{\ell_1})" title="\mathbf{0}\in\partial \left(\frac{1}{2}\|\mathbf{y}-\mathbf{X}\mathbf{\beta}\|^2+\lambda\|\mathbf{\beta}\|_{\ell_1}\right)=\frac{1}{2}\nabla\|\mathbf{y}-\mathbf{X}\mathbf{\beta}\|^2+\partial(\lambda\|\mathbf{\beta}\|_{\ell_1})" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7B0%7D%5Cin%5Cpartial+%5Cleft%28%5Cfrac%7B1%7D%7B2%7D%5C%7C%5Cmathbf%7By%7D-%5Cmathbf%7BX%7D%5Cmathbf%7B%5Cbeta%7D%5C%7C%5E2%2B%5Clambda%5C%7C%5Cmathbf%7B%5Cbeta%7D%5C%7C_%7B%5Cell_1%7D%5Cright%29%3D%5Cfrac%7B1%7D%7B2%7D%5Cnabla%5C%7C%5Cmathbf%7By%7D-%5Cmathbf%7BX%7D%5Cmathbf%7B%5Cbeta%7D%5C%7C%5E2%2B%5Cpartial%28%5Clambda%5C%7C%5Cmathbf%7B%5Cbeta%7D%5C%7C_%7B%5Cell_1%7D%29&bg=ffffff&fg=000&s=0" class="latex"><br>
For the term on the left, we recognize <img alt="\frac{1}{2}\nabla\|\mathbf{y}-\mathbf{X}\mathbf{\beta}\|^2=-\mathbf{X}^T(\mathbf{y}-\mathbf{X}\mathbf{\beta})=-\mathbf{g}" title="\frac{1}{2}\nabla\|\mathbf{y}-\mathbf{X}\mathbf{\beta}\|^2=-\mathbf{X}^T(\mathbf{y}-\mathbf{X}\mathbf{\beta})=-\mathbf{g}" src="https://s0.wp.com/latex.php?latex=%5Cfrac%7B1%7D%7B2%7D%5Cnabla%5C%7C%5Cmathbf%7By%7D-%5Cmathbf%7BX%7D%5Cmathbf%7B%5Cbeta%7D%5C%7C%5E2%3D-%5Cmathbf%7BX%7D%5ET%28%5Cmathbf%7By%7D-%5Cmathbf%7BX%7D%5Cmathbf%7B%5Cbeta%7D%29%3D-%5Cmathbf%7Bg%7D&bg=ffffff&fg=000&s=0" class="latex">so that the previous equation can be writen<img alt="g_k\in\partial(\lambda|\beta_k|)=\begin{cases}\{+\lambda\}\text{ if }\beta_k&gt;0 \\ \{-\lambda\}\text{ if }\beta_k&lt;0 \\ (-\lambda,+\lambda)\text{ if }\beta_k=0\end{cases}[/latex]i.e. if [latex]\beta_k\neq 0[/latex], then [latex]g_k = \text{sign}(\beta_k)\cdot\lambda[/latex].&lt;/p&gt; &lt;p&gt;Then we write the &lt;a href=&quot;https://en.wikipedia.org/wiki/Karush%E2%80%93Kuhn%E2%80%93Tucker_conditions&quot;&gt;KKT conditions&lt;/a&gt; for this formulation and simplify them to produce a set of rules for checking our solution&lt;/p&gt; &lt;p&gt;We can split [latex]\beta_j" title="g_k\in\partial(\lambda|\beta_k|)=\begin{cases}\{+\lambda\}\text{ if }\beta_k&gt;0 \\ \{-\lambda\}\text{ if }\beta_k&lt;0 \\ (-\lambda,+\lambda)\text{ if }\beta_k=0\end{cases}[/latex]i.e. if [latex]\beta_k\neq 0[/latex], then [latex]g_k = \text{sign}(\beta_k)\cdot\lambda[/latex].&lt;/p&gt; &lt;p&gt;Then we write the &lt;a href=&quot;https://en.wikipedia.org/wiki/Karush%E2%80%93Kuhn%E2%80%93Tucker_conditions&quot;&gt;KKT conditions&lt;/a&gt; for this formulation and simplify them to produce a set of rules for checking our solution&lt;/p&gt; &lt;p&gt;We can split [latex]\beta_j" src="https://s0.wp.com/latex.php?latex=g_k%5Cin%5Cpartial%28%5Clambda%7C%5Cbeta_k%7C%29%3D%5Cbegin%7Bcases%7D%5C%7B%2B%5Clambda%5C%7D%5Ctext%7B+if+%7D%5Cbeta_k%3E0+%5C%5C+%5C%7B-%5Clambda%5C%7D%5Ctext%7B+if+%7D%5Cbeta_k%3C0+%5C%5C+%28-%5Clambda%2C%2B%5Clambda%29%5Ctext%7B+if+%7D%5Cbeta_k%3D0%5Cend%7Bcases%7D%26%2391%3B%2Flatex%26%2393%3Bi.e.+if+%26%2391%3Blatex%26%2393%3B%5Cbeta_k%5Cneq+0%26%2391%3B%2Flatex%26%2393%3B%2C+then+%26%2391%3Blatex%26%2393%3Bg_k+%3D+%5Ctext%7Bsign%7D%28%5Cbeta_k%29%5Ccdot%5Clambda%26%2391%3B%2Flatex%26%2393%3B.%3C%2Fp%3E+%3Cp%3EThen+we+write+the+%3Ca+href%3D%22https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FKarush%25E2%2580%2593Kuhn%25E2%2580%2593Tucker_conditions%22%3EKKT+conditions%3C%2Fa%3E+for+this+formulation+and+simplify+them+to+produce+a+set+of+rules+for+checking+our+solution%3C%2Fp%3E+%3Cp%3EWe+can+split+%5Blatex%5D%5Cbeta_j&bg=ffffff&fg=000&s=0" class="latex"> into a sum of its positive and negative parts by replacing <img alt="\beta_j" title="\beta_j" src="https://s0.wp.com/latex.php?latex=%5Cbeta_j&bg=ffffff&fg=000&s=0" class="latex"> with <img alt="\beta_j^+-\beta_j^-" title="\beta_j^+-\beta_j^-" src="https://s0.wp.com/latex.php?latex=%5Cbeta_j%5E%2B-%5Cbeta_j%5E-&bg=ffffff&fg=000&s=0" class="latex"> where <img alt="\beta_j^+,\beta_j^-\geq0" title="\beta_j^+,\beta_j^-\geq0" src="https://s0.wp.com/latex.php?latex=%5Cbeta_j%5E%2B%2C%5Cbeta_j%5E-%5Cgeq0&bg=ffffff&fg=000&s=0" class="latex">. Then the Lasso problem becomes<img alt="-\log\mathcal{L}(\mathbf{\beta})+\lambda\sum_j(\beta_j^+-\beta_j^-)" title="-\log\mathcal{L}(\mathbf{\beta})+\lambda\sum_j(\beta_j^+-\beta_j^-)" src="https://s0.wp.com/latex.php?latex=-%5Clog%5Cmathcal%7BL%7D%28%5Cmathbf%7B%5Cbeta%7D%29%2B%5Clambda%5Csum_j%28%5Cbeta_j%5E%2B-%5Cbeta_j%5E-%29&bg=ffffff&fg=000&s=0" class="latex">with constraints <img alt="\beta_j^+-\beta_j^-" title="\beta_j^+-\beta_j^-" src="https://s0.wp.com/latex.php?latex=%5Cbeta_j%5E%2B-%5Cbeta_j%5E-&bg=ffffff&fg=000&s=0" class="latex">.</p>
<p>Let <img alt="\alpha_j^+,\alpha_j^-" title="\alpha_j^+,\alpha_j^-" src="https://s0.wp.com/latex.php?latex=%5Calpha_j%5E%2B%2C%5Calpha_j%5E-&bg=ffffff&fg=000&s=0" class="latex"> denote the Lagrange multipliers for <img alt="\beta_j^+,\beta_j^-" title="\beta_j^+,\beta_j^-" src="https://s0.wp.com/latex.php?latex=%5Cbeta_j%5E%2B%2C%5Cbeta_j%5E-&bg=ffffff&fg=000&s=0" class="latex">, respectively.</p>
<p><img alt="L({\mathbf{\beta}}) + \lambda \sum_{j} (\beta_{j}^{+} - \beta_{j}^{-}) - \sum_{j}\alpha_{j}^{+}\beta_{j}^{+} - \sum_{j} \alpha_{j}^{-}\beta_{j}^{-}." title="L({\mathbf{\beta}}) + \lambda \sum_{j} (\beta_{j}^{+} - \beta_{j}^{-}) - \sum_{j}\alpha_{j}^{+}\beta_{j}^{+} - \sum_{j} \alpha_{j}^{-}\beta_{j}^{-}." src="https://s0.wp.com/latex.php?latex=L%28%7B%5Cmathbf%7B%5Cbeta%7D%7D%29+%2B+%5Clambda+%5Csum_%7Bj%7D+%28%5Cbeta_%7Bj%7D%5E%7B%2B%7D+-+%5Cbeta_%7Bj%7D%5E%7B-%7D%29+-+%5Csum_%7Bj%7D%5Calpha_%7Bj%7D%5E%7B%2B%7D%5Cbeta_%7Bj%7D%5E%7B%2B%7D+-+%5Csum_%7Bj%7D+%5Calpha_%7Bj%7D%5E%7B-%7D%5Cbeta_%7Bj%7D%5E%7B-%7D.&bg=ffffff&fg=000&s=0" class="latex">To satisfy the stationarity condition, we take the gradient of the Lagrangian with respect to <img alt="\beta_{j}^{+}" title="\beta_{j}^{+}" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%7D%5E%7B%2B%7D&bg=ffffff&fg=000&s=0" class="latex"> and set it to zero to obtain<img alt="\nabla L({\mathbf{\beta}})_{j} + \lambda - \alpha_{j}^{+} = 0" title="\nabla L({\mathbf{\beta}})_{j} + \lambda - \alpha_{j}^{+} = 0" src="https://s0.wp.com/latex.php?latex=%5Cnabla+L%28%7B%5Cmathbf%7B%5Cbeta%7D%7D%29_%7Bj%7D+%2B+%5Clambda+-+%5Calpha_%7Bj%7D%5E%7B%2B%7D+%3D+0&bg=ffffff&fg=000&s=0" class="latex">We do the same with respect to <img alt="\beta_{j}^{-}" title="\beta_{j}^{-}" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%7D%5E%7B-%7D&bg=ffffff&fg=000&s=0" class="latex"> to obtain<img alt="-\nabla L({\mathbf{\beta}})_{j}+\lambda-\alpha_{j}^{-} = 0" title="-\nabla L({\mathbf{\beta}})_{j}+\lambda-\alpha_{j}^{-} = 0" src="https://s0.wp.com/latex.php?latex=-%5Cnabla+L%28%7B%5Cmathbf%7B%5Cbeta%7D%7D%29_%7Bj%7D%2B%5Clambda-%5Calpha_%7Bj%7D%5E%7B-%7D+%3D+0&bg=ffffff&fg=000&s=0" class="latex"></p>
<p>As discussed in <a> Jocelyn Chi’s post</a>, primal feasibility requires that the primal constraints be satisfied so this gives us <img alt="\beta_{j}^{+} \ge 0" title="\beta_{j}^{+} \ge 0" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%7D%5E%7B%2B%7D+%5Cge+0&bg=ffffff&fg=000&s=0" class="latex"> and <img alt="\beta_{j}^{-} \ge 0" title="\beta_{j}^{-} \ge 0" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%7D%5E%7B-%7D+%5Cge+0&bg=ffffff&fg=000&s=0" class="latex">. Then dual feasibility requires non-negativity of the Lagrange multipliers so we get <img alt="\alpha_{j}^{+} \ge 0" title="\alpha_{j}^{+} \ge 0" src="https://s0.wp.com/latex.php?latex=%5Calpha_%7Bj%7D%5E%7B%2B%7D+%5Cge+0&bg=ffffff&fg=000&s=0" class="latex"> and <img alt="\alpha_{j}^{-} \ge 0" title="\alpha_{j}^{-} \ge 0" src="https://s0.wp.com/latex.php?latex=%5Calpha_%7Bj%7D%5E%7B-%7D+%5Cge+0&bg=ffffff&fg=000&s=0" class="latex">. And finally, complementary slackness requires that <img alt="\alpha_{j}^{+}\beta_{j}^{+} = 0" title="\alpha_{j}^{+}\beta_{j}^{+} = 0" src="https://s0.wp.com/latex.php?latex=%5Calpha_%7Bj%7D%5E%7B%2B%7D%5Cbeta_%7Bj%7D%5E%7B%2B%7D+%3D+0&bg=ffffff&fg=000&s=0" class="latex"> and <img alt="\alpha_{j}^{-}\beta_{j}^{-} = 0" title="\alpha_{j}^{-}\beta_{j}^{-} = 0" src="https://s0.wp.com/latex.php?latex=%5Calpha_%7Bj%7D%5E%7B-%7D%5Cbeta_%7Bj%7D%5E%7B-%7D+%3D+0&bg=ffffff&fg=000&s=0" class="latex">We can simplify these conditions to obtain a simple set of rules for checking whether or not our solution is a minimum.</p>
<p>From <img alt="\nabla L(\beta)_{j} + \lambda - \alpha_{j}^{+} = 0" title="\nabla L(\beta)_{j} + \lambda - \alpha_{j}^{+} = 0" src="https://s0.wp.com/latex.php?latex=%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%2B+%5Clambda+-+%5Calpha_%7Bj%7D%5E%7B%2B%7D+%3D+0&bg=ffffff&fg=000&s=0" class="latex">, we have <img alt="\nabla L(\beta)_{j} + \lambda= \alpha_{j}^{+} \ge 0" title="\nabla L(\beta)_{j} + \lambda= \alpha_{j}^{+} \ge 0" src="https://s0.wp.com/latex.php?latex=%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%2B+%5Clambda%3D+%5Calpha_%7Bj%7D%5E%7B%2B%7D+%5Cge+0&bg=ffffff&fg=000&s=0" class="latex">. This gives us <img alt="\nabla L(\beta)_{j} \ge -\lambda" title="\nabla L(\beta)_{j} \ge -\lambda" src="https://s0.wp.com/latex.php?latex=%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%5Cge+-%5Clambda&bg=ffffff&fg=000&s=0" class="latex">. From <img alt="-\nabla L(\beta)_{j} + \lambda - \alpha_{j}^{-} = 0" title="-\nabla L(\beta)_{j} + \lambda - \alpha_{j}^{-} = 0" src="https://s0.wp.com/latex.php?latex=-%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%2B+%5Clambda+-+%5Calpha_%7Bj%7D%5E%7B-%7D+%3D+0&bg=ffffff&fg=000&s=0" class="latex">, we have <img alt="-\nabla L(\beta)_{j} + \lambda = \alpha_{j}^{-} \ge 0" title="-\nabla L(\beta)_{j} + \lambda = \alpha_{j}^{-} \ge 0" src="https://s0.wp.com/latex.php?latex=-%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%2B+%5Clambda+%3D+%5Calpha_%7Bj%7D%5E%7B-%7D+%5Cge+0&bg=ffffff&fg=000&s=0" class="latex">. This gives us <img alt="-\nabla L(\beta)_{j} \ge -\lambda" title="-\nabla L(\beta)_{j} \ge -\lambda" src="https://s0.wp.com/latex.php?latex=-%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%5Cge+-%5Clambda&bg=ffffff&fg=000&s=0" class="latex">, which gives us <img alt="\nabla L(\beta)_{j} \le \lambda" title="\nabla L(\beta)_{j} \le \lambda" src="https://s0.wp.com/latex.php?latex=%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%5Cle+%5Clambda&bg=ffffff&fg=000&s=0" class="latex">. Hence, <img alt="\lvert \nabla L(\beta)_{j} \rvert \le \lambda \; \forall j" title="\lvert \nabla L(\beta)_{j} \rvert \le \lambda \; \forall j" src="https://s0.wp.com/latex.php?latex=%5Clvert+%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%5Crvert+%5Cle+%5Clambda+%5C%3B+%5Cforall+j&bg=ffffff&fg=000&s=0" class="latex"></p>
<p>When <img alt="\beta_{j}^{+} &gt; 0, \lambda &gt; 0" title="\beta_{j}^{+} &gt; 0, \lambda &gt; 0" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%7D%5E%7B%2B%7D+%3E+0%2C+%5Clambda+%3E+0&bg=ffffff&fg=000&s=0" class="latex">, complementary slackness requires <img alt="\alpha_{j}^{+} = 0" title="\alpha_{j}^{+} = 0" src="https://s0.wp.com/latex.php?latex=%5Calpha_%7Bj%7D%5E%7B%2B%7D+%3D+0&bg=ffffff&fg=000&s=0" class="latex">. So <img alt="\nabla L(\beta)_{j} + \lambda = \alpha_{j}^{+} = 0" title="\nabla L(\beta)_{j} + \lambda = \alpha_{j}^{+} = 0" src="https://s0.wp.com/latex.php?latex=%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%2B+%5Clambda+%3D+%5Calpha_%7Bj%7D%5E%7B%2B%7D+%3D+0&bg=ffffff&fg=000&s=0" class="latex">. Hence, <img alt="\nabla L(\beta)_{j} = -\lambda &lt; 0[/latex] since [latex]\lambda &gt; 0" title="\nabla L(\beta)_{j} = -\lambda &lt; 0[/latex] since [latex]\lambda &gt; 0" src="https://s0.wp.com/latex.php?latex=%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%3D+-%5Clambda+%3C+0%26%2391%3B%2Flatex%26%2393%3B+since+%26%2391%3Blatex%26%2393%3B%5Clambda+%3E+0&bg=ffffff&fg=000&s=0" class="latex">. At the same time, <img alt="-\nabla L(\beta)_{j} + \lambda = \alpha_{j}^{-} \ge 0" title="-\nabla L(\beta)_{j} + \lambda = \alpha_{j}^{-} \ge 0" src="https://s0.wp.com/latex.php?latex=-%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%2B+%5Clambda+%3D+%5Calpha_%7Bj%7D%5E%7B-%7D+%5Cge+0&bg=ffffff&fg=000&s=0" class="latex"> so <img alt="2 \lambda = \alpha_{j}^{-} &gt; 0" title="2 \lambda = \alpha_{j}^{-} &gt; 0" src="https://s0.wp.com/latex.php?latex=2+%5Clambda+%3D+%5Calpha_%7Bj%7D%5E%7B-%7D+%3E+0&bg=ffffff&fg=000&s=0" class="latex"> since <img alt="\lambda &gt; 0" title="\lambda &gt; 0" src="https://s0.wp.com/latex.php?latex=%5Clambda+%3E+0&bg=ffffff&fg=000&s=0" class="latex">. Then complementary slackness requires <img alt="\beta_{j}^{-} = 0" title="\beta_{j}^{-} = 0" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%7D%5E%7B-%7D+%3D+0&bg=ffffff&fg=000&s=0" class="latex">. Hence, when <img alt="\beta_{j}^{+} &gt; 0" title="\beta_{j}^{+} &gt; 0" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%7D%5E%7B%2B%7D+%3E+0&bg=ffffff&fg=000&s=0" class="latex">, we have <img alt="\beta_{j}^{-}=0" title="\beta_{j}^{-}=0" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%7D%5E%7B-%7D%3D0&bg=ffffff&fg=000&s=0" class="latex"> and <img alt="\nabla L(\beta)_{j} = -\lambda" title="\nabla L(\beta)_{j} = -\lambda" src="https://s0.wp.com/latex.php?latex=%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%3D+-%5Clambda&bg=ffffff&fg=000&s=0" class="latex"></p>
<p>Similarly, when <img alt="\beta_{j}^{-} &gt; 0, \lambda &gt; 0" title="\beta_{j}^{-} &gt; 0, \lambda &gt; 0" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%7D%5E%7B-%7D+%3E+0%2C+%5Clambda+%3E+0&bg=ffffff&fg=000&s=0" class="latex">, complementary slackness requires <img alt="\alpha_{j}^{-}=0" title="\alpha_{j}^{-}=0" src="https://s0.wp.com/latex.php?latex=%5Calpha_%7Bj%7D%5E%7B-%7D%3D0&bg=ffffff&fg=000&s=0" class="latex">. So <img alt="-\nabla L(\beta)_{j} + \lambda = \alpha_{j}^{-} = 0" title="-\nabla L(\beta)_{j} + \lambda = \alpha_{j}^{-} = 0" src="https://s0.wp.com/latex.php?latex=-%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%2B+%5Clambda+%3D+%5Calpha_%7Bj%7D%5E%7B-%7D+%3D+0&bg=ffffff&fg=000&s=0" class="latex"> and <img alt="\nabla L(\beta)_{j}=\lambda&gt;0" title="\nabla L(\beta)_{j}=\lambda&gt;0" src="https://s0.wp.com/latex.php?latex=%5Cnabla+L%28%5Cbeta%29_%7Bj%7D%3D%5Clambda%3E0&bg=ffffff&fg=000&s=0" class="latex"> since <img alt="\lambda &gt; 0" title="\lambda &gt; 0" src="https://s0.wp.com/latex.php?latex=%5Clambda+%3E+0&bg=ffffff&fg=000&s=0" class="latex">. Then from <img alt="\nabla L(\beta)_{j} + \lambda = \alpha_{j}^{+} \ge 0" title="\nabla L(\beta)_{j} + \lambda = \alpha_{j}^{+} \ge 0" src="https://s0.wp.com/latex.php?latex=%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%2B+%5Clambda+%3D+%5Calpha_%7Bj%7D%5E%7B%2B%7D+%5Cge+0&bg=ffffff&fg=000&s=0" class="latex"> and the above, we get <img alt="2 \lambda = \alpha_{j}^{+} &gt; 0" title="2 \lambda = \alpha_{j}^{+} &gt; 0" src="https://s0.wp.com/latex.php?latex=2+%5Clambda+%3D+%5Calpha_%7Bj%7D%5E%7B%2B%7D+%3E+0&bg=ffffff&fg=000&s=0" class="latex">. Then complementary slackness requires <img alt="\beta_{j}^{+} = 0" title="\beta_{j}^{+} = 0" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%7D%5E%7B%2B%7D+%3D+0&bg=ffffff&fg=000&s=0" class="latex">. Hence, when <img alt="\beta_{j}^{-} &gt; 0" title="\beta_{j}^{-} &gt; 0" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%7D%5E%7B-%7D+%3E+0&bg=ffffff&fg=000&s=0" class="latex">, we have <img alt="\beta_{j}^{+}=0" title="\beta_{j}^{+}=0" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%7D%5E%7B%2B%7D%3D0&bg=ffffff&fg=000&s=0" class="latex"> and <img alt="\nabla L(\beta)_{j} = \lambda" title="\nabla L(\beta)_{j} = \lambda" src="https://s0.wp.com/latex.php?latex=%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%3D+%5Clambda&bg=ffffff&fg=000&s=0" class="latex">.</p>
<p>Since <img alt="\beta_{j} = \beta_{j}^{+} - \beta_{j}^{-}" title="\beta_{j} = \beta_{j}^{+} - \beta_{j}^{-}" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%7D+%3D+%5Cbeta_%7Bj%7D%5E%7B%2B%7D+-+%5Cbeta_%7Bj%7D%5E%7B-%7D&bg=ffffff&fg=000&s=0" class="latex">, this means that when <img alt="\beta_{j} &gt; 0" title="\beta_{j} &gt; 0" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%7D+%3E+0&bg=ffffff&fg=000&s=0" class="latex">, <img alt="\nabla L(\beta)_{j} = -\lambda" title="\nabla L(\beta)_{j} = -\lambda" src="https://s0.wp.com/latex.php?latex=%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%3D+-%5Clambda&bg=ffffff&fg=000&s=0" class="latex">. And when <img alt="\beta_{j} &lt;0[/latex], [latex]\nabla L(\beta)_{j} = \lambda[/latex]. Combining this with [latex]\lvert \nabla L(\beta)_{j} \rvert \le \lambda \; \forall j[/latex], we arrive at the same convergence requirements that we obtained before using subdifferential calculus.&lt;/p&gt; &lt;p&gt;For conveniency, introduce the soft-thresholding function[latex display=&quot;true&quot;]S(z,\gamma)=\text{sign}(z)\cdot(|z|-\gamma)_+=\begin{cases}z-\gamma&amp;\text{ if }\gamma&gt;|z|\text{ and }z&lt;0\\z+\gamma&amp;\text{ if }\gamma&lt;|z|\text{ and }z&lt;0 \\0&amp;\text{ if }\gamma\geq|z|\end{cases}[/latex]&lt;br /&gt; Noticing that the optimization problem [latex display=&quot;true&quot;]\frac{1}{2}\|\mathbf{y}-\mathbf{X}\mathbf{\beta}\|_{\ell_2}^2+\lambda\|\mathbf{\beta}\|_{\ell_1}" title="\beta_{j} &lt;0[/latex], [latex]\nabla L(\beta)_{j} = \lambda[/latex]. Combining this with [latex]\lvert \nabla L(\beta)_{j} \rvert \le \lambda \; \forall j[/latex], we arrive at the same convergence requirements that we obtained before using subdifferential calculus.&lt;/p&gt; &lt;p&gt;For conveniency, introduce the soft-thresholding function[latex display=&quot;true&quot;]S(z,\gamma)=\text{sign}(z)\cdot(|z|-\gamma)_+=\begin{cases}z-\gamma&amp;\text{ if }\gamma&gt;|z|\text{ and }z&lt;0\\z+\gamma&amp;\text{ if }\gamma&lt;|z|\text{ and }z&lt;0 \\0&amp;\text{ if }\gamma\geq|z|\end{cases}[/latex]&lt;br /&gt; Noticing that the optimization problem [latex display=&quot;true&quot;]\frac{1}{2}\|\mathbf{y}-\mathbf{X}\mathbf{\beta}\|_{\ell_2}^2+\lambda\|\mathbf{\beta}\|_{\ell_1}" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%7D+%3C0%26%2391%3B%2Flatex%26%2393%3B%2C+%26%2391%3Blatex%26%2393%3B%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%3D+%5Clambda%26%2391%3B%2Flatex%26%2393%3B.+Combining+this+with+%26%2391%3Blatex%26%2393%3B%5Clvert+%5Cnabla+L%28%5Cbeta%29_%7Bj%7D+%5Crvert+%5Cle+%5Clambda+%5C%3B+%5Cforall+j%26%2391%3B%2Flatex%26%2393%3B%2C+we+arrive+at+the+same+convergence+requirements+that+we+obtained+before+using+subdifferential+calculus.%3C%2Fp%3E+%3Cp%3EFor+conveniency%2C+introduce+the+soft-thresholding+function%5Blatex+display%3D%22true%22%5DS%28z%2C%5Cgamma%29%3D%5Ctext%7Bsign%7D%28z%29%5Ccdot%28%7Cz%7C-%5Cgamma%29_%2B%3D%5Cbegin%7Bcases%7Dz-%5Cgamma%26%5Ctext%7B+if+%7D%5Cgamma%3E%7Cz%7C%5Ctext%7B+and+%7Dz%3C0%5C%5Cz%2B%5Cgamma%26%5Ctext%7B+if+%7D%5Cgamma%3C%7Cz%7C%5Ctext%7B+and+%7Dz%3C0+%5C%5C0%26%5Ctext%7B+if+%7D%5Cgamma%5Cgeq%7Cz%7C%5Cend%7Bcases%7D%26%2391%3B%2Flatex%26%2393%3B%3Cbr+%2F%3E+Noticing+that+the+optimization+problem+%5Blatex+display%3D%22true%22%5D%5Cfrac%7B1%7D%7B2%7D%5C%7C%5Cmathbf%7By%7D-%5Cmathbf%7BX%7D%5Cmathbf%7B%5Cbeta%7D%5C%7C_%7B%5Cell_2%7D%5E2%2B%5Clambda%5C%7C%5Cmathbf%7B%5Cbeta%7D%5C%7C_%7B%5Cell_1%7D&bg=ffffff&fg=000&s=0" class="latex">can also be written<br>
<img alt="\min\left\lbrace\sum_{j=1}^p -\widehat{\beta}_j^{ols}\cdot\beta_j+\frac{1}{2}\beta_j^2+\lambda|\beta_j|\right\rbrace" title="\min\left\lbrace\sum_{j=1}^p -\widehat{\beta}_j^{ols}\cdot\beta_j+\frac{1}{2}\beta_j^2+\lambda|\beta_j|\right\rbrace" src="https://s0.wp.com/latex.php?latex=%5Cmin%5Cleft%5Clbrace%5Csum_%7Bj%3D1%7D%5Ep+-%5Cwidehat%7B%5Cbeta%7D_j%5E%7Bols%7D%5Ccdot%5Cbeta_j%2B%5Cfrac%7B1%7D%7B2%7D%5Cbeta_j%5E2%2B%5Clambda%7C%5Cbeta_j%7C%5Cright%5Crbrace&bg=ffffff&fg=000&s=0" class="latex">observe that<img alt="\widehat{\beta}_{j,\lambda}=S(\widehat{\beta}_j^{ols},\lambda)" title="\widehat{\beta}_{j,\lambda}=S(\widehat{\beta}_j^{ols},\lambda)" src="https://s0.wp.com/latex.php?latex=%5Cwidehat%7B%5Cbeta%7D_%7Bj%2C%5Clambda%7D%3DS%28%5Cwidehat%7B%5Cbeta%7D_j%5E%7Bols%7D%2C%5Clambda%29&bg=ffffff&fg=000&s=0" class="latex">which is a coordinate-wise update.</p>
<p>Now, if we consider a (slightly) more general problem, with weights in the first part<img alt="\min\left\lbrace\frac{1}{2n}\sum_{i=1}^n{\color{red}{\omega_i}} [y_i-(\beta_0+\mathbf{x}_i^T\mathbf{\beta})]^2+\lambda \sum_j |\beta_j|\right\rbrace" title="\min\left\lbrace\frac{1}{2n}\sum_{i=1}^n{\color{red}{\omega_i}} [y_i-(\beta_0+\mathbf{x}_i^T\mathbf{\beta})]^2+\lambda \sum_j |\beta_j|\right\rbrace" src="https://s0.wp.com/latex.php?latex=%5Cmin%5Cleft%5Clbrace%5Cfrac%7B1%7D%7B2n%7D%5Csum_%7Bi%3D1%7D%5En%7B%5Ccolor%7Bred%7D%7B%5Comega_i%7D%7D+%5By_i-%28%5Cbeta_0%2B%5Cmathbf%7Bx%7D_i%5ET%5Cmathbf%7B%5Cbeta%7D%29%5D%5E2%2B%5Clambda+%5Csum_j+%7C%5Cbeta_j%7C%5Cright%5Crbrace&bg=ffffff&fg=000&s=0" class="latex">the coordinate-wise update becomes<br>
<img alt="\widehat{\beta}_{j,\lambda,{\color{red}{\omega}}}=S(\widehat{\beta}_j^{{\color{red}{\omega-}}ols},\lambda)" title="\widehat{\beta}_{j,\lambda,{\color{red}{\omega}}}=S(\widehat{\beta}_j^{{\color{red}{\omega-}}ols},\lambda)" src="https://s0.wp.com/latex.php?latex=%5Cwidehat%7B%5Cbeta%7D_%7Bj%2C%5Clambda%2C%7B%5Ccolor%7Bred%7D%7B%5Comega%7D%7D%7D%3DS%28%5Cwidehat%7B%5Cbeta%7D_j%5E%7B%7B%5Ccolor%7Bred%7D%7B%5Comega-%7D%7Dols%7D%2C%5Clambda%29&bg=ffffff&fg=000&s=0" class="latex"><br>
An alternative is to set<img alt="\mathbf{r}_j=\mathbf{y} - \left(\beta_0\mathbf{1}+\sum_{k\neq j}\beta_k\mathbf{x}_k\right)=\mathbf{y}-\widehat{\mathbf{y}}^{(j)}" title="\mathbf{r}_j=\mathbf{y} - \left(\beta_0\mathbf{1}+\sum_{k\neq j}\beta_k\mathbf{x}_k\right)=\mathbf{y}-\widehat{\mathbf{y}}^{(j)}" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7Br%7D_j%3D%5Cmathbf%7By%7D+-+%5Cleft%28%5Cbeta_0%5Cmathbf%7B1%7D%2B%5Csum_%7Bk%5Cneq+j%7D%5Cbeta_k%5Cmathbf%7Bx%7D_k%5Cright%29%3D%5Cmathbf%7By%7D-%5Cwidehat%7B%5Cmathbf%7By%7D%7D%5E%7B%28j%29%7D&bg=ffffff&fg=000&s=0" class="latex"><br>
so that the optimization problem can be written, equivalently<br>
<img alt="\min\left\lbrace\frac{1}{2n}\sum_{j=1}^p [\mathbf{r}_j-\beta_j\mathbf{x}_j]^2+\lambda |\beta_j|\right\rbrace" title="\min\left\lbrace\frac{1}{2n}\sum_{j=1}^p [\mathbf{r}_j-\beta_j\mathbf{x}_j]^2+\lambda |\beta_j|\right\rbrace" src="https://s0.wp.com/latex.php?latex=%5Cmin%5Cleft%5Clbrace%5Cfrac%7B1%7D%7B2n%7D%5Csum_%7Bj%3D1%7D%5Ep+%5B%5Cmathbf%7Br%7D_j-%5Cbeta_j%5Cmathbf%7Bx%7D_j%5D%5E2%2B%5Clambda+%7C%5Cbeta_j%7C%5Cright%5Crbrace&bg=ffffff&fg=000&s=0" class="latex"><br>
hence<img alt="\min\left\lbrace\frac{1}{2n}\sum_{j=1}^p \beta_j^2\|\mathbf{x}_j\|-2\beta_j\mathbf{r}_j^T\mathbf{x}_j+\lambda |\beta_j|\right\rbrace" title="\min\left\lbrace\frac{1}{2n}\sum_{j=1}^p \beta_j^2\|\mathbf{x}_j\|-2\beta_j\mathbf{r}_j^T\mathbf{x}_j+\lambda |\beta_j|\right\rbrace" src="https://s0.wp.com/latex.php?latex=%5Cmin%5Cleft%5Clbrace%5Cfrac%7B1%7D%7B2n%7D%5Csum_%7Bj%3D1%7D%5Ep+%5Cbeta_j%5E2%5C%7C%5Cmathbf%7Bx%7D_j%5C%7C-2%5Cbeta_j%5Cmathbf%7Br%7D_j%5ET%5Cmathbf%7Bx%7D_j%2B%5Clambda+%7C%5Cbeta_j%7C%5Cright%5Crbrace&bg=ffffff&fg=000&s=0" class="latex"><br>
and one gets<br>
<img alt="\beta_{j,\lambda} = \frac{1}{\|\mathbf{x}_j\|^2}S(\mathbf{r}_j^T\mathbf{x}_j,n\lambda)" title="\beta_{j,\lambda} = \frac{1}{\|\mathbf{x}_j\|^2}S(\mathbf{r}_j^T\mathbf{x}_j,n\lambda)" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%2C%5Clambda%7D+%3D+%5Cfrac%7B1%7D%7B%5C%7C%5Cmathbf%7Bx%7D_j%5C%7C%5E2%7DS%28%5Cmathbf%7Br%7D_j%5ET%5Cmathbf%7Bx%7D_j%2Cn%5Clambda%29&bg=ffffff&fg=000&s=0" class="latex"><br>
or, if we develop<br>
<img alt="\beta_{j,\lambda} = \frac{1}{\sum_i x_{ij}^2}S\left(\sum_ix_{i,j}[y_i-\widehat{y}_i^{(j)}],n\lambda\right)" title="\beta_{j,\lambda} = \frac{1}{\sum_i x_{ij}^2}S\left(\sum_ix_{i,j}[y_i-\widehat{y}_i^{(j)}],n\lambda\right)" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%2C%5Clambda%7D+%3D+%5Cfrac%7B1%7D%7B%5Csum_i+x_%7Bij%7D%5E2%7DS%5Cleft%28%5Csum_ix_%7Bi%2Cj%7D%5By_i-%5Cwidehat%7By%7D_i%5E%7B%28j%29%7D%5D%2Cn%5Clambda%5Cright%29&bg=ffffff&fg=000&s=0" class="latex"><br>
Again, if there are weights <img alt="\mathbf{\omega}=(\omega_i)" title="\mathbf{\omega}=(\omega_i)" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7B%5Comega%7D%3D%28%5Comega_i%29&bg=ffffff&fg=000&s=0" class="latex">, the coordinate-wise update becomes<br>
<img alt="\beta_{j,\lambda,{\color{red}{\omega}}} = \frac{1}{\sum_i {\color{red}{\omega_i}}x_{ij}^2}S\left(\sum_i{\color{red}{\omega_i}}x_{i,j}[y_i-\widehat{y}_i^{(j)}],n\lambda\right)" title="\beta_{j,\lambda,{\color{red}{\omega}}} = \frac{1}{\sum_i {\color{red}{\omega_i}}x_{ij}^2}S\left(\sum_i{\color{red}{\omega_i}}x_{i,j}[y_i-\widehat{y}_i^{(j)}],n\lambda\right)" src="https://s0.wp.com/latex.php?latex=%5Cbeta_%7Bj%2C%5Clambda%2C%7B%5Ccolor%7Bred%7D%7B%5Comega%7D%7D%7D+%3D+%5Cfrac%7B1%7D%7B%5Csum_i+%7B%5Ccolor%7Bred%7D%7B%5Comega_i%7D%7Dx_%7Bij%7D%5E2%7DS%5Cleft%28%5Csum_i%7B%5Ccolor%7Bred%7D%7B%5Comega_i%7D%7Dx_%7Bi%2Cj%7D%5By_i-%5Cwidehat%7By%7D_i%5E%7B%28j%29%7D%5D%2Cn%5Clambda%5Cright%29&bg=ffffff&fg=000&s=0" class="latex"><br>
The code to compute this componentwise descent is</p>
<div><table><tr><td><pre>1
2
3
4
</pre></td><td><pre>soft_thresholding <span>=</span> <span>function</span><span>(</span>x,a<span>)</span><span>{</span>
  result  a<span>)</span><span>]</span>  a<span>)</span><span>]</span> <span>-</span> a
  result<span>[</span><span>which</span><span>(</span>x <span>&amp;</span>lt<span>;</span> <span>-</span>a<span>)</span><span>]</span> <span>=</span> x<span>[</span><span>which</span><span>(</span>x <span>&amp;</span>lt<span>;</span> <span>-</span>a<span>)</span><span>]</span> <span>+</span> a
  <span>return</span><span>(</span>result<span>)</span><span>}</span></pre></td></tr></table></div>
<p>and the code</p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
</pre></td><td><pre>lasso_coord_desc <span>=</span> <span>function</span><span>(</span>X,y,<span>beta</span>,lambda,tol<span>=</span>1e<span>-</span>6,maxiter<span>=</span><span>1000</span><span>)</span><span>{</span>
  <span>beta</span> <span>=</span> <span>as.<span>matrix</span></span><span>(</span><span>beta</span><span>)</span>
  X <span>=</span> <span>as.<span>matrix</span></span><span>(</span>X<span>)</span>
  omega <span>=</span> <span>rep</span><span>(</span><span>1</span><span>/</span><span>length</span><span>(</span>y<span>)</span>,<span>length</span><span>(</span>y<span>)</span><span>)</span>
  obj <span>=</span> <span>numeric</span><span>(</span><span>length</span><span>=</span><span>(</span>maxiter<span>+</span><span>1</span><span>)</span><span>)</span>
  betalist <span>=</span> <span>list</span><span>(</span><span>length</span><span>(</span>maxiter<span>+</span><span>1</span><span>)</span><span>)</span>
  betalist<span>[</span><span>[</span><span>1</span><span>]</span><span>]</span> <span>=</span> <span>beta</span>
    beta0list <span>=</span> <span>numeric</span><span>(</span><span>length</span><span>(</span>maxiter<span>+</span><span>1</span><span>)</span><span>)</span>
    beta0 <span>=</span> <span>sum</span><span>(</span>y<span>-</span>X<span>%*%</span><span>beta</span><span>)</span><span>/</span><span>(</span><span>length</span><span>(</span>y<span>)</span><span>)</span>
    beta0list<span>[</span><span>1</span><span>]</span> <span>=</span> beta0
    <span>for</span> <span>(</span>j <span>in</span> <span>1</span><span>:</span>maxiter<span>)</span><span>{</span>
      <span>for</span> <span>(</span>k <span>in</span> <span>1</span><span>:</span><span>length</span><span>(</span><span>beta</span><span>)</span><span>)</span><span>{</span>
        r <span>=</span> y <span>-</span> X<span>[</span>,<span>-</span>k<span>]</span><span>%*%</span><span>beta</span><span>[</span><span>-</span>k<span>]</span> <span>-</span> beta0<span>*</span><span>rep</span><span>(</span><span>1</span>,<span>length</span><span>(</span>y<span>)</span><span>)</span>
        <span>beta</span><span>[</span>k<span>]</span> <span>=</span> <span>(</span><span>1</span><span>/</span><span>sum</span><span>(</span>omega<span>*</span>X<span>[</span>,k<span>]</span><span>^</span><span>2</span><span>)</span><span>)</span><span>*</span>soft_thresholding<span>(</span><span>t</span><span>(</span>omega<span>*</span>r<span>)</span><span>%*%</span>X<span>[</span>,k<span>]</span>,<span>length</span><span>(</span>y<span>)</span><span>*</span>lambda<span>)</span>
      <span>}</span>
      beta0 <span>=</span> <span>sum</span><span>(</span>y<span>-</span>X<span>%*%</span><span>beta</span><span>)</span><span>/</span><span>(</span><span>length</span><span>(</span>y<span>)</span><span>)</span>
      beta0list<span>[</span>j<span>+</span><span>1</span><span>]</span> <span>=</span> beta0
      betalist<span>[</span><span>[</span>j<span>+</span><span>1</span><span>]</span><span>]</span> <span>=</span> <span>beta</span>
      obj<span>[</span>j<span>]</span> <span>=</span> <span>(</span><span>1</span><span>/</span><span>2</span><span>)</span><span>*</span><span>(</span><span>1</span><span>/</span><span>length</span><span>(</span>y<span>)</span><span>)</span><span>*</span>norm<span>(</span>omega<span>*</span><span>(</span>y <span>-</span> X<span>%*%</span><span>beta</span> <span>-</span> 
beta0<span>*</span><span>rep</span><span>(</span><span>1</span>,<span>length</span><span>(</span>y<span>)</span><span>)</span><span>)</span>,<span>'F'</span><span>)</span><span>^</span><span>2</span> <span>+</span> lambda<span>*</span><span>sum</span><span>(</span><span>abs</span><span>(</span><span>beta</span><span>)</span><span>)</span>
      <span>if</span> <span>(</span>norm<span>(</span><span>rbind</span><span>(</span>beta0list<span>[</span>j<span>]</span>,betalist<span>[</span><span>[</span>j<span>]</span><span>]</span><span>)</span> <span>-</span> <span>rbind</span><span>(</span>beta0,<span>beta</span><span>)</span>,<span>'F'</span><span>)</span> <span>&amp;</span>lt<span>;</span> tol<span>)</span> <span>{</span> <span>break</span> <span>}</span> 
    <span>}</span> 
<span>return</span><span>(</span><span>list</span><span>(</span>obj<span>=</span>obj<span>[</span><span>1</span><span>:</span>j<span>]</span>,<span>beta</span><span>=</span><span>beta</span>,intercept<span>=</span>beta0<span>)</span><span>)</span> <span>}</span></pre></td></tr></table></div>
<p>Let’s keep that one warm, and let’s get back to our initial problem.</p>
<h2>The lasso logistic regression</h2>
<p>The trick here is that the logistic problem can be formulated as a quadratic programming problem. Recall that the log-likelihood is here <img alt="\log\mathcal{L}=\frac{1}{n}\sum_{i=1}^n y_i\cdot(\beta_0+\mathbf{x}_i^T\mathbf{\beta})-\log[1+\exp(\beta_0+\mathbf{x}_i^T\mathbf{\beta})]" title="\log\mathcal{L}=\frac{1}{n}\sum_{i=1}^n y_i\cdot(\beta_0+\mathbf{x}_i^T\mathbf{\beta})-\log[1+\exp(\beta_0+\mathbf{x}_i^T\mathbf{\beta})]" src="https://s0.wp.com/latex.php?latex=%5Clog%5Cmathcal%7BL%7D%3D%5Cfrac%7B1%7D%7Bn%7D%5Csum_%7Bi%3D1%7D%5En+y_i%5Ccdot%28%5Cbeta_0%2B%5Cmathbf%7Bx%7D_i%5ET%5Cmathbf%7B%5Cbeta%7D%29-%5Clog%5B1%2B%5Cexp%28%5Cbeta_0%2B%5Cmathbf%7Bx%7D_i%5ET%5Cmathbf%7B%5Cbeta%7D%29%5D&bg=ffffff&fg=000&s=0" class="latex"><br>
which is a concave function of the parameters. Hence, one can use a quadratic approximation of the log-likelihood – using Taylor expansion,<img alt="\log\mathcal{L}\approx\log\mathcal{L}'=\frac{1}{n}\sum_{i=1}^n \omega_i\cdot[z_i-(\beta_0+\mathbf{x}_i^T\mathbf{\beta})]^2" title="\log\mathcal{L}\approx\log\mathcal{L}'=\frac{1}{n}\sum_{i=1}^n \omega_i\cdot[z_i-(\beta_0+\mathbf{x}_i^T\mathbf{\beta})]^2" src="https://s0.wp.com/latex.php?latex=%5Clog%5Cmathcal%7BL%7D%5Capprox%5Clog%5Cmathcal%7BL%7D%27%3D%5Cfrac%7B1%7D%7Bn%7D%5Csum_%7Bi%3D1%7D%5En+%5Comega_i%5Ccdot%5Bz_i-%28%5Cbeta_0%2B%5Cmathbf%7Bx%7D_i%5ET%5Cmathbf%7B%5Cbeta%7D%29%5D%5E2&bg=ffffff&fg=000&s=0" class="latex"><br>
where <img alt="z_i" title="z_i" src="https://s0.wp.com/latex.php?latex=z_i&bg=ffffff&fg=000&s=0" class="latex"> is the working response<br>
<img alt="z_i=(\beta_0+\mathbf{x}_i^T\mathbf{\beta})+\frac{y_i-p_i}{p_i[1-p_i]}" title="z_i=(\beta_0+\mathbf{x}_i^T\mathbf{\beta})+\frac{y_i-p_i}{p_i[1-p_i]}" src="https://s0.wp.com/latex.php?latex=z_i%3D%28%5Cbeta_0%2B%5Cmathbf%7Bx%7D_i%5ET%5Cmathbf%7B%5Cbeta%7D%29%2B%5Cfrac%7By_i-p_i%7D%7Bp_i%5B1-p_i%5D%7D&bg=ffffff&fg=000&s=0" class="latex"><br>
<img alt="p_i" title="p_i" src="https://s0.wp.com/latex.php?latex=p_i&bg=ffffff&fg=000&s=0" class="latex"> is the prediction<img alt="p_i = \frac{\exp[\beta_0+\mathbf{x}_i^T\mathbf{\beta}]}{1+\exp[\beta_0+\mathbf{x}_i^T\mathbf{\beta}]}" title="p_i = \frac{\exp[\beta_0+\mathbf{x}_i^T\mathbf{\beta}]}{1+\exp[\beta_0+\mathbf{x}_i^T\mathbf{\beta}]}" src="https://s0.wp.com/latex.php?latex=p_i+%3D+%5Cfrac%7B%5Cexp%5B%5Cbeta_0%2B%5Cmathbf%7Bx%7D_i%5ET%5Cmathbf%7B%5Cbeta%7D%5D%7D%7B1%2B%5Cexp%5B%5Cbeta_0%2B%5Cmathbf%7Bx%7D_i%5ET%5Cmathbf%7B%5Cbeta%7D%5D%7D&bg=ffffff&fg=000&s=0" class="latex">and <img alt="\omega_i" title="\omega_i" src="https://s0.wp.com/latex.php?latex=%5Comega_i&bg=ffffff&fg=000&s=0" class="latex"> are weights <img alt="\omega_i = p_i[1-p_i]" title="\omega_i = p_i[1-p_i]" src="https://s0.wp.com/latex.php?latex=%5Comega_i+%3D+p_i%5B1-p_i%5D&bg=ffffff&fg=000&s=0" class="latex">.</p>
<p>Thus, we obtain a penalized least-square problem. And we can use what was done previously</p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
</pre></td><td><pre>lasso_coord_desc <span>=</span> <span>function</span><span>(</span>X,y,<span>beta</span>,lambda,tol<span>=</span>1e<span>-</span>6,maxiter<span>=</span><span>1000</span><span>)</span><span>{</span>
  <span>beta</span> <span>=</span> <span>as.<span>matrix</span></span><span>(</span><span>beta</span><span>)</span>
  X <span>=</span> <span>as.<span>matrix</span></span><span>(</span>X<span>)</span>
  obj <span>=</span> <span>numeric</span><span>(</span><span>length</span><span>=</span><span>(</span>maxiter<span>+</span><span>1</span><span>)</span><span>)</span>
  betalist <span>=</span> <span>list</span><span>(</span><span>length</span><span>(</span>maxiter<span>+</span><span>1</span><span>)</span><span>)</span>
  betalist<span>[</span><span>[</span><span>1</span><span>]</span><span>]</span> <span>=</span> <span>beta</span>
  beta0 <span>=</span> <span>sum</span><span>(</span>y<span>-</span>X<span>%*%</span><span>beta</span><span>)</span><span>/</span><span>(</span><span>length</span><span>(</span>y<span>)</span><span>)</span>
  p <span>=</span> <span>exp</span><span>(</span>beta0<span>*</span><span>rep</span><span>(</span><span>1</span>,<span>length</span><span>(</span>y<span>)</span><span>)</span> <span>+</span> X<span>%*%</span><span>beta</span><span>)</span><span>/</span><span>(</span><span>1</span><span>+</span><span>exp</span><span>(</span>beta0<span>*</span><span>rep</span><span>(</span><span>1</span>,<span>length</span><span>(</span>y<span>)</span><span>)</span> <span>+</span> X<span>%*%</span><span>beta</span><span>)</span><span>)</span>
  z <span>=</span> beta0<span>*</span><span>rep</span><span>(</span><span>1</span>,<span>length</span><span>(</span>y<span>)</span><span>)</span> <span>+</span> X<span>%*%</span><span>beta</span> <span>+</span> <span>(</span>y<span>-</span>p<span>)</span><span>/</span><span>(</span>p<span>*</span><span>(</span><span>1</span><span>-</span>p<span>)</span><span>)</span>
  omega <span>=</span> p<span>*</span><span>(</span><span>1</span><span>-</span>p<span>)</span><span>/</span><span>(</span><span>sum</span><span>(</span><span>(</span>p<span>*</span><span>(</span><span>1</span><span>-</span>p<span>)</span><span>)</span><span>)</span><span>)</span>
    beta0list <span>=</span> <span>numeric</span><span>(</span><span>length</span><span>(</span>maxiter<span>+</span><span>1</span><span>)</span><span>)</span>
    beta0 <span>=</span> <span>sum</span><span>(</span>y<span>-</span>X<span>%*%</span><span>beta</span><span>)</span><span>/</span><span>(</span><span>length</span><span>(</span>y<span>)</span><span>)</span>
    beta0list<span>[</span><span>1</span><span>]</span> <span>=</span> beta0
    <span>for</span> <span>(</span>j <span>in</span> <span>1</span><span>:</span>maxiter<span>)</span><span>{</span>
      <span>for</span> <span>(</span>k <span>in</span> <span>1</span><span>:</span><span>length</span><span>(</span><span>beta</span><span>)</span><span>)</span><span>{</span>
        r <span>=</span> z <span>-</span> X<span>[</span>,<span>-</span>k<span>]</span><span>%*%</span><span>beta</span><span>[</span><span>-</span>k<span>]</span> <span>-</span> beta0<span>*</span><span>rep</span><span>(</span><span>1</span>,<span>length</span><span>(</span>y<span>)</span><span>)</span>
       <span>beta</span><span>[</span>k<span>]</span> <span>=</span> <span>(</span><span>1</span><span>/</span><span>sum</span><span>(</span>omega<span>*</span>X<span>[</span>,k<span>]</span><span>^</span><span>2</span><span>)</span><span>)</span><span>*</span>soft_thresholding<span>(</span><span>t</span><span>(</span>omega<span>*</span>r<span>)</span><span>%*%</span>X<span>[</span>,k<span>]</span>,<span>length</span><span>(</span>y<span>)</span><span>*</span>lambda<span>)</span>
      <span>}</span>
      beta0 <span>=</span> <span>sum</span><span>(</span>y<span>-</span>X<span>%*%</span><span>beta</span><span>)</span><span>/</span><span>(</span><span>length</span><span>(</span>y<span>)</span><span>)</span>
      beta0list<span>[</span>j<span>+</span><span>1</span><span>]</span> <span>=</span> beta0
      betalist<span>[</span><span>[</span>j<span>+</span><span>1</span><span>]</span><span>]</span> <span>=</span> <span>beta</span>
      obj<span>[</span>j<span>]</span> <span>=</span> <span>(</span><span>1</span><span>/</span><span>2</span><span>)</span><span>*</span><span>(</span><span>1</span><span>/</span><span>length</span><span>(</span>y<span>)</span><span>)</span><span>*</span>norm<span>(</span>omega<span>*</span><span>(</span>z <span>-</span> X<span>%*%</span><span>beta</span> <span>-</span> 
beta0<span>*</span><span>rep</span><span>(</span><span>1</span>,<span>length</span><span>(</span>y<span>)</span><span>)</span><span>)</span>,<span>'F'</span><span>)</span><span>^</span><span>2</span> <span>+</span> lambda<span>*</span><span>sum</span><span>(</span><span>abs</span><span>(</span><span>beta</span><span>)</span><span>)</span>
  p <span>=</span> <span>exp</span><span>(</span>beta0<span>*</span><span>rep</span><span>(</span><span>1</span>,<span>length</span><span>(</span>y<span>)</span><span>)</span> <span>+</span> X<span>%*%</span><span>beta</span><span>)</span><span>/</span><span>(</span><span>1</span><span>+</span><span>exp</span><span>(</span>beta0<span>*</span><span>rep</span><span>(</span><span>1</span>,<span>length</span><span>(</span>y<span>)</span><span>)</span> <span>+</span> X<span>%*%</span><span>beta</span><span>)</span><span>)</span>
  z <span>=</span> beta0<span>*</span><span>rep</span><span>(</span><span>1</span>,<span>length</span><span>(</span>y<span>)</span><span>)</span> <span>+</span> X<span>%*%</span><span>beta</span> <span>+</span> <span>(</span>y<span>-</span>p<span>)</span><span>/</span><span>(</span>p<span>*</span><span>(</span><span>1</span><span>-</span>p<span>)</span><span>)</span>
  omega <span>=</span> p<span>*</span><span>(</span><span>1</span><span>-</span>p<span>)</span><span>/</span><span>(</span><span>sum</span><span>(</span><span>(</span>p<span>*</span><span>(</span><span>1</span><span>-</span>p<span>)</span><span>)</span><span>)</span><span>)</span>
      <span>if</span> <span>(</span>norm<span>(</span><span>rbind</span><span>(</span>beta0list<span>[</span>j<span>]</span>,betalist<span>[</span><span>[</span>j<span>]</span><span>]</span><span>)</span> <span>-</span> 
<span>rbind</span><span>(</span>beta0,<span>beta</span><span>)</span>,<span>'F'</span><span>)</span> <span>&amp;</span>lt<span>;</span> tol<span>)</span> <span>{</span> <span>break</span> <span>}</span> 
        <span>}</span> 
<span>return</span><span>(</span><span>list</span><span>(</span>obj<span>=</span>obj<span>[</span><span>1</span><span>:</span>j<span>]</span>,<span>beta</span><span>=</span><span>beta</span>,intercept<span>=</span>beta0<span>)</span><span>)</span> <span>}</span></pre></td></tr></table></div>
<p><img alt="" src="https://i2.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-04-a%CC%80-15.32.01.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>It looks like what can get when calling glmnet… and here, we do have null components for some <img alt="\lambda" title="\lambda" src="https://s0.wp.com/latex.php?latex=%5Clambda&bg=ffffff&fg=000&s=0" class="latex"> large enough ! Really null… and that’s cool actually.</p>
<h2>Application on our second dataset</h2>
<p>Consider now the second dataset, with two covariates. The code to get lasso estimates is</p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
</pre></td><td><pre>df0 <span>=</span> <span>df</span>
df0$y <span>=</span> <span>as.<span>numeric</span></span><span>(</span><span>df</span>$y<span>)</span><span>-</span><span>1</span>
plot_lambda <span>=</span> <span>function</span><span>(</span>lambda<span>)</span><span>{</span>
m <span>=</span> <span>apply</span><span>(</span>df0,<span>2</span>,<span>mean</span><span>)</span>
s <span>=</span> <span>apply</span><span>(</span>df0,<span>2</span>,<span>sd</span><span>)</span>
<span>for</span><span>(</span>j <span>in</span> <span>1</span><span>:</span><span>2</span><span>)</span> df0<span>[</span>,j<span>]</span> <span>&amp;</span>lt<span>;-</span> <span>(</span>df0<span>[</span>,j<span>]</span><span>-</span>m<span>[</span>j<span>]</span><span>)</span><span>/</span>s<span>[</span>j<span>]</span>
reg <span>=</span> glmnet<span>(</span><span>cbind</span><span>(</span>df0$x1,df0$x2<span>)</span>, df0$y<span>==</span><span>1</span>, alpha<span>=</span><span>1</span>,lambda<span>=</span>lambda<span>)</span>
u <span>=</span> <span>seq</span><span>(</span><span>0</span>,<span>1</span>,<span>length</span><span>=</span><span>101</span><span>)</span>
p <span>=</span> <span>function</span><span>(</span>x,y<span>)</span><span>{</span>
  xt <span>=</span> <span>(</span>x<span>-</span>m<span>[</span><span>1</span><span>]</span><span>)</span><span>/</span>s<span>[</span><span>1</span><span>]</span>
  yt <span>=</span> <span>(</span>y<span>-</span>m<span>[</span><span>2</span><span>]</span><span>)</span><span>/</span>s<span>[</span><span>2</span><span>]</span>
  <span>predict</span><span>(</span>reg,newx<span>=</span><span>cbind</span><span>(</span>x1<span>=</span>xt,x2<span>=</span>yt<span>)</span>,type<span>=</span><span>&quot;response&quot;</span><span>)</span><span>}</span>
v <span>=</span> <span>outer</span><span>(</span>u,u,p<span>)</span>
<span>image</span><span>(</span>u,u,v,<span>col</span><span>=</span>clr10,breaks<span>=</span><span>(</span><span>0</span><span>:</span><span>10</span><span>)</span><span>/</span><span>10</span><span>)</span>
<span>points</span><span>(</span><span>df</span>$x1,<span>df</span>$x2,pch<span>=</span><span>19</span>,cex<span>=</span><span>1.5</span>,<span>col</span><span>=</span><span>&quot;white&quot;</span><span>)</span>
<span>points</span><span>(</span><span>df</span>$x1,<span>df</span>$x2,pch<span>=</span><span>c</span><span>(</span><span>1</span>,<span>19</span><span>)</span><span>[</span><span>1</span><span>+</span>z<span>]</span>,cex<span>=</span><span>1.5</span><span>)</span>
<span>contour</span><span>(</span>u,u,v,<span>levels</span> <span>=</span> .5,add<span>=</span>TRUE<span>)</span><span>}</span></pre></td></tr></table></div>
<p>Consider some small values, for [\lambda], so that we only have some sort of shrinkage of parameters,</p>
<div><table><tr><td><pre>1
2
3
4
5
</pre></td><td><pre>reg <span>=</span> glmnet<span>(</span><span>cbind</span><span>(</span>df0$x1,df0$x2<span>)</span>, df0$y<span>==</span><span>1</span>, alpha<span>=</span><span>1</span><span>)</span>
<span>par</span><span>(</span>mfrow<span>=</span><span>c</span><span>(</span><span>1</span>,<span>2</span><span>)</span><span>)</span>
<span>plot</span><span>(</span>reg,xvar<span>=</span><span>&quot;lambda&quot;</span>,<span>col</span><span>=</span><span>c</span><span>(</span><span>&quot;blue&quot;</span>,<span>&quot;red&quot;</span><span>)</span>,lwd<span>=</span><span>2</span><span>)</span>
<span>abline</span><span>(</span>v<span>=</span><span>exp</span><span>(</span><span>-</span><span>2.8</span><span>)</span><span>)</span>
plot_lambda<span>(</span><span>exp</span><span>(</span><span>-</span><span>2.8</span><span>)</span><span>)</span></pre></td></tr></table></div>
<p><img alt="" src="https://i2.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-02-a%CC%80-22.55.12.png?w=456&ssl=1" data-recalc-dims="1"><br>
But with a larger <img alt="\lambda" title="\lambda" src="https://s0.wp.com/latex.php?latex=%5Clambda&bg=ffffff&fg=000&s=0" class="latex">, there is variable selection: here <img alt="\widehat{\beta}_{1,\lambda}=0" title="\widehat{\beta}_{1,\lambda}=0" src="https://s0.wp.com/latex.php?latex=%5Cwidehat%7B%5Cbeta%7D_%7B1%2C%5Clambda%7D%3D0&bg=ffffff&fg=000&s=0" class="latex"></p>
<div><table><tr><td><pre>1
2
3
4
5
</pre></td><td><pre>reg <span>=</span> glmnet<span>(</span><span>cbind</span><span>(</span>df0$x1,df0$x2<span>)</span>, df0$y<span>==</span><span>1</span>, alpha<span>=</span><span>1</span><span>)</span>
<span>par</span><span>(</span>mfrow<span>=</span><span>c</span><span>(</span><span>1</span>,<span>2</span><span>)</span><span>)</span>
<span>plot</span><span>(</span>reg,xvar<span>=</span><span>&quot;lambda&quot;</span>,<span>col</span><span>=</span><span>c</span><span>(</span><span>&quot;blue&quot;</span>,<span>&quot;red&quot;</span><span>)</span>,lwd<span>=</span><span>2</span><span>)</span>
<span>abline</span><span>(</span>v<span>=</span><span>exp</span><span>(</span><span>-</span><span>2.1</span><span>)</span><span>)</span>
plot_lambda<span>(</span><span>exp</span><span>(</span><span>-</span><span>2.1</span><span>)</span><span>)</span></pre></td></tr></table></div>
<p><img alt="" src="https://i2.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-02-a%CC%80-22.54.26.png?w=456&ssl=1" data-recalc-dims="1"></p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://freakonometrics.hypotheses.org/52894"> R-english – Freakonometrics</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/7Lj4oRVrGME" height="1"></a></div>