How prior distribution over neural networks parameters is implemented in practice? I have read some papers that mentioned about using standard normal distribution as prior distribution over deep neural networks parameters or use normal distribution with some configurable variance. But so far I can't find information on how this is actually implemented in practice. In my naive understanding, for every parameters in the DNN, I would draw samples from normal distribution in the form of tensors with dimensionality matching the dimensionality of the DNN parameters. Is this how it is commonly done in practice, if not how the prior over DNN parameters supposed to be implemented?
Because I can't find any example implementation that sample normal distribution to set DNN parameters so far. For example, this paper mentioned that "the prior on the parameters is set to $p(\theta) = N(0, \sigma^2I)$", where $\sigma^2$ is set to 1 by default and 100 in one of the experiment scenario. But I can't find anything that sample from normal distribution to set the DNN parameters anywhere in the official source code for that paper.
I also vaguely remember reading about prior in DNN implemented in the form of regularization term or weight decay, but couldn't find any explanation on how regularization or weight decay maps to implementing normal distribution prior over the parameters i.e how the regularization or weight decay formula supposed to be changed if I would increase the variance of the normal distribution prior.

A zero-mean, isotropic multivariate gaussian prior on the network weights $\theta$ reduces to a penalty on the $L^2$ norm on the parameter vector $\theta$.
Finding the MAP estimate of the posterior reduces to maximizing the probability $p(y|x) = p(y|x,\theta)p(\theta) 
$, which is equivalent to minimizing the negative logarithm of the same:
$$\begin{align}
p(y|x) &= p(y|x,\theta)p(\theta) \\
-\ln p(y|x) &= L(y|x,\theta) - \ln p(\theta) \\
&= L(y|x,\theta) - 
\ln \left[
(2\pi)^{ -\frac{k}{2} } \det(\Sigma)^{ -\frac{1}{2} } \exp \left(
-\frac{1}{2} (\theta - \mu)^T \Sigma^{-1} (\theta - \mu)
\right)
\right] \\
&=  L(y|x,\theta) + \frac{1}{2}\theta^T \left(\sigma^2 I \right)^{-1}\theta +C \\
&=  L(y|x,\theta) + \frac{\sigma^{-2}}{2}\theta^T \theta \\
&= L(y|x,\theta) + \frac{\lambda}{2} \| \theta \|_2^2
\end{align}$$
where $ L(y|x,\theta)=-\ln p(y|x,\theta)$ is your loss function (e.g. mean square error or categorical cross-entropy loss), the negative log likelihood given the model, the parameters $\theta$, and the data $(x,y)$.
Some notes about this derivation:

*

*The last line makes the substitution $\sigma^{-2}=\lambda$ and writes the penalty as a norm to make the connection to ridge regression more apparent.


*We can neglect the constant additive terms $C=-\frac{1}{2}\left(k\ln(2\pi)+\ln|\Sigma|\right)$ because they do not depend on $\theta$; including them will change the value of the extrema, but not its location.


*This is given as a generic statement about any loss $L$ which can be expressed as the negative log of the probability, so if you're working on a classification problem, a regression problem, or any problem formulated as a probability model, you can just substitute the appropriate expression for $L$.
Of course, if you're interested in Bayesian methods, you might not wish to be constrained solely to MAP estimates of the model. Radford Neal looks at some methods to utilize the posterior distribution of $\theta$ in his book Bayesian Learning for Neural Networks, including MCMC to estimate neural networks. Since publication, there are probably many more works which have taken these concepts even further.
One could optimize this augmented loss function directly. Alternatively, it could be implemented as weight decay during training; PyTorch does it this way, for instance. The reason you might want to implement weight decay as a component of the optimizer (as opposed to just using on autograd on the regularized loss) is that the gradient update looks like
$$\begin{align}
\theta_{i+1} &= \theta_i - \eta \frac{\partial}{\partial \theta} \left[L + \frac{\lambda}{2}\| \theta \|_2^2 \right]\\
 &= \theta_i - \eta \frac{\partial L}{\partial \theta} - \eta \lambda  \theta_i \\
&= (1 - \eta \lambda) \theta_i - \eta \frac{\partial L}{\partial \theta}
\end{align}$$
where $\theta_i$ is the parameter vector at the $i$th optimization step and $\eta$ is the learning rate. But when using adaptive optimizers (e.g. Adam), the effect of weight decay is slightly different; see "Decoupled Weight Decay Regularization" by
Ilya Loshchilov and Frank Hutter.
