Increasing multicollinearity in multilevel/hierarchical modeling?

Since you say :

However, the relationship between $textbf{y}$ and $textbf{X}$ is not strictly linear over the entire domain, but could be better modeled as such within several subgroups $g$ (and coefficients are more meaningful if defined for each subgroup):

it sounds to me very much like a mixed effects model (of which multilevel and hierarchical models are special cases) with random slopes for $textbf{X}$ in subgroups $g$. This will have the general form:

$$y = textbf{X}beta+textbf{Z}u+e$$

where $beta$ is a vector of fixed effects, $X$ and $Z$ are model matrices for the fixed effects and random effects respectively and $u$ and $e$ are vectors of random effects such that $E(u) = E(e) = 0$

In R you could fit such a model with, for example:

y = func(y ~ X1 + X2 + (X1 + X2 | g ), ...)

where func will be the relevant function from whatever package you choose, eg lme4 or GLMMAdaptive. Note that some packages, eg nlme use different syntax. This will estimate fixed effects (slopes) and random slopes for X1 and X2 and random intercepts for each group. If you do not want random intercepts - ie. you wish to allow the slopes to vary by group, but all pass through the same point on the y axis, then you would use:

y = func(y ~ X1 + X2 + (X1 + X2 + 0 | g ), ...)