Graphical path Coordinate Descent in case of semi-differentiable functions such as Lasso I am trying to understand how the graphical solution path to the optimum would look in the case of Lasso Regression. I can find only Pictures for the differentiable or non differentiable case. The picture below is the solution path for a convex and differentiable function.
How would the solution path look for the lasso case? Would there be a way to produce such a graph including the solution path in R?

The contours would look like this:

 A: Context and background
Recall that the Lasso minimization problem can be expressed as:
$$ \hat \theta_{lasso} = argmin_{\theta \in \mathbb{R}^n} \sum_{i=1}^m (y_i - \mathbf{x_i}^T \theta)^2 + \lambda \sum_{j=1}^n | \theta_j| \ $$
Which can be viewed as the minimization of two terms: $OLS + L_1$. 


*

*The first OLS term can be written as $(y - X \theta)^T(y - X \theta) $ which gives rise to an elipse contour plot centered around the Maximum Likelihood Estimator. 

*The second $L_1$ term is the equation of a diamond centered around 0 (or a romboid in higher dimensions) 

*The solution to the constrained optimization lies at the intersection between the contours of the two functions, and this intersection varies as a function of $\lambda$. For $\lambda = 0$ the solution is the MLE (as usual) and for $\lambda = \infty$ the solution is at $[0,0]$. 

*Since at the vertices of the diamond, one or many of the variables have value 0, there is a probability that one or many of the features will have a value exactly equal to 0. 

*In fact, I have read on this website that in higher dimensions, the probability of the intersection taking place at the vertices increases, but I will leave it to someone else to explain this intuition. 


Vsualizing the solution as a function of $\lambda$
Intuitively, as we increase $\lambda$ from 0 to $\infty$ we would expect the lasso solution to move from the OLS solution to the $L_1$ solution, which is $0$. You will find below the result of a simulation I have run for Ridge and Lasso regression on simulated data. 

Some comments on this simulation
The data is simulated from a sine curve with uniform noise added. The data is normalized, to facilitate the lasso coordinate descent, and centered, to remove the need for an intercept coefficient. I have added a polynomial term as the second feature:
$$ y = \theta_1 x + \theta_2 x^2$$
This means that due to the polynomial feature, the $x_1$ and $x_2$ features are highly correlated, in fact a quick check gives the following:


*

*Condition number of matrix $X$: $= 7.87$

*Covariance matrix $X.T @ X =  \begin{bmatrix}1. & 0.968 \\ 0.968 & 1. \end{bmatrix}$


This will translate into a "flattened" OLS solution in the shape of a valley.  
As you see, the coefficient path moves towards the $(0,0)$ point until hitting a vertex, at which point, one of the two coefficient is equal to zero. This can be seen by using a lasso path diagram

Visualizing the coordinate descent algorithm steps
Coordinate descent minimizes along coordinate directions. At each iteration, the algorithm determines a coordinate, then minimizes over the corresponding hyperplane while fixing all other coordinates. 
Here are two versions of coordinate descent on the same dataset as previously, the first for a low value of $\lambda$, the second for a high value of $\lambda$


Source
For anyone interested, you can find most of the code and associated mathematical derivations on my blog and at this page
