Is there a way to write Elastic Net in expanded matrix form? I am working through a regression problem for a matrix of data that isn't full rank and has more features than observations. For these reasons, I'd like to use elastic net because of its $L1$ and $L2$ regularization parameters.
While I've already found tutorials on how to use elastic net in python in practice, I'm still trying to fully understand how the regularization parameters fit into the estimation of $\hat{\beta}_i$.
I know that in an OLS model the expanded matrix form of regression looks like the following:
\begin{equation} 
{\begin{bmatrix}
 y_1  \\
 \vdots  \\
 y_n 
 \end{bmatrix}} = {\begin{bmatrix}
                   x_{11} & \ldots & x_{1k} \\
                   \vdots & \ddots & \vdots \\
                   x_{n1} & \ldots & x_{nk}
                   \end{bmatrix}} \,\cdot\, {\begin{bmatrix}
                                             \beta_1  \\
                                             \vdots  \\
                                             \beta_k 
                                             \end{bmatrix}} \,+\, {\begin{bmatrix}
                                                                   \epsilon_1  \\
                                                                   \vdots  \\
                                                                   \epsilon_n 
                                                                   \end{bmatrix}}
\end{equation}
This notation is tremendously helfpul in trying to understand how we derive $\hat{\beta_i}$ under OLS. Is it possible to write elastic net in a similarly expanded form, or does that not work because of the $L1$ penalty on the loss? If it is possible, where and how would I add the $L1$ and $L2$ penalties to the OLS equation above?

In short: To the best of my knowledge, no.
In detail:
The general elastic net solution is
$$\min_{\beta}{\{ \| y-X\beta\|_2^2+\gamma \|\beta\|_1+\lambda \|\beta\|_2^2 \}}$$
Ridge solution is defined as $\hat{\beta}^R=(X^TX+\lambda I_p)^{-1}X^Ty$ (contrary to the simple OLS solution $\hat{\beta}=(X^TX)^{-1}X^Ty$), so we can work something out:
$$\hat{\beta}^R=(X^TX+\lambda I_p)^{-1}X^Ty\\
(X^TX+\lambda I_p)\hat{\beta}^R=X^Ty\\
(X^T)^{-1}(X^TX+\lambda I_p)\hat{\beta}^R=y\\
y=(X+\lambda (X^T)^{-1})\hat{\beta}^R$$
That is, $y=(X+\lambda (X^T)^{-1})\beta+\epsilon$ assuming that $X^T$ is invertible. Note that $\hat{\beta}^R$ is not an unbiased estimator of the OLS problem $y=X\beta+\epsilon$.
While the $L_2$ penalty can somehow be expressed in matrix form, the closed form LASSO solution is given by $$\beta_j^{LASSO}=sign(\beta_j^{LS})(|\beta_j^{LS}|-\gamma)^+$$ where $\beta_j^{LS}$ is the jth coordinate of the OLS solution. An alternative representation is $\beta$ which minimizes $(Y-X\beta)^T(Y-X\beta)+\gamma\sum_j{|\beta_j|}$. As far as I know, there is no representation of this problem in matrix form, and thus we cannot represent the general elastic net in expanded matrix form.
