A common formulation of multilevel/hierarchical regression models is $y = Xb + Zc + e$, where $X$ is an $n \times p$ matrix of $p$ individual level predictors, $Z$ is an $n \times q$ matrix of $q$ group level predictors, $y$ is an $n \times 1$ vector of observations at the individual level and $e$ is the error term. Suppose there are $J$ groups and the $q$ group level predictors are continuous. I would like to understand how to construct a full design matrix from this formulation, for example to apply a gradient descent optimization algorithm.
Is this formulation equivalent to a design matrix constructed as: $p$ columns for the individual level predictors, $J$ columns of indicator variables for the $J$ groups (varying intercepts), $pJ$ columns for the interaction between the group indicators and the individual level predictors (varying slopes), and $q$ columns for the group level predictors (i.e., columns with repeated values for within-group observations).
Are there alternative ways to construct the design matrix when the number of groups is large ($J>10000$) to reduce the number of parameters?