After equation 2 in the paper, the author writes:

> where $\mathbf{w}^{(i)}$ denotes the $i^{th}$ Monte Carlo sample drawn from the **variational posterior** $q(\mathbf{w}^{(i)}|\theta)$.

In section 3.2 titled "Gaussian variational posterior":

> Suppose that the **variational posterior is a diagonal Gaussian distribution,  then a sample of the weights $\mathbf{w}$ can be obtained by sampling a unit Gaussian, shifting it by a mean $\mu$ and scaling by a standard deviation $\sigma$**. We parameterise the standard deviation pointwise as $\sigma = \log(1 + \exp(\rho))$ and so $\sigma$ is always non-negative. The variational posterior parameters are $\theta = (\mu,\rho)$. Thus the transform from a sample of parameter-free noise and the variational posterior parameters that yields a posterior sample of the weights $\mathbf{w}$ is: $\mathbf{w} = t(\theta,\epsilon) = \mu + \log(1 + \exp(\rho)) \circ \epsilon$ where $\circ$ is point-wise multiplication.

However, in your question you wrote
\begin{align}
\text{prior} &= \log(q(\mathbf{w}|\mu,\rho)) = \sum_i \log(p(w_i | 0, 1)) \\
\text{posterior} &= \log(P(\mathbf{w})) = \sum_i \log(p(w_i | \mu, \sigma^2)) \\
\text{likelihood} &= \log(P(\mathcal{D}|\mathbf{w})) = y \cdot \log(\text{softmax}(\hat{y}))
\end{align}
which does not seem to be correct, since $q(\mathbf{w}|\mu,\rho)$ is the posterior according to the author. Here is what I think the authors meant
\begin{align}
\text{posterior} &= q(\mathbf{w}|\theta) \\
\text{prior} &= P(\mathbf{w}) \\
\text{likelihood} &= P(\mathcal{D}|\mathbf{w})
\end{align}
Also, the authors did not specify the prior $P(\mathbf{w})$ until equation 7 in section 3.3, which is a mixture of Gaussians and not a standard Gaussian as you wrote. Just wanted to point this out.

Here is a rough outline of what you can do in torch to implement this:
1. Create the dataset $\mathcal{D}$ consisting of the vectors $\mathbf{x}_1,\mathbf{x}_2,\dots,\mathbf{x}_N$ and the labels $y_1,y_2,...,y_N$. This can be done by sampling $\mathbf{x}$ vectors from some multivariate distribution, and sampling $y$ labels from a Bernoulli distribution.
2. Sample $\epsilon$ from the multivariate normal distribution $N(0,I)$.
3. Sample initial values for $\mu$ and $\rho$ <s>from the multivariate normal distribution $N(0,I)$ </s>. The authors do not say anything about whether $\mu$ and $\rho$ are themselves random variables with an associated distribution $p(\theta)$ or are treated as unknown constants. For simplicity, I will assume they are treated as unknown constants. In that case, their initial values can be sampled from any distribution. For simplicity, you can sample them using the `torch_randn` function. You only need to do this once to start performing gradient descent. Note that $\mu$ is a column vector and $\rho$ is <s>a diagonal matrix</s> also a column vector because the authors use the element-wise multiplication operation $\circ$ in their paper to multiply $\epsilon$ by $\log(1 + \exp(\rho))$.
4. Compute $\mathbf{w} = \mu + \log(1 + \exp(\rho)) \circ \epsilon$.
5. Compute $\log(q(\mathbf{w}|\mu,\rho))$ by inputting the $\mathbf{w}$,$\mu$, and $\rho$ that you obtained in steps 3 and 4 into the multivariate normal probability density function.
6. Compute $\log(P(\mathbf{w}))$ by inputting the $\mathbf{w}$ that you obtained in step 3 into the multivariate normal probability density function with mean $\mathbf{0}$ and covariance $I$. **Note** again that this does not follow what the authors did in the paper, and that you would instead need to implement equation 7 in the paper.
7. Compute $\log(P(\mathcal{D}|\mathbf{w}))$. Note that you did not specify what $\hat{y}$ is. For the sake of simplicity, I will assume that $\hat{y} = \mathbf{w}^T \mathbf{x}$. This means that $\log(P(\mathcal{D}|\mathbf{w})) = -\sum_{i=1}^N y_i \cdot \log(\hat{y}_i)$. This is just the [cross-entropy](https://en.wikipedia.org/wiki/Cross_entropy) between the $y$ labels and the $\hat{y}$ labels. Note that each $y$ and $\hat{y}$ **must** range between 0 and 1.
8. Compute $f(\mathbf{w},\theta) = \log(q(\mathbf{w}|\theta)) - \log(P(\mathbf{w})) - \log(P(\mathcal{D}|\mathbf{w}))$.
9. Compute the gradients as given in steps 5 and 6 of the algorithm.
10. Update $\mu$ and $\rho$ as given in step 7 of the algorithm.
11. Repeat the following aforementioned steps until convergence: step 2 (sampling $\epsilon$) $\rightarrow$ step 4 (computing $\mathbf{w}$) $\rightarrow$ steps 5,6,7,8,9, and 10.

I am aware that I did not mention how to compute the gradients in torch. If that is still something that you are not sure about, let me know.

## Response to comments

> What is still unclear to me is how to design the hidden layers of neurons. Am i correct in understanding that if i have 100 observations, for (2) this means i sample 100 $\epsilon$'s from a mv normal distribution with $\mu=0$ and a 100x100 identity matrix? for (3) i sample 100 $\mu$'s as a column vector and $\rho$ is 100x100 diagonal matrix of 100 samples of the same mv normal distribution? and for (4) i will get 100 samples $\mathbf{w}$? For my hidden layer i have specified 32 neurons. Where do these come into play?

Suppose you have $N$ observations in your dataset $\mathcal{D}$ of $(\mathbf{x},y)$ pairs, as discussed in step 1 above, and suppose that each $\mathbf{x}$ is $K$-dimensional such that
$$
\mathbf{x} =
\begin{bmatrix}
x_1 \\
x_2 \\
\vdots \\
x_K
\end{bmatrix}
$$
Also, for the sake of simplicity, suppose that your network consists of a single layer and the output of your network is
$$
\hat{y} = \phi(\mathbf{W}^T \mathbf{x} + \mathbf{b})
$$
where $\mathbf{W}$ is a $K \times M$ matrix of weights, where $M$ is the number of possible classes, $\mathbf{b}$ is a $M \times 1$ column vector, and $\phi(\cdot)$ is the softmax function. Recall that the authors wrote
$$
\mathbf{w} = \mu + \log(1 + \exp(\rho)) \circ \epsilon
$$
This be generalized to a $K \times M$ matrix $\mathbf{W}$ by sampling a $K \times M$ matrix $\epsilon$ using `torch_randn`, sampling **initial values** (only need to sample once) of $K \times M$ matrices $\mu$ and $\rho$ using `torch_randn`. Since $\epsilon \sim N(0,I)$, then you can also sample a $K \times M$ matrix $\epsilon$ using `torch_randn`. You can then obtain $\mathbf{W}$ as shown above. Each row of $\mathbf{W}$ represents a single weight vector. You can repeat this process to obtain $\mathbf{b}$:
$$
\mathbf{b} = \mu + \log(1 + \exp(\rho)) \circ \epsilon
$$
where these are now column vectors instead of matrices. You can then compute $\hat{y}$ as shown above.

> For step 7, $\hat{y}$ is supposed to represent the predicted values, which are compared to the ground truth $y$. So i think we're talking about the same thing here. You mention both $y$ and $\hat{y}$ have to range between 0 and 1, does that mean i should apply a sigmoid activation function after the output layer? 

I think we have slightly different definitions of $\hat{y}$. I am defining $\hat{y}$ as
$$
\hat{y} = \phi(\mathbf{W}^T \mathbf{x} + \mathbf{b})
$$