Constructing multilevel regression design matrix 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?
 A: Not knowing about the type of data you are using, but having been stuck in similar situations in the past: 
1) If your groups are sub-members (as would be mesh blocks in a census), then it may be cheap to define a hierarchy on larger groups. To use the analogy of a census, this will depend on whether you think there are sufficiently informative differences in the mesh blocks within the same local government area, or whether the differences between local government areas would be sufficient. 
2) Your problem may not be all that bad (apart from the computation), so long as each group has enough observations. Chapter 18 of Gelman and Hill's book goes through this quite intuitively. Essentially, if each of your groups has only a few (or zero) observations, then you just get pooled estimate; the more observations the group has, the further from the pooled estimate its parameter estimate could possibly be. 
3) Recall you need a weighting vector as well, to account for the possible differences in the number of observations in each group.
4) I would consider cutting down the number of varying slopes, as this adds a lot of columns. Though, again, this depends on your research question. 
Good luck!
