The solution we propose to overcome these challenges is twofold: (1) employ a diagonal covariance structure for both the parameters \({\theta }\) , and hidden units \({Z}^{(j)}\) , and (2) use the inherent conditional independence of hidden units between layers, that is, \({Z}^{(j-1)}\protect {\protect \mathrel {\unknown. {\perp }\hspace{1.111pt}{\perp }}{Z}^{(j+1)}|{z}^{(j)}, in order to perform recursive layer-wise inference. As depicted in Figure \ref {FIG:FNN_inference_long}b, the first step consists in inferring the posterior mean vector and diagonal covariance for the output layer following\begin{equation}\begin{array}{rcl}f({z}^{(\mathtt {O})}|{y})&=&\mathcal {N}({z}^{(\mathtt {O})};{\mu }_{{Z}^{(\mathtt {O})}|{y}},{\Sigma }_{{Z}^{(\mathtt {O})}|{y}})\\[6pt]{\mu }_{{Z}^{(\mathtt {O})}|{y}}&=&{\mu }_{{Z}^{(\mathtt {O})}}+\mathbf {\Sigma }_{{Y\!Z}^{(\mathtt {O})}}^{\intercal }\mathbf {\Sigma }_{{Y}}^{-1}({y}-{\mu }_{{Y}})\\[4pt]\mathbf {\Sigma }_{{Z}^{(\mathtt {O})}|{y}}&=&\mathbf {\Sigma }_{{Z}^{(\mathtt {O})}}-\mathbf {\Sigma }_{{Y\!Z}^{(\mathtt {O})}}^{\intercal }\mathbf {\Sigma }_{{Y}}^{-1}\mathbf {\Sigma }_{{Y\!Z}^{(\mathtt {O})}}.\end{array}\end{equation}Note that for classification problems, because of the hierarchical formulation described in §\ref {SS:FFNN}, even if there are \mathtt {Y} classes, only \mathtt {H}=\lceil \log _2(\mathtt {Y})\rceil  hidden units from the output layer are updated for each observation.}\) In order to perform layer-wise inference for hidden units and parameters, we define the short-hand notation \(\lbrace {\theta }^{\text{+}},{Z}^{\text{+}}\rbrace \equiv \lbrace {\theta }^{(j+1)},{Z}^{(j+1)}\rbrace \)  and \(\lbrace {\theta },{Z}\rbrace \equiv \lbrace {\theta }^{(j)},{Z}^{(j)}\rbrace \) . As depicted in Figure REF b, the inference for the \(j^\text{th}\)  layer's weights \({\theta }\)  and hidden units \({Z}\)  is done using the Rauch-Tung-Striebel recursive procedure [1]} where
\(\begin{array}{rcl}f({z}|{y})&=&\mathcal {N}({z};{\mu }_{{Z}|{y}},{\Sigma }_{{Z}|{y}})\\[4pt]{\mu }_{{Z}|{y}}&=&{\mu }_{{Z}}+\mathbf {J}_{{Z}}\left({\mu }_{{Z}^{\text{+}}|{y}}-{\mu }_{{Z}^{\text{+}}}\right)\\[4pt]{\Sigma }_{{Z}|{y}}&=&{\Sigma }_{{Z}}+\mathbf {J}_{{Z}}\left({\Sigma }_{{Z}^{\text{+}}|{y}}-{\Sigma }_{{Z}^{\text{+}}}\right)\mathbf {J}_{{Z}}^{\intercal }\\[4pt]\mathbf {J}_{{Z}}&=&{\Sigma }_{{Z}\!{Z}^{\text{+}}}{\Sigma }_{{Z}^{\text{+}}}^{-1},\end{array}\) 
\(\begin{array}{rcl}f({\theta }|{y})&=&\mathcal {N}({\theta };{\mu }_{{\theta }|{y}},{\Sigma }_{{\theta }|{y}})\\[4pt]{\mu }_{{\theta }|{y}}&=&{\mu }_{{\theta }}+\mathbf {J}_{{\theta }}\left({\mu }_{{Z}^{\text{+}}|{y}}-{\mu }_{{Z}^{\text{+}}}\right)\\[4pt]{\Sigma }_{{\theta }|{y}}&=&{\Sigma }_{{\theta }}+\mathbf {J}_{{\theta }}\left({\Sigma }_{{Z}^{\text{+}}|{y}}-{\Sigma }_{{Z}^{\text{+}}}\right)\mathbf {J}_{{\theta }}^{\intercal }\\[4pt]\mathbf {J}_{{\theta }}&=&{\Sigma }_{{\theta }{Z}^{\text{+}}}{\Sigma }_{{Z}^{\text{+}}}^{-1}.\end{array}\) 
