1) Yes
2) No (or Maybe, depending on what you mean).
Hierarchical Bayesian Regression Models are parametric, generative models. Regression is the technique of determining the parameters. The fact that the models are Bayesian implies that (in the model) the parameters of the lowest level data generating distributions are themselves drawn from higher-level (prior) distributions. Some of the parameters of these distributions are typically referred to as hyper-parameters.
In the regression procedure these hyper-parameters are held fixed, and cross validation is used to determine optimal values for these hyper-parameters. The most common form of estimation used in these contexts is Maximum a Posteriori (MAP) regression. In these instances, many commonly used priors lead to relatively simple terms in the objective function. Adding these terms to the optimization function is called regularization. In the following examples $w$ is the vector of parameters in the model.
- Normal Priors with fixed Variance -> L2 regularization (add $\lambda * \|w \|_2$) to the objective function.
- Laplacian Priors with fixed Variance -> L1 regularization (add $\lambda * \|w \|_1$) to the objective function.
Cross validation is used to determine $\lambda$. Now, back to the "maybe." L1 regularization is well known for the fact that it tends to lead to sparse parameter vectors ($w_i = 0$ for many $i$). In theory, one would love to try regressing with every possible subset of features, testing the results on a validation set. Unfortunately, this is computationally very challenging (impossible in contexts with thousands of variables).
I hope this answers your question and gives you the terminology necessary to find any other information needed.
EDIT
To be clear, there are two issues here. Setting up the (parametric) model, and estimating the parameters. Consider this very simple model
$$
w_i \sim \mathcal{N}(0, \sigma^2), i \in \{1, \dots, n \},
$$
$$
y_{i,j} \sim \mathcal{N}(w_i, \tau^2), j \in \{1, \dots, n_i \}.
$$
Our data are the "y"s and we know their corresponding "i"s. So we've drawn many samples from normal distributions, whose means were themselves drawn from one normal distribution. I will refer to our data as $Y$ and to some vector of possible "w"s as $w$. Our goal is to find probable values for $w$, given $Y$.
$$
P(w | Y) = \frac{P(Y|w) \cdot P(w)}{P(Y)}.
$$
The numerator is easy to compute from our model, as we will see. The denominator is more difficult to compute (in general) because we would have to marginalize $w$ out of the joint distribution (for which we have an equation). But since the denominator is always positive we can maximize the whole thing by maximizing the numerator, or the log of the numerator (since log is monotone increasing). We will actually minimize the negative log (convention). So,
\begin{eqnarray}
P(w | Y) &\propto& P(Y|w) \cdot P(w) \\
&\propto& \left[\prod_i\prod_j exp\left(-\frac{(y_{i,j} - w_i)^2}{2\tau^2}\right) \right] \cdot \prod_i exp\left(-\frac{w_i^2}{2\sigma^2}\right).
\end{eqnarray}
Taking the negative log, multiplying by $2\tau^2$, and thinking of it as a function of $w$, we have
\begin{eqnarray}
f(w) &=& \left[ \Sigma_i \Sigma_j (y_{i,j} - w_i)^2 \right] + \frac{\tau^2}{\sigma^2} \cdot \left[ \Sigma_i w_i^2 \right] \\
&=& \left[ \Sigma_i \Sigma_j (y_{i,j} - w_i)^2 \right] + \lambda \cdot \left[ \Sigma_i w_i^2 \right].
\end{eqnarray}
This is what we minimize. The first term is the log-loss and the second term is the regularization. If we minimized only the first term out estimates would be $w_i = \Sigma_j y_{i,j} / n_{i}$. The effect of the regularization will be to pull all of these $w_i$ toward 0. It turned out that variances of our distributions only mattered relative to one another, so I replaced their quotient with $\lambda$. The larger $\lambda$ is, the more $w_i$ will be pulled toward 0. In order to determine $\lambda$ we can train our model on a subset of our data, by minimizing the below function, then evaluate our model on the remainder of our data by computing only the logloss term.
