\widetildebW
\begin{bmatrix}
1 \\
\bz_n
\end{bmatrix}-\ba_n
\right\Vert^2
\underbrace{\lambda_z \norm{\widetildebz_n}^2}_{=\lambda_z \norm{\bz_n}^2+\lambda_z}
\sum_{i\neq n} \norm{\widetildebW\widetildebz_i-\ba_i}^2  + \lambda_z \sum_{i\neq n}\norm{\widetildebz_i}^2   +    \lambda_w \norm{\widetildebW}^2 \\
&= \left\Vert
\begin{bmatrix}
\widebarbw_0 & \widebarbW
\end{bmatrix}
\begin{bmatrix}
1 \\
\bz_n
\end{bmatrix}-\ba_n
\right\Vert^2
+ \lambda_z \norm{\bz_n}^2 + C_{z_n}
= \left\Vert
\widebarbW \bz_n -
\underbrace{(\ba_n-\widebarbw_0)}_{\widebarba_n}
\right\Vert^2
+ \lambda_z \norm{\bz_n}^2 + C_{z_n},
\end{aligned}
\end{equation}

where $\widebarbw_0$ represents the first column of $\widetildebW$, $\widebarbW$ denotes the last $K$ columns of $\widetildebW$, and $C_{z_n}$ is a constant with respect to $\bz_n$. Let $\widebarba_n = \ba_n-\widebarbw_0$, the update for $\bz_n$ is just similar to the one in Equation (<ref>), with the gradient given by

\nabla L(\bz_n) = 2\widebarbW^\top\widebarbW\bz_n - 2\widebarbW^\top\widebarba_n + 2\lambda_z\bz_n.

Therefore, the update for $\bz_n$ is given by determining the root of the  gradient above:

\text{update for $\widetildebz_n$ is }
\left\{
\begin{aligned}
\bz_n &= (\widebarbW^\top\widebarbW+ \lambda_z\bI)^{-1} \widebarbW^\top \widebarba_n, \gap \text{for $n\in \{1,2,\ldots, N\}$};\\
\widetildebz_n &= \begin{bmatrix}
\end{bmatrix}.
\end{aligned}
\right.

Similarly, following the loss with respect to each row of $\bW$ in Equation (<ref>), suppose $\widetildebw_m =
\begin{bmatrix}
\bw_m \\
\end{bmatrix}$ is the $m$-th row of $\widetildebW$ (or $m$-th column of $\widetildebW^\top$), we have

\begin{equation}
\begin{aligned}
&\gap L(\bw_m ) \\
&=\norm{\widetildebZ^\top\widetildebW-\bA^\top}^2 +\lambda_w \norm{\widetildebW^\top}^2 + \lambda_z \norm{\widetildebZ}^2\\\
\norm{\widetildebZ^\top\widetildebw_m-\bb_m}^2 +
\underbrace{\lambda_w \norm{\widetildebw_m}^2}_{=\lambda_w \norm{\bw_m}^2+\lambda_w}
\sum_{i\neq m} \norm{\widetildebZ^\top\widetildebw_i-\bb_i}^2  + \lambda_w \sum_{i\neq m}\norm{\widetildebw_i}^2   +    \lambda_z \norm{\widetildebZ}^2 \\
\left\Vert
\begin{bmatrix}
\widebarbZ^\top&
\widebarbz_0
\end{bmatrix}
\begin{bmatrix}
\bw_m \\
\end{bmatrix}
-\bb_m\right\Vert^2 +
\lambda_w \norm{\bw_m}^2
\left\Vert
\widebarbZ^\top\bw_m
-(\bb_m-\widebarbz_0) \right\Vert^2+
\lambda_w \norm{\bw_m}^2
\end{aligned}
\end{equation}

where $\widebarbz_0$ represents the last column of $\widetildebZ^\top$, $\widebarbZ^\top$ contains the remaining columns of $\widetildebZ^\top$,

$C_{w_m}$ is a constant with respect to $\bw_m$, and $\bW^\top=[\bw_1, \bw_2, \ldots, \bw_M], \bA^\top=[\bb_1,\bb_2, \ldots, \bb_M]$ are the column partitions of $\bW^\top, \bA^\top$, respectively. Let $\widebarbb_m = \bb_m-\widebarbz_0$, the update for $\bw_m$ is again just similar to  the one in Equation (<ref>), with the gradient given by

\nabla L(\bw_md) = 2\widebarbZ\cdot \widebarbZ^\top\bw_m - 2\widebarbZ\cdot \widebarbb_m + 2\lambda_w\bw_m.

Therefore, the update for $\bw_m$ is given by the root of the  gradient above:

\text{update for $\widetildebw_m$ is }
\left\{
\begin{aligned}
\bw_m&=(\widebarbZ\cdot \widebarbZ^\top+\lambda_w\bI)^{-1}\widebarbZ\cdot \widebarbb_m,  \gap \text{for $m\in \{1,2,\ldots, M\}$} ;\\
\widetildebw_m &= \begin{bmatrix}
\bw_m \\ 1
\end{bmatrix}.
\end{aligned}
\right.

Similar updates through gradient descent with the consideration of the bias terms or treatment on missing entries can be deduced, and we shall not repeat the details (see Section <ref> and <ref> for a reference).

A gray flag image to be compressed. The size of the image is $600\times 1200$ with a rank of 402.

[All 5 singular values: $\sum_{i=1}^{5}\sigma_i\bu_i\bv_i^\top$,$F=$44379]

[Pseudoskeleton Rank 5 $\sum_{i=1}^{5}\bc_i\br_i^\top$,$G=45905$.]

[ALS Rank 5 $\sum_{i=1}^{5}\bw_i\bz_i^\top$,$S=52157$.]

Image compression for gray flag image into a rank-5 matrix via the SVD, and decompose into 5 parts where $\sigma_1 \geq \sigma_2 \geq \ldots \geq \sigma_{5}$, i.e., $F_1\leq F_2\leq \ldots \leq F_5$ with $F_i = \norm{\sigma_i\bu_i\bv^\top - \bA}_F$ for $i\in \{1,2,\ldots, 5\}$. And reconstruct images by single singular value and its corresponding left and right singular vectors, $\bc_i\br_i^\top$, $\bw_i\bz_i^\top$ respectively. Upper: SVD; Middle: Pseudoskeleton; Lower: ALS.

§ APPLICATIONS

§.§ Low-Rank Approximation

Figure <ref> shows an example of a gray image to be compressed. The size of the image is $600\times 1200$ with a rank of 402.
Low-rank approximation via the ALS decomposition, singular value decomposition (SVD), and Pseudoskeleton decomposition [Lu, 2021] can be applied to find the compression.

In Figure <ref> to <ref>, we approximate the image into a rank-5 matrix by truncated SVD: $\bA\approx \sum_{i=1}^{5}\sigma_i\bu_i\bv_i^\top$ [Lu, 2021]. It is known that the singular values contain the spectrum information with higher singular values containing lower-frequency information. And low-frequency contains more useful information [Leondes, 1995]. We find that the image, $\sigma_1\bu_1\bv_1^\top$, reconstructed by the first singular value $\sigma_1$, first left singular vector $\bu_1$, and first right singular vector $\bv_1$ is very close to the original flag image; and second to the fifth images reconstructed by the corresponding singular values and singular vectors contains more details of the flag to reconstruct the raw image.

Similar results can be observed for the low-rank approximation via the pseudoskeleton decomposition (See Lu, 2021 for more details). At a high level, the pseudoskeleton decomposition finds  the low-rank approximation by $\bA\approx \bC\bR$, where $\bC\in \real^{m\times \gamma}$ and $\bR\in \real^{\gamma\times n}$ if $\bA\in \real^{m\times n}$ such that $\bC$ and $\bR$ are rank-$\gamma$ matrices. Suppose $\gamma=5$, and

\bC=[\bc_1, \bc_2, \ldots, \bc_5],
\qquad
\text{and}
\qquad
\bR =
\begin{bmatrix}
\br_1^\top \\
\br_2^\top \\
\vdots \\
\br_5^\top
\end{bmatrix},

are the column and row partitions of $\bC$ and $\bR$, respectively. Then $\bA$ can be approximated by $\sum_{i=1}^{5}\bc_i\br_i^\top$. The partitions are ordered such that

\underbrace{\norm{\bc_1\br_1^\top-\bA}_F}_{G_1} \leq
\underbrace{\norm{\bc_2\br_2^\top-\bA}_F}_{G_2}
\leq \ldots \leq
\underbrace{\norm{\bc_5\br_5^\top-\bA}_F}_{G_5}.

We observe (in Figure <ref> to <ref>) that $\bc_1\br_1^\top$ works similarly to that of $\sigma_1\bu_1\bv^\top$, where the reconstruction errors measured by the Frobenius norm are  very close (60,464 in the pseudoskeleton case compared to that of 60,217 in the SVD case). This is partly because the pseudoskeleton decomposition relies on the SVD such that $\bc_1\br_1^\top$ internally has the largest “singular value" meaning in this sense [Lu, 2021].

Similarly, the ALS approximation is given by $\bA\approx\bW\bZ$, where $\bW\in \real^{m\times \gamma}$ and $\bZ\in \real^{\gamma\times n}$ if $\bA\in \real^{m\times n}$ such that $\bW$ and $\bZ$ are rank-$\gamma$ matrices. Suppose

\bW=[\bw_1, \bw_2, \ldots, \bw_5]
\qquad
\text{and}
\qquad
\bZ =
\begin{bmatrix}
\bz_1^\top \\
\bz_2^\top \\
\vdots \\
\bz_5^\top
\end{bmatrix},

are the column and row partitions of $\bW$ and $\bZ$, respectively [For simplicity, note that this definition is different from what we have defined in Section <ref>, where we define $\bw_i$'s as the rows of $\bW$. ]. Then $\bA$ can be approximated by $\sum_{i=1}^{5}\bw_i\bz_i^\top$. The partitions are ordered such that

\underbrace{\norm{\bw_1\bz_1^\top-\bA}_F}_{S_1} \leq
\underbrace{\norm{\bw_2\bz_2^\top-\bA}_F}_{S_2}
\leq \ldots \leq
\underbrace{\norm{\bw_5\bz_5^\top-\bA}_F}_{S_5}.

We observe (in Figure <ref> to <ref>) that $\bw_1\bz_1^\top$ works slightly differently to that of $\sigma_1\bu_1\bv^\top$, where the reconstruction errors measured by Frobenius norm are not close as well (82,727 in the ALS case compared to that of 60,217 in the SVD case). As we mentioned previously, $\bc_1\br_1^\top$ works similarly to that of $\sigma_1\bu_1\bv^\top$ since the pseudoskeleton relies on the SVD. However, in ALS, the reconstruction is all from least squares optimization.
The key difference between ALS and SVD is in the fact  that, in SVD, the importance of each vector in the basis is relative to the value of the singular value associated with that vector. This usually means that the first vector of the basis dominates and is the most used vector to reconstruct data; then the second vector and so on. So the bases in SVD have an implicit hierarchy and that doesn't happen in ALS, where we find the second component $\bw_2\bz_2^\top$ via ALS in Figure <ref> plays an important role in the reconstruction of the original figure; whereas the second component $\sigma_2\bu_2\bv_2^\top$ via SVD in Figure <ref> plays a small role in the reconstruction.

Comparison of reconstruction errors measured by Frobenius norm among the SVD, pseudoskeleton, and ALS, where the approximated rank ranges from $3$ to 100. ALS with well-selected parameters works similarly to SVD.

We finally compare low-rank approximation among the SVD, pseudoskeleton, and ALS with different ranks (3 to 100). Figure <ref> shows the difference of each compression with rank 90, 60, 30, 10. We observe that the SVD reconstructs well with rank 90, 60, 30. The pseudoskeleton approximation compresses well in the black horizontal and vertical lines in the image. But it performs poorly in the details of the flag. ALS works similarly to the SVD in terms of visual expression and reconstruction errors measured by the Frobenius norm.
Figure <ref> shows the comparison of the reconstruction errors among the SVD, the pseudoskeleton, and the ALS approximations measured by the Frobenius norm ranging from rank $3$ to $100$, where we find in all cases, the truncated SVD performs best in terms of Frobenius norm.
Similar results can be observed when applied to the spectral norm. The ALS works better than the pseudoskeleton decomposition when $\lambda_w=\lambda_z=0.15$. An interesting cutoff happens when $\lambda_w=\lambda_z=\{0.03, 0.08, 0.15\}$. That is, when the value of rank increases, the ALS will be very close to the SVD in the sense of low-rank approximation.

Truncated SVD

[SVD with rank 90Frobenius norm=6,498]

[Pseudoskeleton with rank 90Frobenius norm=13,751]

[ALS with rank 90Frobenius norm=6,622]

[SVD with rank 60Frobenius norm=8,956]

[Pseudoskeleton with rank 60Frobenius norm=14,217]

[ALS with rank 60Frobenius norm=9,028]

[SVD with rank 30Frobenius norm=14,586]

[Pseudoskeleton with rank 30Frobenius norm=17,853]

[ALS with rank 30Frobenius norm=18,624]

[SVD with rank 10Frobenius norm=31,402]

[Pseudoskeleton with rank 10Frobenius norm=33,797]

[ALS with rank 10Frobenius norm=33,449]

Image compression for gray flag image with different ranks.

§.§ Movie Recommender

The ALS is extensively developed for the movie recommender system. To see this, we obtain
the “MovieLens 100K" data set from MovieLens [Harper and Konstan, 2015] [http://grouplens.org]. It consists of 100,000 ratings from 943 users on 1,682 movies. The rating values go from 0 to 5. The data was collected through the MovieLens website
during the seven-month period from September 19th,
1997 through April 22nd, 1998. This data has been cleaned up—users
who had less than 20 ratings or did not have complete demographic
information were removed from this data set such that simple demographic info for the users (age, gender, occupation, zip) can be obtained. However, we will only work on the trivial rating matrix [In the next chapters, the data set is further cleaned by removing movies with less than 3 users such that 1,473 movies are kept accordingly.].

The data set is split into training and validation data, around 95,015 and 4,985 ratings, respectively.
The error is measured by root mean squared error (RMSE). The RMSE is frequently used as a measure of the difference between values. For a set of values $\{x_1, x_2, \ldots, x_n\}$ and its predictions $\{\hat{x}_1, \hat{x}_2, \ldots, \hat{x}_n\}$, the RMSE can be described as

\text{RMSE}(\bx, \hat{\bx}) = \sqrt{\frac{1}{n} \sum_{i=1}^{n}(x_i-\hat{x}_i)}.

The minimal RMSE for validation is obtained when $K=185$ and $\lambda_w=\lambda_z=0.15$, and it is equal to $0.806$ as shown in Figure <ref>. Therefore, when the rating ranges from 0 to 5, the ALS at least can predict whether the user likes to watch the movie (e.g., ranges 4 to 5) or not (e.g., ranges 0 to 2).

Comparison of training and validation error for “MovieLens 100K" data set with different reduction dimensions and regularization parameters.

Recommender 1.
A recommender system can work simply by suggesting the movie $m$ when $a_{mn}\geq4$ if user $n$ has not rated the movie $m$.

Recommender 2.
Or in rare cases, it happens that the user $n$ has rated all the movies he likes (say rates $\geq4$). Then  a partial solution is to find out similar movies to the high-rated movies to recommend. Suppose user $n$ likes movie $m$ very much and he has rated the movie $m$ with $5$: $a_{mn}=5$. Under the ALS approximation $\bA=\bW\bZ$, where each row of $\bW$ represents the hidden features of each movie (see Section <ref>  on the vector inner product), the solution is given by finding the most similar movies to movie $m$ that user $n$ has not rated (or watched). In mathematical language,

\mathop{\arg \max}_{\bw_i} \gap \text{similarity}(\bw_i, \bw_m), \qquad \text{for all} \gap i \notin \bo_n,

where $\bw_i$'s are the rows of $\bW$ representing the hidden feature of movie $i$, and $\bo_n$ is a mask vector indicating the movies that user $n$ has rated.

The method above relies on the similarity function between two vectors. The cosine similarity is the most commonly used measure. It is defined as the cosine of the angle between the two vectors:

\cos(\bx, \by) = \frac{\bx^\top\by}{\norm{\bx}\cdot \norm{\by}},

where the value ranges from $-1$ to 1, with $-1$ representing perfectly dissimilar and 1 being perfectly similar. Based on this definition, it follows that the cosine similarity depends only on the angle between the two nonzero vectors, but not on their magnitudes since it can be regarded as the inner product between the normalized  versions of these vectors.
Another measure for calculating similarity is known as the Pearson similarity:

\text{Pearson}(\bx,\by) =\frac{\Cov(\bx,\by)}{\sigma_x \cdot \sigma_y}
= \frac{\sum_{i=1}^{n} (x_i - \bar{x} ) (y_i -\bar{y})}{ \sqrt{\sum_{i=1}^{n} (x_i-\bar{x})^2}\sqrt{ \sum_{i=1}^{n} (y_i-\bar{y})^2 }},

whose range varies between $-1$ and 1, where $-1$ is perfectly dissimilar, 1 is perfectly similar, and 0 indicates no linear relationship. The Pearson similarity is usually used to measure the linear correlation between two sets of data. It is calculated as the ratio between the covariance of two variables and the product of their standard deviations.

Both Pearson correlation and cosine similarity are used in various fields, including machine learning and data analysis. Pearson correlation is commonly used in regression analysis, while cosine similarity is commonly used in recommendation systems and information retrieval, as we will see the cosine similarity performs better in the precision-recall (PR) curve analysis.

[Cosine Bin Plot.]

[Pearson Bin Plot.]

[PR Curve.]

Distribution of the insample and outsample using cosine and Pearson similarity, and the Precision-Recall curves for them.

Building upon the previous example  using the MovieLens 100K data set, we choose $\lambda_w=\lambda_z=0.15$ for  regularization and a rank of $62$ to minimize  RMSE.
We want to look at the similarity between different movie hidden vectors, and the goal is to see whether the matrix factorization can help differentiate high-rated from low-rated movies, so that the system can recommend the  movies correlated to the existing high-rated movies for each user.
Define further the term “insample" as the similarity between the movies having rates $5$ for each user, and “outsample" as the similarity between the movies having rates $5$ and $1$ for each user. Figure <ref> and <ref> depict the bin plots of the distributions of insample and outsample under cosine and Pearson similarities, respectively.
In both scenarios, a clear distinction is observed between the distributions of the “insample" and “outsample" data, showing that the ALS decomposition can actually find the hidden features of different movies for each user.
Figure <ref> shows the precision-recall (PR) curve for these scenarios, where we find the cosine similarity works better such that it can find out more than $73\%$ of the potential high-rated movies with a precision of $90\%$. However, Pearson similarity can  identify only about $64\%$ of the high-rated movies with the same precision.
In practice, other measures can also be explored, such as the negative Euclidean distance, in which case the Euclidean distance can measure the “dissimilarity" between two vectors; and a negative one thus represents the similarity between them.

* Prove Lemma <ref>.

* Find the update of column-by-column fashion for Algorithm <ref>.

* Constrained (Regularized) least squares (CLS).
Given $\bA\in\real^{m\times n}, \bb\in\real^{n}, \bB\in\real^{p\times n}$, and $\lambda\in\real_{++}$, we consider the regularized least squares problem

\mathop{\min}_{\bx\in\real^n} \norm{\bA\bx-\bb}^2 + \lambda\norm{\bB\bx}^2.

Show that the constrained least squares (CLS) problem has a unique solution if and only if $\nspace(\bA)\cap \nspace(\bB) = \{\bzero\}$.

* Weighted least squares (WLS).
Going further from the assumptions in Lemma <ref>, we consider further that each data point $i\in\{1,2,\ldots, m\}$ (i.e., each row of $\bA$) has a weight $w_i$.
This means some  data points may carry greater significance than others and there are ways to produce approximate minimzers that reflect this.
Show that the value $\bx_{WLS} = (\bA^\top\bW^2\bA)^{-1}\bA^\top\bW^2\bb$ serves as the weighted least squares (WLS)  estimator of $\bx$, where $\bW=\diag([w_1, w_2, \ldots, w_m])\in\real^{m\times m}$. Hint: find the normal equation for this problem.

* Restricted least squares (RLS).
Going further from the assumptions in Lemma <ref>, we consider further the restriction $\bx=\bC\bgamma+\bc$, where $\bC\in\real^{n\times k}$ is a known matrix such that $\bA\bC$ has full rank, $\bc$ is a known vector, and $\bgamma$ is an unknown vector.
Show that the value $\bx_{RLS}=\bC(\bC^\top\bA^\top\bA\bC)^{-1}(\bC^\top\bA^\top)(\bb-\bA\bc) +\bc$ serves as the restricted least squares (RLS) estimator of $\bx$.

* Find the restricted weighted least squares estimator.

*  First-order optimality condition for local optima points.
Consider the Fermat's theorem: for a one-dimensional function $g(\cdot)$ defined and differentiable over an interval ($a, b$), if a point $x^\star\in(a,b)$ is a local maximum or minimum, then $g^\prime(x^\star)=0$.
Prove the first-order optimality conditions for multivariate functions based on this Fermat's theorem for one-dimensional functions.
That is, consider function $f: \sS\rightarrow \real$ as a function defined on a set $\sS\subseteq \real^n$. Suppose that $\bx^\star\in\text{int}(\sS)$, i.e., in the interior point of the set, is a local optimum point and that all the partial derivatives (Definition <ref>, p. definition:partial_deri) of $f$ exist at $\bx^\star$. Then $\nabla f(\bx^\star)=\bzero$, i.e., the gradient vanishes at all local optimum points. (Note that, this optimality condition is a necessary condition; however, there could be vanished points which are not local maximum or minimum point.)

*  Global minimum point. Let function $f$ be a twice continuously differentiable function defined over $\real^n$. Suppose that the Hessian $\nabla^2f(\bx) \geq 0$ for any $\bx\in\real^n$. Then $\bx^\star$ is a global minimum point of $f$ if $\nabla f(\bx^\star)=\bzero$. Hint: use linear approximation theorem in Theorem <ref> (p. theorem:linear_approx).

* Two-sided matrix least squares. Let $\bB$ be an $M\times K$ matrix and $\bC$ be a $P\times N$ matrix. Find the $K\times P$ matrix $\bX$ such that $L(\bX)=\norm{\bA - \bB\bX\bC}_F^2$ is minimized, where $\bA\in\real^{M\times N}$ is known.

* Derive the derivative of $L$ with respect to $\bX$ and the optimality conditions.

* Show that one possible solution to the optimality conditions is $\bX=\bB^+\bA\bC^+$, where $\bB^+$ and $\bC^+$ are the pseudo-inverses of $\bB$ and $\bC$, respectively.

CHAPTER: NONNEGATIVE MATRIX FACTORIZATION (NMF)

Decomposition: NMF
Nonnegativity constraint

§ NONNEGATIVE MATRIX FACTORIZATION

In the era of big data, extracting meaningful patterns and latent structures from high-dimensional data sets has become a central challenge in various scientific and technological domains. Nonnegative matrix factorization (NMF) has emerged as a powerful and interpretable tool for dimensionality reduction, feature extraction, and discovering latent structures within complex data.
Early consideration of the NMF problem is due to Paatero and Tapper, 1994, Cohen and Rothblum, 1993, where they called it positive matrix factorization.
While Lee and Seung, 2001 later made the problem famous by the multiplicative update.

Following  the matrix factorization via the ALS, we now consider algorithms for solving the nonnegative matrix factorization (NMF) problem:

* Given a nonnegative matrix $\bA\in \real_+^{M\times N}$ with rank $R$, find nonnegative matrix factors $\bW\in \real_+^{M\times K}$ and $\bZ\in \real_+^{K\times N}$ such that:

\bA\approx\bW\bZ.

The nonnegativity constraint imposes an inherent sparsity, enabling the factorization to capture additive features, which is especially advantageous in applications where parts-based representations are meaningful.
This constraint make it find applications in diverse fields such as image processing, document analysis, and bioinformatics, where the identified components often correspond to distinct parts or features.
For example, in image processing, NMF has proven valuable for tasks such as  object detection,  image segmentation, and facial recognition. The decomposition into nonnegative components aligns with the intuitive notion that images are composed of identifiable parts.

To measure the quality of the  approximation, we evaluate the loss by computing  the Frobenius norm of the difference between the original matrix and the approximation:

L(\bW,\bZ) = \norm{\bW\bZ-\bA}^2.

See also the applications of the NMF in the survey paper Berry et al., 2007. When we want to find two nonnegative matrices $\bW\in\real^{M\times R}_+$ and $\bZ\in\real_+^{R\times N}$ such that $\bA=\bW\bZ$, the problem is known as the Exact NMF of $\bA$ of size $R$. Exact NMF is NP-hard [Gillis, 2020].
Thus, we only consider the approximation of NMF here.

Alternating update
Kullback-Leibler divergence
Multiplicative update

§ NMF VIA MULTIPLICATIVE UPDATE (MU)

We consider the NMF via an alternating update. The hidden features in $\bW$ and $\bZ$ are modeled as nonnegative vectors in low-dimensional space. These latent vectors are randomly initialized and iteratively updated via an alternating multiplicative update rule to minimize the Kullback-Leibler divergence between the observed and modeled matrices.
Following Section <ref> (p. section:als-netflix), we consider the low-rank with $K$ problem; given $\bW\in \real_+^{M\times K}$, we want to update $\bZ\in \real_+^{K\times N}$, the gradient with respect to $\bZ$ is given by Equation (<ref>) (p. equation:givenw-update-z-allgd):

\begin{aligned}
\frac{\partial L(\bZ|\bW)}{\partial \bZ} =2  \bW^\top(\bW\bZ-\bA) \in \real^{K\times N}.
\end{aligned}

Applying the gradient descent idea discussed in Section <ref> (p. section:als-gradie-descent), the trivial update for $\bZ$ can be achieved by

(\text{GD on $\bZ$})\gap \bZ \leftarrow \bZ - \eta \left(\frac{\partial L(\bZ|\bW)}{\partial \bZ}\right)=\bZ - \eta \left(2  \bW^\top\bW\bZ-2\bW^\top\bA\right),

where $\eta$ represents a small positive step size. Now if we impose a different step size for each entry of $\bZ$ and incorporate the constant 2 into the step size, the update can be obtained by

(\text{GD$^\prime$ on $\bZ$})\gap
\begin{aligned}
z_{kn} &\leftarrow z_{kn} - \frac{\eta_{kn}}{2} \left(\frac{\partial L(\bZ|\bW)}{\partial \bZ}\right)_{kn}\\
&=z_{kn} - \eta_{kn}(\bW^\top\bW\bZ-\bW^\top\bA)_{kn}, \gap k\in \{1,\ldots,K\}, n\in\{1,\ldots,N\},
\end{aligned}

where $z_{kn}$ denotes the $(k,n)$-th entry of $\bZ$.
To proceed, we further rescale the step size:

\eta_{kn} = \frac{z_{kn}}{(\bW^\top\bW\bZ)_{kn}}.

Then we obtain the update rule:

(\text{MU on $\bZ$})\gap
z_{kn} \leftarrow z_{kn}\frac{(\bW^\top\bA)_{kn}}{(\bW^\top\bW\bZ)_{kn}}, \gap k\in \{1,\ldots,K\}, n\in\{1,\ldots,N\},

which is known as the multiplicative update (MU), and is first developed in Lee and Seung, 2001 and further discussed in Pauca et al., 2006. Analogously, the multiplicative update for $\bW$ can be obtained by

\begin{equation}\label{equation:multi-update-w}
(\text{MU on $\bW$})\,\,
w_{mk} \leftarrow w_{mk} \frac{(\bA\bZ^\top)_{mk}}{(\bW\bZ\bZ^\top)_{mk}}, \,\, m\in\{1,\ldots,M\}, k\in\{1,\ldots,K\}.
\end{equation}

The loss $L(\bW,\bZ)=\norm{\bW\bZ-\bA}^2$ remains non-increasing under the following multiplicative update rules:

\left\{
\begin{aligned}
z_{kn} &\leftarrow z_{kn}\frac{(\bW^\top\bA)_{kn}}{(\bW^\top\bW\bZ)_{kn}}, \gap k\in \{1,\ldots,K\}, n\in\{1,\ldots,N\};\\
w_{mk} &\leftarrow w_{mk} \frac{(\bA\bZ^\top)_{mk}}{(\bW\bZ\bZ^\top)_{mk}}, \gap  m\in\{1,\ldots,M\}, k\in\{1,\ldots,K\}.
\end{aligned}
\right.

We refer the proof of the  theorem above to Lee and Seung, 2001.
Clearly, the approximations $\bW$ and $\bZ$ remain nonnegative during the updates.
It is generally better to update $\bW$ and $\bZ$ “simultaneously” rather than “sequentially," i.e., updating each matrix completely before the other. In this case, after updating a row of $\bZ$, we update the corresponding column of $\bW$.
In the implementation, it is advisable to introduce a small positive quantity, say the square root of the machine precision, to the denominators in the approximations of $\bW$ and $\bZ$
at each iteration.
And a trivial value like $\epsilon=10^{-9}$  suffices. The full procedure is shown in Algorithm <ref>.

Machine precision
NMF via Multiplicative Updates

[1]
Matrix $\bA\in \real^{M\times N}$;
Initialize $\bW\in \real^{M\times K}$, $\bZ\in \real^{K\times N}$ randomly with nonnegative entries;
Choose a stop criterion on the approximation error $\delta$;
Choose maximal number of iterations $C$;
$iter=0$; Count for the number of iterations
$\norm{\bA- (\bW\bZ)}^2>\delta $ and $iter<C$
$k=1$ to $K$
$n=1$ to $N$ update $k$-th row of $\bZ$
$z_{kn} \leftarrow z_{kn}\frac{(\bW^\top\bA)_{kn}}{(\bW^\top\bW\bZ)_{kn}+\textcolor{mylightbluetext}{\epsilon}}$;
$m=1$ to $M$ update $k$-th column of $\bW$
$w_{mk} \leftarrow w_{mk} \frac{(\bA\bZ^\top)_{mk}}{(\bW\bZ\bZ^\top)_{mk}+\textcolor{mylightbluetext}{\epsilon}}$;
Output $\bW,\bZ$;

§ REGULARIZATION

NMF via Regularized Multiplicative Updates

[1]
Matrix $\bA\in \real^{M\times N}$;
Initialize $\bW\in \real^{M\times K}$, $\bZ\in \real^{K\times N}$ randomly with nonnegative entries;
Choose a stop criterion on the approximation error $\delta$;
Choose maximal number of iterations $C$;
Choose regularization parameter $\lambda_z, \lambda_w$;
$iter=0$;Count for the number of iterations
$\norm{\bA- (\bW\bZ)}^2>\delta $ and $iter<C$
$k=1$ to $K$
$n=1$ to $N$ update $k$-th row of $\bZ$
$z_{kn} \leftarrow z_{kn}\frac{(\bW^\top\bA)_{kn}-
\textcolor{mylightbluetext}{\lambda_z z_{kn}}}{(\bW^\top\bW\bZ)_{kn}+\textcolor{mylightbluetext}{\epsilon}}$;
$m=1$ to $M$ update $k$-th column of $\bW$
$w_{mk} \leftarrow w_{mk} \frac{(\bA\bZ^\top)_{mk}	-\textcolor{mylightbluetext}{\lambda_w w_{mk}}}{(\bW\bZ\bZ^\top)_{mk}+\textcolor{mylightbluetext}{\epsilon}}$;
Output $\bW,\bZ$;

Similar to the ALS with regularization discussed in Section <ref> (p. section:regularization-extention-general), recall that the regularization can help extend the applicability of ALS to general matrices.
Additionally, a regularization term can be incorporated into the NMF framework to enhance its performance:

L(\bW,\bZ)  =\norm{\bW\bZ-\bA}^2 +\lambda_w \norm{\bW}^2 + \lambda_z \norm{\bZ}^2, \qquad \lambda_w>0, \lambda_z>0,

where the employed matrix norm is still the Frobenius norm. The gradient with respect to $\bZ$ given $\bW$ is the same as that in Equation (<ref>) (p. equation:als-regulari-gradien):

\begin{aligned}
\frac{\partial L(\bZ|\bW)}{\partial \bZ} =2\bW^\top(\bW\bZ-\bA) + \textcolor{mylightbluetext}{2\lambda_z\bZ}  \in \real^{K\times N}.
\end{aligned}

The gradient descent update can be obtained by

(\text{GD on }\bZ) \gap \bZ \leftarrow \bZ - \eta \left(\frac{\partial L(\bZ|\bW)}{\partial \bZ}\right)=\bZ - \eta \left(2  \bW^\top\bW\bZ-2\bW^\top\bA+\textcolor{mylightbluetext}{2\lambda_z\bZ}\right),

Analogously, if we assume a different step size for each entry of $\bZ$ and incorporate the constant 2 into the step size, the update can be obtained by

(\text{GD$^\prime$ on $\bZ$})\gap
\begin{aligned}
z_{kn} &\leftarrow z_{kn} - \frac{\eta_{kn}}{2} \left(\frac{\partial L(\bZ|\bW)}{\partial \bZ}\right)_{kn}\\
&=z_{kn} - \eta_{kn}(\bW^\top\bW\bZ-\bW^\top\bA+\textcolor{mylightbluetext}{\lambda_z\bZ})_{kn},
\,  k\in \{1,\ldots,K\}, n\in\{1,\ldots,N\}.
\end{aligned}

Now if we rescale the step size:

\eta_{kn} = \frac{z_{kn}}{(\bW^\top\bW\bZ)_{kn}},

then we obtain the multiplicative update rule:

\begin{aligned}
(\text{MU on $\bZ$})\gap
z_{kn} &\leftarrow z_{kn}\frac{(\bW^\top\bA)_{kn}-
\textcolor{mylightbluetext}{\lambda_z z_{kn}}
\gap  k\in \{1,\ldots,K\}, n\in\{1,\ldots,N\}.
\end{aligned}

Similarly, the multiplicative update on $\bW$ can be obtained by

\begin{aligned}
(\text{MU on $\bW$})\gap
w_{mk} &\leftarrow w_{mk} \frac{(\bA\bZ^\top)_{mk}
-\textcolor{mylightbluetext}{\lambda_w w_{mk}}
\gap   m\in\{1,\ldots,M\}, k\in\{1,\ldots,K\}.
\end{aligned}

The procedure is then formulated in Algorithm <ref>.

A nonnegative matrix factorization $\bA\approx \bW\bZ$ can be applied for clustering algorithm.
Specifically, the data vector $\ba_j$ is assigned to cluster $i$ if $z_{ij}$ is the largest element in column $j$ of $\bZ$ [Brunet et al., 2004, Gao and Church, 2005].

In the context of collaborative filtering, it is recognized  that the NMF via multiplicative update can result in overfitting  despite favorable convergence results.
The overfitting can be partially mitigated through regularization, but its out-of-sample performance remains low.
Bayesian optimization through the use of generative models, on the other hand, can effectively prevent overfitting for nonnegative matrix factorization [Brouwer et al., 2017, Lu and Ye, 2022]; see Chapter <ref> (p. chapter:bnmf).

Readers are advised to consult the survey of Berry et al., 2007 for other issues in the NMF.

Bayesian inference
Bayesian optimization
Bayesian matrix decomposition

§ INITIALIZATION

A significant challenge in NMF is the absence of guaranteed convergence to a global minimum.
It often happens that convergence is slow and a suboptimal approximation is reached.
In the preceding discussion, we initialized $\bW$ and $\bZ$ randomly.
To mitigate this issue, there are also alternative strategies designed to obtain better initial estimates in the hope of converging more rapidly to a good solution [Boutsidis and Gallopoulos, 2008, Gillis, 2014]. We sketch the methods as follows for reference:

* Clustering techniques. Apply some clustering methods to the columns of $\bA$, set the cluster means of the top $K$ clusters as the columns of $\bW$, and initialize $\bZ$ as a proper scaling of
the cluster indicator matrix (that is, $z_{kn}\neq 0$ indicates $\ba_n$ belongs to the $k$-th cluster);

* Subset selection. Pick $K$ columns of $\bA$, and set those as the initial columns for $\bW$. And analogously, $K$ rows of $\bA$ are selected to form the rows of $\bZ$;

* SVD-based approach. Suppose the optimal rank-$K$ approximation  of $\bA$ via SVD is $\bA=\sum_{i=1}^{K}\sigma_i\bu_i\bv_i^\top$, where each factor $\sigma_i\bu_i\bv_i^\top$ is a rank-one matrix with possible negative values in $\bu_i$ and $\bv_i$, and nonnegative $\sigma_i$. Denote $[x]_+=\max(x, 0)$, we notice

\bu_i\bv_i^\top = [\bu_i]_+[\bv_i]_+^\top+[-\bu_i]_+[-\bv_i]_+^\top-[-\bu_i]_+[\bv_i]_+^\top-[\bu_i]_+[-\bv_i]_+^\top,

where the first two rank-one factors in this decomposition are nonnegative.
Then, either $[\bu_i]_+[\bv_i]_+^\top$ or $[-\bu_i]_+[-\bv_i]_+^\top$ can be selected to replace the factor $\bu_i\bv_i^\top$. Boutsidis and Gallopoulos, 2008 suggests to replace each rank-one factor in $\sum_{i=1}^{K}\sigma_i\bu_i\bv_i^\top$ with  either $[\bu_i]_+[\bv_i]_+^\top$ or $[-\bu_i]_+[-\bv_i]_+^\top$, selecting the one with larger norm and scaling it properly.
[Note here we consider general matrix $\bA$. If $\bA$ is nonnegative, then $\bu_i$ and $\bv_i$ are nonnegative as well.]
In other words, if we select $[\bu_i]_+[\bv_i]_+^\top$, then $\sigma_i\cdot [\bu_i]_+  $ can be initialized as the $i$-th column of $\bW$, and $[\bv_i]_+^\top$ can be chosen as the $i$-th row of $\bZ$.

However, these techniques are not guaranteed to yield better performance theoretically.
We recommend referring to the aforementioned papers for more information.

Implicit hierarchy

§ MOVIE RECOMMENDER CONTEXT

Both the NMF and the ALS methods approximate the matrix and reconstruct the entries in the matrix with a set of basis vectors.
The basis  in the NMF is composed of vectors with nonnegative elements while the basis vectors in the ALS can have positive or negative values.
The key distinction lies in how the approximation is carried out.
In  NMF, each vector is reconstructed  as a positive summation of the basis vectors with a “relative" small component in the direction of each basis vector.
Whereas, in the ALS approximation, the data is modeled as a linear combination of the basis such that we can add or subtract vectors as needed; and the components in the direction of each basis vector can be large positive values or negative values. Therefore, depending on the application, one or the other factorization can be utilized to describe the data with different meanings.

In the context of a movie recommender system,  the rows of $\bW$ represent the hidden features of  movies, while the columns of $\bZ$ represent the hidden features of users.
In the NMF method, we can say that a movie is 0.5 comedy, 0.002 action, and 0.09 romantic. However, in the ALS approach, we can get combinations such as 4 comedy, $-0.05$ comedy, and $-3$ drama, i.e., a positive or negative component on that feature.

The ALS and NMF are similar in the sense that the importance of each basis vector is not ranked in a hierarchical manner. Whereas, the key difference between the ALS (or NMF) and the SVD is that, in the SVD, the importance of each vector in the basis is relative to the value of the singular value associated with that vector.
In the SVD representation of $\bA=\sum_{i=1}^{r}\sigma_i\bu_i\bv_i^\top$,
this usually means that the reconstruction $\sigma_1\bu_1\bv_1^\top$ via the first set of basis vectors dominates and is the most used set to reconstruct data, followed by the second set, and so on.
So the basis in the SVD has an implicit hierarchy and that doesn't happen in the ALS or the NMF approaches.
Recall the low-rank approximation of the flag image in Section <ref> (p. section:als-low-flag), where we find the second component $\bw_2\bz_2^\top$ via the ALS in Figure <ref> (p. fig:als52) plays an important role in the reconstruction of the original figure.
Conversely, the second component $\sigma_2\bu_2\bv_2^\top$ via the SVD in Figure <ref> (p. fig:svd22) plays a less important role in the reconstruction.

* * Suppose $\bC=P(\bA)$ is defined as the matrix obtained by setting each entry of $\bC$ as the absolute value of $\bA\in\real^{n\times n}$. Suppose further that $\bB - P(\bA)\in\real_+^{n\times n}$, prove that

\rho(\bA) \leq \rho(P(\bA)) \leq \rho(\bB),

where $\rho(\bX)$ represents the spectral radius of matrix $\bX$ (Definition <ref>, p. definition:spectrum).

* * Given $\bA\in\real_+^{n\times n}$, show that

\begin{aligned}
\mathop{\min}_{1\leq i \leq n} \sum_{j=1}^{n} a_{ij}
&\leq \rho(\bA)
\leq
\mathop{\max}_{1\leq i \leq n} \sum_{j=1}^{n} a_{ij}; \\
\mathop{\min}_{1\leq j \leq n} \sum_{i=1}^{n} a_{ij}
&\leq \rho(\bA)
\leq
\mathop{\max}_{1\leq j \leq n} \sum_{i=1}^{n} a_{ij}. \\
\end{aligned}

Bayesian Matrix Decomposition

CHAPTER: BAYESIAN REAL MATRIX FACTORIZATION

§ INTRODUCTION

The explosion of data from advancements in sensor technology and computer hardware has presented new challenges for data analysis.
The large volume of data often contains noise and other distortions, requiring pre-processing for deductive science to be applied.
For example, signals received by antenna arrays often are contaminated by noise and other degradations.
To effectively analyze the data, it is necessary to reconstruct or represent it in a way that reduces inaccuracies while maintaining certain feasibility conditions.

Additionally, in many cases, the data collected from complex systems is the result of multiple interrelated variables acting in unison. When these variables are not well defined, the information contained in the original data can be overlapping and unclear. By creating a reduced system model, we can achieve a level of accuracy that is close to the original system. The common approach in removing noise, reducing the model, and reconstructing feasibility, is to replace the original data with a lower-dimensional representation obtained through subspace approximation.
Therefore, low-rank approximations or low-rank matrix decompositions play a crucial role in
a wide range of applications.

Low-rank matrix decomposition is a powerful technique   used in machine learning and data mining  to represent a given matrix as the product of two or more matrices with lower dimensions.
It is used to capture the essential structure of a matrix while ignoring noise and redundancies. The most common methods for low-rank matrix decomposition include singular value decomposition (SVD), principal component analysis (PCA), and multiplicative update nonnegative matrix factorization (NMF).

Bayesian low-rank decomposition is a variant of low-rank matrix decomposition that incorporates Bayesian modeling. It models the observed data as a low-rank matrix, where the low-rank approximation is assumed to be generated from a prior distribution.
This incorporation  allows for the inclusion of prior knowledge and uncertainty about the low-rank matrix into the decomposition.
The use of priors can also lead to more interpretable results by providing a probabilistic representation of the uncertainty associated with the factor matrices.
In addition, Bayesian methods allow for the quantification of uncertainty in the results, providing a measure of the confidence in the estimated factor matrices.
Therefore, it can help to mitigate overfitting and produce more robust results, establishing it as a powerful method for modeling both predictive and explanatory data.

Given an observed data set, represented as an $M\times N$ matrix of results, $\bA$, where rows represent the number of observations and columns represent the variables of interest.
Following Chapter <ref> (p. section:als),
the real matrix factorization (RMF) problem, a common bilinear decomposition problem, can be stated as $\bA=\bW\bZ+\bE$,  where $\bA=[\ba_1, \ba_2, \ldots, \ba_N]\in \real^{M\times N}$ is approximately factorized into an $M\times K$  matrix $\bW\in \real^{M\times K}$ and a $K\times N$ matrix  $\bZ\in \real^{K\times N}$. The data set $\bA$ needs not to be complete such that the indices of observed entries can be represented by the mask matrix $\bM\in \real^{M\times N}$, where a value of 1 indicates the entry is observed and a value of 0 indicates the entry is missing.
Matrices $\bW$ and $\bZ$ represent the values of explanatory variables which, when multiplied, give a predictor of the values in $\bA$.
If entries in $\bA$ are missing, then $\bW$ and $\bZ$ can be used to give predictions of their values.
When one of $\bW$ and $\bZ$ is observed, the factorization becomes a regression problem.

The factorization of the original data matrix $\bA$ is achieved by finding two such  real matrices, one representing the basis or dictionary components and the other representing the activations or coefficients.
Let $\bz_n$ denote the $n$-th column of $\bZ$. Then the matrix multiplication of $\bW\bZ$ can be implemented as computing the column vectors of $\bA$ as linear combinations of the columns in $\bW$ using coefficients supplied by columns of $\bZ$:

\ba_n = \bW\bz_n.

In the Netflix context, the value of $a_{mn}$, the $(m,n)$-th element of $\bA$, denotes the rating of the $n$-th user for the $m$-th movie (the larger the user more likes the movie). Then $\bw_m$ can represent the hidden features of the $m$-th movie, and $\bz_n$ contains the features of user $n$ (Section <ref>, p. section:als-vector-product).

To simplify the problem, let us assume that there are no missing entries firstly.
Project data vectors $\ba_n$ to a smaller dimension $\bz_n \in \real^K$  with $K<M$, such that the reconstruction error measured by Frobenius norm is minimized (assume $K$ is known):

\begin{equation}\label{equation:als-per-example-loss_real}
\mathop{\min}_{\bW,\bZ}  \sum_{n=1}^N \sum_{m=1}^{M} \left(a_{mn} - \bw_m^\top\bz_n\right)^2,
\end{equation}

where $\bW=[\bw_1^\top; \bw_2^\top; \ldots; \bw_M^\top]\in \real^{M\times K}$ and $\bZ=[\bz_1, \bz_2, \ldots, \bz_N] \in \real^{K\times N}$ containing $\bw_m$'s and $\bz_n$'s as rows and columns, respectively [Note in some contexts, $\bZ$ represents an $N\times K$ matrix such that $\bA$ is decomposed into $\bA=\bW\bZ\textcolor{mylightbluetext}{^\top} +\bE$.]. The loss form in Equation (<ref>) is known as the per-example loss. It can be equivalently written as

\begin{equation}\label{equation:frob_loss_brmf}
L(\bW,\bZ) = \sum_{n=1}^N \sum_{m=1}^{M} \left(a_{mn} - \bw_m^\top\bz_n\right)^2
= \norm{\bW\bZ-\bA}^2
=\trace\left\{(\bW\bZ-\bA)^\top (\bW\bZ-\bA)\right\},
\end{equation}

where $\trace(\cdot)$ represents the trace of the quantity in the brackets.
This matrix factorization problem is similar to a standard “inverse" problem, except that in the “inverse" problem, one of the factored components is known, and thus ordinary least squares or similar methods can be applied to find the other component, which minimizes the residuals between the reconstruction and the data.
When neither $\bW$ nor $\bZ$ is known, the factorization problem is difficult even if the latent dimension $K$ is only 2 or 3.
We sample the space of potential solutions using the Markov chain Monte Carlo  (MCMC) procedure to determine its properties because there are a vast number of possible solutions and no analytical approach to identify them.

We have discussed the Bayesian approach in Section <ref> (p. section:bayes_approach). In the Bayesian matrix factorization context, the model leads to the specific form of Bayes' equation,

\begin{equation}\label{equation:bmf_bayes}
p(\bW,\bZ \mid \bA) \propto p(\bA \mid \bW, \bZ) \times p(\bW, \bZ),
\end{equation}

where $p(\bW, \bZ)$ captures the prior beliefs encoding the knowledge of the solution independent of the data, and $p(\bA\mid \bW, \bZ)$ denotes the likelihood comparing the model to the data.

Terminology. There are three types of choices we make that determine
the specific type of matrix decomposition model we use, namely, the likelihood function, the priors we place over the factored matrices $\bW$ and $\bZ$, and whether we use any further hierarchical priors.
We will call the model by the density function in the order of the types of likelihood and priors. For example, if the likelihood function for the model is chosen to be a Gaussian density, and the two prior density functions are selected to be exponential density and Gaussian density functions, respectively, then the model will be denoted by the Gaussian Exponential-Gaussian (GEG) model. Sometimes, we will put a hyperprior over the parameters of the prior density functions, e.g., we put a Gamma prior over the exponential density, then it will further be termed as a Gaussian Exponential-Gaussian Gamma (GEGA) model (“A" is short for the Gamma density to avoid confusion with the Gaussian density). Table <ref> summarizes the Bayesian models for real matrix factorization in this chapter.

Name     Likelihood    Prior $\bW$     Prior $\bZ$      Hierarchical prior

GGG     $\normal(a_{mn}|\bw_m^\top\bz_n, \sigma^2)$
$\normal(w_{mk}|0, (\lambda_{mk}^W)^{-1})$     $\normal(z_{kn}|0, (\lambda_{kn}^Z)^{-1})$
GGGM     $\normal(a_{mn}|\bw_m^\top\bz_n, \sigma^2)$
$\normal(\bw_m|\bzero, \lambda^{-1}\bI)$     $\normal(\bz_{n}|\bzero, \lambda^{-1}\bI)$
GGGA    $\normal(a_{mn}|\bw_m^\top\bz_n, \sigma^2)$
$\normal(w_{mk}|0, (\lambda_k)^{-1})$     $\normal(z_{kn}|0, (\lambda_k)^{-1})$     $\gammadist(\lambda_k|\alpha_\lambda, \beta_\lambda)$
GGGW    $\normal(a_{mn}|\bw_m^\top\bz_n, \sigma^2)$
$\normal(\bw_m|\bmu_w, \bSigma_w)$    $\normal(\bz_n|\bmu_z, \bSigma_z)$
$\{\bmu_w, \bSigma_w\}$, $\{\bmu_z, \bSigma_z\}\sim$

$\niw(\bmm_0, \kappa_0, \nu_0, \bS_0)$
GVG     $\normal(a_{mn}|\bw_m^\top\bz_n, \sigma^2)$

$\exp\{-\gamma \bW^\top\bW\}$
$\normal(z_{kn}|0, (\lambda_{kn}^Z)^{-1})$

Overview of Bayesian real matrix factorization models.

Decomposition: GGG

§ ALL GAUSSIAN (GGG) MODEL AND MARKOV BLANKET

Graphical model representation of the GGG model.  Green circles denote prior variables, orange circles represent observed and latent variables, and plates represent repeated variables. The slash “/" in the variable represents “or."

The all Gaussian (GGG) model is perhaps the simplest approach for Bayesian RMF,  where Gaussian priors are applied over factored matrices [Salakhutdinov and Mnih, 2008, Gönen, 2012, Virtanen et al., 2011, Virtanen et al., 2012]. The model involves using Gaussian likelihood and Gaussian priors.

We view the data $\bA$ as being produced according to the probabilistic generative process shown in Figure <ref>. The observed $(m,n)$-th data  entry $a_{mn}$ of matrix $\bA$ is modeled using a Gaussian likelihood function with variance $\sigma^2$ and a mean determined by the latent decomposition $\bw_m^\top\bz_n$ (Equation (<ref>)):

\begin{equation}\label{equation:ggg_data_entry_likelihood}
p(a_{mn} \mid \bw_m^\top\bz_n, \sigma^2) = \normal(a_{mn}\mid\bw_m^\top\bz_n, \sigma^2).
\end{equation}

This is equivalent to assuming the residuals, $e_{mn}$, are i.i.d. drawn from a zero-mean normal with variance $\sigma^2$, which gives rise to the following likelihood function:

\begin{equation}\label{equation:ggg_likelihood}
\begin{aligned}
p(\bA\mid \btheta) &= \prod_{m,n=1}^{M,N} \normal \left(a_{mn}\mid (\bW\bZ)_{mn}, \sigma^2 \right)\\
&= \prod_{m,n=1}^{M,N} \normal \left(a_{mn}\mid (\bW\bZ)_{mn}, \tau^{-1} \right),
\end{aligned}
\end{equation}

where $\btheta=\{\bW,\bZ,\sigma^2\}$ denotes all parameters in the model, $\sigma^2$ is the variance, $\tau^{-1}=\sigma^2$ is the precision, and

\normal(x\mid \mu,\sigma^2)=\frac{1}{(2\pi\sigma^2)^{1/2}} \exp \left\{ -\frac{1}{2\sigma^2 } (x-\mu)^2\right\}
=\sqrt{\frac{\tau}{2\pi}}\exp \left\{ -\frac{\tau}{2}(x-\mu)^2 \right\}

is the normal density (Definition <ref>, p. definition:gaussian_distribution).

Prior. We assume $\bW$ and $\bZ$ are independently Gaussian distributed with precisions $\lambda_{mk}^W$ and $\lambda_{kn}^Z$, respectively,

\begin{equation}\label{equation:ggg_prior_density_gaussian}
\begin{aligned}
w_{mk} &\sim \normal(w_{mk}\mid 0, (\lambda_{mk}^W)^{-1}),
\gap
&z_{kn}\sim&  \normal(z_{kn}\mid 0, (\lambda_{kn}^Z)^{-1});\\
p(\bW) &=\prod_{m,k=1}^{M,K} \normal(w_{mk}\mid 0, (\lambda_{mk}^W)^{-1}),
\gap
&p(\bZ) =&\prod_{k,n=1}^{K,N} \normal(z_{kn}\mid 0,( \lambda_{kn}^Z)^{-1}).
\end{aligned}
\end{equation}

The prior for the noise variance $\sigma^2$ is chosen as an inverse-Gamma density with shape ${\alpha_\sigma}$ and scale ${\beta_\sigma}$ (Definition <ref>, p. definition:inverse_gamma_distribution),

p(\sigma^2)= \inversegammadist(\sigma^2\mid \alpha_\sigma, \beta_\sigma) = \frac{{\beta_\sigma}^{\alpha_\sigma}}{\Gamma({\alpha_\sigma})} (\sigma^2)^{-\alpha_\sigma-1} \exp\left( -\frac{{\beta_\sigma}}{\sigma^2} \right).

By Bayes' rule (Equation (<ref>), p. equation:posterior_abstract_for_mcmc), the posterior is proportional to the product of likelihood and prior, it can be maximized to yield an estimate of $\bW$ and $\bZ$.

Inverse-Gamma distribution

Markov blanket.

The most widely used posterior inference methods in Bayesian inference models are Markov chain Monte Carlo (MCMC) methods as described in Section <ref> (p. sec:monte_carlo_methods). The concept behind MCMC methods is to define a Markov chain on the hidden variables that have the posterior as its equilibrium distribution [Andrieu et al., 2003]. By drawing samples from this Markov chain, one eventually obtains samples from the posterior distribution. A simple form of MCMC sampling is  Gibbs sampling, where the Markov chain is constructed by sampling the conditional distribution of each hidden variable given the values of  other hidden variables and the observations.
Gibbs sampling is widely used when these conditional distributions can be sampled from easily.

The Markov blanket of a directed acyclic graphical (DAG) model. In a Bayesian network, the Markov blanket of node $A$ includes its parents, children, and the other parents of all of its children. That is, the nodes in the cycle are in the Markov blanket of node $A$. The figure is due to wikipedia page of Markov blanket.

To do Gibbs sampling, we need to derive the conditional posterior distributions for each parameter conditioned on all the other parameters $p(\theta_i \mid  \btheta_{-i}, \mathcalX)$, where $\mathcalX$ is again the set of  data points (here, the observed matrix $\bA$), and $\theta_i$'s are the variables for which we want to sample the distributions. But for a graphical model, this conditional distribution is a function only of the nodes in the Markov blanket.
For the GGG model shown in Figure <ref>, which is  a directed acyclic graphical (DAG) model, the Markov blanket of a node includes the parents, the children, and the coparents [Jordan and Bishop, 2004], as shown in Figure <ref>. The Markov blanket of node A is all nodes in the cycle.

An example for the Markov blanket. The idea of the Markov blanket might be mysterious at first glance. Suppose we want to sample the $(m,k)$-th element $w_{mk}$ of $\bW$ for its distribution.
From Figure <ref>, we find its parents, children, and coparents are $\{\lambda_{mk}^W\}$, $\{a_{mn}\}$, and $\{\sigma^2, \bZ, \bW_{-mk}\}$, respectively. Therefore, the conditional distribution of $w_{mk}$ only depends on the three pairs of parameters:

p(w_{mk} \mid  -)  = p(w_{mk}\mid  \bA,\bW_{-mk},\bZ,  \sigma^2,  \lambda_{mk}^W).

More specifically, from this graphical representation, we can find the Markov blanket for each parameter in the GGG model, and then figure out their conditional posterior distributions to be derived:

\begin{align}
p(w_{mk} \mid  -)
&= p(w_{mk}\mid  \bA,   \bW_{-mk}, \bZ,  \sigma^2, \lambda_{mk}^W), \label{equation:ggg_1} \\
p(z_{kn} \mid  -)
&= p(z_{kn}\mid  \bA,  \bW, \bZ_{-kn}, \sigma^2,\lambda_{kn}^Z), \label{equation:ggg_2}\\
p(\sigma^2 \mid  -)
&= p(\sigma^2\mid   \bA,\bW, \bZ, \alpha_\sigma, \beta_\sigma,)  \label{equation:ggg_3}.
\end{align}

We sequentially draw samples from the posterior of each parameter, conditioned on all other parameters.
It can be shown that the sequence of samples computed constitutes a Markov chain, for which the stationary distribution is the posterior of interest.
In other words, Gibbs sampler moves the chain forward by one step as follows:

* Sample the first component $w_{mk}$ for each observation from Equation (<ref>), which is known as the conditional distribution of $w_{mk}$;

* Sample the second component  $z_{kn}$ for each observation from Equation (<ref>), which is known as the conditional distribution of $z_{kn}$;

* Sample the parameter $\sigma^2$ from Equation (<ref>), which is known as the conditional distribution of data variance.

Given the observed matrix $\bA$, we want to estimate the conditional distribution of the latent structure $p(\bW,\bZ \mid \bA)$, termed the posterior.
This estimation is the key to matrix decomposition-based applications. For example, in the Netflix context, we estimate the posterior expectation of each user's hidden feature/preferences and each movie's hidden attributes to perform predictions of which  unconsumed movies each user will like.
In this book, we apply Gibbs sampling to perform  optimization since it tends to be very accurate at finding the true posterior. The advantage of the MCMC or Gibbs sampling methods lies in their ability to produce exact results asymptotically.
Other than this method, variational Bayesian inference can be an alternative way, but we shall not go into the details.
For RMF, following the Bayes' rule and MCMC, this means we need to be able to draw from distributions (by Markov blanket):

\begin{aligned}
p(w_{mk} &\mid  \bA, \bW_{-mk}, \bZ, \sigma^2, \lambda_{mk}^W), \\
p(z_{kn}&\mid \bA,  \bW, \bZ_{-kn},\sigma^2,   \lambda_{kn}^Z ), \\
p(\sigma^2 &\mid  \bA,\bW, \bZ, \alpha_\sigma, \beta_\sigma), \\
\end{aligned}

where $\bW_{-{mk}}$ denotes all elements of $\bW$ except $w_{mk}$, and $\bZ_{-kn}$ denotes all elements of $\bZ$ except $\bz_{kn}$.
Using Bayes' theorem,
the conditional density of $w_{mk}$ depends on its parents ($\lambda_{mk}^W$), children ($a_{mn}$), and coparents ($\tau$ or $\sigma^2$, $\bW_{-mk}, \bZ$) [See Figure <ref> and Section <ref>.]. And it can be obtained by

\begin{equation}\label{equation:ggg_poster_wmk1}
\begin{aligned}
%&\gap p(w_{mk}\mid  \sigma^2, \bW_{-mk}, \bZ, \blambda^W,\cancel{\blambda^Z}, \bA)=
&\gap p(w_{mk} \mid  \bA ,   \bW_{-mk}, \bZ, \sigma^2, \lambda_{mk}^W )
\propto p(\bA\mid  \bW, \bZ, \sigma^2) \times p(w_{mk}\mid  \lambda_{mk}^W)\\
&=\prod_{i,j=1}^{M,N} \normal \left(a_{ij}\mid  \bw_i^\top\bz_j, \sigma^2 \right)\times \normal(w_{mk}\mid 0, (\lambda_{mk}^W)^{-1}) \\
&\propto \exp\left\{   -\frac{1}{2\sigma^2}  \sum_{i,j=1}^{M,N}(a_{ij} - \bw_i^\top\bz_j  )^2\right\}  \times \exp\left\{ -\frac{w_{mk}^2}{2} \lambda_{mk}^W \right\} \\
&\propto \exp\left\{   -\frac{1}{2\sigma^2}  \sum_{j=1}^{N}(a_{mj} - \bw_m^\top\bz_j  )^2\right\}
\times \exp\left\{ -\frac{w_{mk}^2}{2}\lambda_{mk}^W \right\} \\
&\propto \exp\left\{   -\frac{1}{2\sigma^2}  \sum_{j=1}^{N}
\left[ w_{mk}^2z_{kj}^2 + 2w_{mk} z_{kj}\bigg(\sum_{i\neq k}^{K}w_{mi}z_{ij} - a_{mj}\bigg)  \right]
\right\}  \cdot \exp\left\{ -\frac{w_{mk}^2}{2}\lambda_{mk}^W \right\}\\
&\stackrel{\star}{\propto} \exp\left\{
\left(\frac{\sum_{j=1}^{N} z_{kj}^2 }{2\sigma^2} +\frac{\lambda_{mk}^W}{2} \right) }_{\textcolor{mylightbluetext}{1/(2\widetilde{\sigma_{mk}^{2}})}}
w_{mk}\underbrace{\left(  \frac{1}{\sigma^2} \sum_{j=1}^{N} z_{kj}\bigg( a_{mj} - \sum_{i\neq k}^{K}w_{mi}z_{ij}\bigg)  \right)}_{\textcolor{mylightbluetext}{\widetilde{\sigma_{mk}^{2}}^{-1} \widetilde{\mu_{mk}}}}
\right\}  \\
&\propto   \normal(w_{mk} \mid  \widetilde{\mu_{mk}}, \widetilde{\sigma_{mk}^{2}}),
\end{aligned}
\end{equation}

where the equality ($\star$) conforms to Equation (<ref>) (p. equation:gaussian_form_conform) such that it follows from the normal distribution with variance $\widetilde{\sigma_{mk}^{2}}$,

\begin{equation}\label{equation:ggg_posterior_variance}
\widetilde{\sigma_{mk}^{2}} = 1\bigg/ \left(	\frac{1 }{\sigma^2} \sum_{j=1}^{N} z_{kj}^2 +\lambda_{mk}^W\right)
\end{equation}

and  mean $\widetilde{\mu_{mk}}$,

\begin{equation}\label{equation:ggg_posterior_mean}
\widetilde{\mu_{mk}} = \frac{\widetilde{\sigma_{mk}^{2}}}{\sigma^2} \cdot  \sum_{j=1}^{N} z_{kj}\bigg( a_{mj} - \sum_{i\neq k}^{K}w_{mi}z_{ij}\bigg) .
\end{equation}

In this case, the posterior precision $1/ \widetilde{\sigma_{mk}^{2}}$ is the sum of prior precision $\lambda_{mk}^W$ and “data precision" $\frac{1 }{\sigma^2} \sum_{j=1}^{N} z_{kj}^2$.
This also shows, in the Netflix context, that the conditional distribution over the movie feature vector $\bw_m$, conditioned on the user features, observed rating matrix $\bA$, and the values of the hyper-parameters, follows a Gaussian density.

Similarly, due to symmetry,
a completely analogous derivation  for parameter $z_{kn}$, where $k\in\{1,2,\ldots, K\}$ and $n\in\{1,2,\ldots, N\}$, is used for the factor $\bZ$. See Problem <ref>.<ref>

The conditional density of $\sigma^2$ depends on its parents ($\alpha_\sigma$, $\beta_\sigma$), children ($\bA$), and coparents ($\bW$, $\bZ$). And it follows an inverse-Gamma distribution, as demonstrated by the conjugacy presented in  Equation (<ref>) (p. equation:inverse_gamma_conjugacy_general):

\begin{equation}\label{equation:ggg_posterior_sigma2}
\begin{aligned}
&p(\sigma^2 \mid \bA, \bW,\bZ, \alpha_\sigma, \beta_\sigma ) = \inversegammadist (\sigma^2\mid  \widetilde{\alpha_{\sigma}}, \widetilde{\beta_{\sigma}}),  \\
& \widetilde{\alpha_{\sigma}} = \frac{MN}{2} +{\alpha_\sigma}, \gap
\widetilde{\beta_{\sigma}}  =  \frac{1}{2} \sum_{m,n=1}^{M,N} (\bA-\bW\bZ)_{mn}^2 + {\beta_\sigma}.
\end{aligned}
\end{equation}

Missing entries. In many cases, e.g., in the Netflix context, some entries of $\bA$ are missing. Let  $\Omega$ denote the set containing all observed entries in data $\bA$. Denote further $\Omega_m = \{n \mid  (m,n)\in \Omega\}$, i.e., the observed entries in the $m$-th row; $\Omega_n = \{m \mid  (m,n)\in \Omega\}$, i.e., the observed entries in the $n$-th column. The   posterior density of $w_{mk}$ can be obtained by

\begin{equation}
w_{mk}\sim \normal(w_{mk} \mid  \widetilde{\mu_{mk}}, \widetilde{\sigma_{mk}^{2}}),
\end{equation}

\begin{equation}
\begin{aligned}
\widetilde{\sigma_{mk}^{2}} &= 1\bigg/ \left(	\frac{1}{\sigma^2}\sum_{\textcolor{mylightbluetext}{j\in\Omega_m}} z_{kj}^2 +\lambda_{mk}^W\right),\\
\widetilde{\mu_{mk}} &= \frac{\widetilde{\sigma_{mk}^{2}}}{\sigma^2} \cdot  \sum_{\textcolor{mylightbluetext}{j\in\Omega_m}} z_{kj}\bigg( a_{mj} - \sum_{i\neq k}^{K}w_{mi}z_{ij}\bigg) .
\end{aligned}
\end{equation}

In the following discussion, for the sake of simplicity, we only consider   the data matrix $\bA$ with full observations. Analogous results for instances involving missing entries can be obtained in a similar manner.

Decomposition: GGGM

GGG with shared prior (GGGM).
When we set $\lambda=\lambda_{mk}^W$ for all $m\in\{1,2,\ldots, M\}$ and $k\in\{1,2,\ldots,K\}$, the conditional posterior distributions we obtain in the Gibbs sampling algorithm can be written as a multivariate Gaussian density (Definition <ref>, p. definition:multivariate_gaussian). In this case, we can place a multivariate Gaussian prior over each row $\bw_m$ of $\bW$ and each column $\bz_n$ of $\bZ$:

\bw_m \sim \normal(\bw_m \mid  \bzero, \lambda^{-1}\bI), \gap
\bz_n \sim \normal(\bz_n \mid  \bzero, \lambda^{-1}\bI),

i.e., each of the $M$ or $N$ item factors follows a multivariate normal distribution.
Again by Bayes' rule, let $\bW_{-m}$ denote all elements of $\bW$ except the $m$-th row, the conditional posterior distribution we obtain in the Gibbs sampling algorithm can be obtained by

\begin{equation}\label{equation:gggm_wm_post}
\begin{aligned}
& \gap p(\bw_m \mid  \sigma^2, \bW_{-m}, \bZ, \lambda, \bA)
\propto p(\bA \mid  \bW, \bZ, \sigma^2) \times \normal(\bw_m \mid  \bzero ,\lambda^{-1}\bI) \\
&\propto \normal(\bA \mid  \bW\bZ, \sigma^2\bI) \times \normal(\bw_m \mid  \bzero ,\lambda^{-1}\bI) \\
&\propto \exp\left\{   -\frac{1}{2\sigma^2}  \sum_{j=1}^{N}(a_{mj} - \bw_m^\top\bz_j  )^2\right\}
\times \exp\left\{ -\frac{\lambda}{2} \bw_m^\top\bw_m \right\} \\
&\stackrel{\star}{\propto} \exp\left\{-\frac{1}{2} \bw_m^\top
\underbrace{\left[\lambda\bI + \frac{1}{\sigma^2}\sum_{j=1}^{N}\bz_j\bz_j^\top \right]}_{\textcolor{mylightbluetext}{\widetilde{\bSigma}^{-1}}}
\bw_m
+ \bw_m^\top
\underbrace{\frac{1}{\sigma^2} \sum_{j=1}^{N} a_{mj}\bz_j }_{\textcolor{mylightbluetext}{\widetilde{\bSigma}^{-1} \widetilde{\bmu} }}
\right\}
\propto \normal(\bw_m \mid  \widetilde{\bmu}, \widetilde{\bSigma}),
\end{aligned}
\end{equation}

where the equality ($\star$) conforms to Equation (<ref>) (p. equation:multi_gaussian_form_conform) such that it follows from the multivariate Gaussian distribution with covariance $\widetilde{\bSigma}$,

\widetilde{\bSigma}= \left[\lambda\bI + \frac{1}{\sigma^2}\sum_{j=1}^{N}\bz_j\bz_j^\top \right]^{-1}

and  mean vector $\widetilde{\bmu}$,

\widetilde{\bmu} = \frac{1}{\sigma^2}  \widetilde{\bSigma}\cdot \sum_{j=1}^{N} a_{mj}\bz_j.

In this case, the posterior precision matrix $\widetilde{\bSigma}^{-1}$ is the sum of the prior precision matrix $\lambda\bI$ and “data precision matrix" $\frac{1}{\sigma^2}\sum_{j=1}^{N}\bz_j\bz_j^\top$.
And due to symmetry, a analogous expression can be derived for each column $\bz_n$ of factor $\bZ$, where $n\in\{1,2,\ldots, N\}$.
See Problem <ref>.<ref>

Gibbs sampling.
Because conjugate priors for both the parameters and hyper-parameters are used in the Bayesian matrix factorization model, it is easy to sample from the conditional distributions derived from the posterior distribution.
By this Gibbs sampling method introduced in Section <ref> (p. section:gibbs-sampler), we can construct a Gibbs sampler for the GGG model as formulated in Algorithm <ref>.
Due to our choice of priors, we can sample from all conditional distributions directly using standard methods, which obviates slower sampling procedures such as rejection sampling.
And also in practice, all the parameters of Gaussian priors are set to be the same value $\lambda=\{\lambda_{mk}^W\}'s = \{\lambda_{kn}^Z\}'s$ for all $m,k,n$.
By default, uninformative hyper-parameters are $\alpha_\sigma=\beta_\sigma=1$, $\{\lambda_{mk}^W\} = \{\lambda_{kn}^Z\} = 0.1$.

Gibbs sampler for GGG model in one iteration (prior on variance $\sigma^2$ here, similarly for the precision $\tau$). The procedure presented
here may not be efficient but is explanatory. A more efficient one can be implemented in a
vectorized manner. By default, uninformative hyper-parameters are $\alpha_\sigma=\beta_\sigma=1$, $\{\lambda_{mk}^W\} = \{\lambda_{kn}^Z\} = 0.1$.

[1]
Choose initial $\alpha_\sigma, \beta_\sigma, \lambda_{mk}^W, \lambda_{kn}^Z$;
$k=1$ to $K$
$m=1$ to $M$
Sample $w_{mk}$ from $p(w_{mk}\mid  \bA ,   \bW_{-mk}, \bZ, \sigma^2, \lambda_{mk}^W )$;
Equation (<ref>)
$n=1$ to $N$
Sample $z_{kn}$ from $p(z_{kn}\mid \bA ,   \bW, \bZ_{-kn}, \sigma^2, \lambda_{kn}^Z)$;
Symmetry of Eq. (<ref>)
Sample $\sigma^2$ from $p(\sigma^2\mid \bA, \bW,\bZ, \alpha_\sigma, \beta_\sigma)$;
Equation (<ref>)
Report loss in Equation (<ref>), stop if it converges.

Prior by Gamma distribution. We also notice that placing an inverse-Gamma prior on the variance parameter of a Gaussian density is equivalent to assigning a Gamma prior on the precision parameter (the inverse of variance).
In representing the precision $\tau=\sigma^{-2}$, we use a Gamma distribution with shape $\alpha_\tau>0$ and rate $\beta_\tau>0$ (Definition <ref>, p. definition:gamma-distribution):

\begin{equation}\label{equation:ggg_gamma_prior}
p(\tau)\sim \gammadist (\tau\mid  \alpha_\tau, \beta_\tau) = \frac{\beta_{\tau}^{\alpha_{\tau}}}{\Gamma(\alpha_{\tau})} \tau^{\alpha_{\tau}-1} \exp({-\beta_{\tau}\cdot \tau}),
\end{equation}

The posterior is obtained similarly (Equation (<ref>), p. equation:gamma_conjugacy_general),

\begin{aligned}
p(\tau \mid  \bW,\bZ, \bA) = \gammadist (\tau; \widetilde{\alpha_\tau}, \widetilde{\beta_\tau}), \\
&\gap \widetilde{\alpha_\tau} = \frac{MN}{2} +{\alpha_\tau}, \gap
\widetilde{\beta_\tau}  =  \frac{1}{2} \sum_{m,n=1}^{M,N} (\bA-\bW\bZ)_{mn}^2 + {\beta_\tau}.
\end{aligned}

In practical applications, the prior parameters $\alpha_\tau$ and $\beta_\tau$ can be chosen to be equal to $\alpha_\sigma$ and $\beta_\sigma$, respectively.

Graphical model representation of GGGA and GGGW models. Green circles denote prior variables, orange circles represent observed and latent variables, and plates represent repeated variables.
The slash “/" in the variable represents “or," and the comma “," in the variable represents “and."

Decomposition: GGGA

§ ALL GAUSSIAN MODEL WITH ARD HIERARCHICAL PRIOR (GGGA)

The all Gaussian with hierarchical Gamma prior (GGGA) model is proposed by Virtanen et al., 2011, Virtanen et al., 2012 based on the GGG model.
The distinction lies in the incorporation of a hyperprior over the Gaussian prior in the GGGA model.
Moreover, the GGGA model  favors an automatic relevance determination (ARD) that helps perform automatic model selection (Figure <ref>).

Automatic relevance determination
Linear model
Linear regression
Bayesian linear regression

Automatic relevance determination (ARD).
Automatic relevance determination stands as a Bayesian technique employed in machine learning and statistical modeling to automatically discern the relevance or significance of input features. ARD is particularly applied in the context of Bayesian linear regression models.

In the Bayesian approach, we model the relationship between input features and the output using a probabilistic framework. In linear regression, this relationship is often represented as:

\bb = \bA\bx +\bepsilon,

where $\bA\in\real^{M\times N}$ is the input data matrix, $\bx\in\real^N$ is a vector of weights, and $\bepsilon\in\real^N$ is a noise term. ARD extends this framework by introducing a probabilistic prior distribution over the weights.
The ARD method assigns a separate precision parameter (inverse variance) to each weight in the model.

The fundamental concept behind ARD is to automatically determine the relevance of input features by allowing the model to adaptively assign large precision (low variance) to irrelevant features and small precision (high variance) to relevant features. Features with small precision effectively become “shrunk" and contribute less to the model, allowing the algorithm to automatically select a subset of features.

In a Bayesian context, the ARD model involves placing a hyperprior over the precision parameters, and the optimization process aims to find the most probable values for these parameters given the observed data. This allows the model to effectively estimate the relevance of each feature while considering the uncertainty associated with the estimation.

The advantages  of ARD include its ability to perform feature selection automatically, thereby reducing overfitting and improving the interpretability of the model. However, it's essential to carefully choose appropriate prior distributions and hyper-parameters to achieve meaningful results.

Going further from the GGG model, we consider the ARD prior, in which case we place a further  Gamma prior (i.e., a hierarchical prior):

\begin{equation}\label{equation:hyper_eq_ggga}
\begin{aligned}
w_{mk} &\sim \normal( 0, (\lambda_{k})^{-1}),
\gap
&z_{kn}\sim&  \normal( 0, (\lambda_{k})^{-1}), \gap
\lambda_{k} \sim \gammadist( \alpha_{\lambda}, \beta_{\lambda}),\\
%\bw_{m} &\sim \normal(\bw_{m}\mid \bzero, \diag(\blambda^{-1})),
%&\bz_{n}\sim&  \normal(\bz_{n}\mid 0, \diag(\blambda^{-1})),
\end{aligned}
\end{equation}

where $k\in\{1,2,\ldots, K\}$, and $\lambda_k$ is shared by all entries in the same column of $\bW$ and the same row of $\bZ$.
The entire factor $k$ is then either activated if $\lambda_k$ has a low value or “turned off" if $\lambda_k$ has a high value (see Figure <ref>, p. fig:dists_gaussian).

Markov blanket

For RMF, following the Bayes' rule and MCMC, this means we need to be able to draw from distributions (again by Markov blanket, Section <ref>):

\begin{aligned}
&p(\sigma^2\mid \bA, \bW, \bZ, \alpha_\sigma, \beta_\sigma ),
&\gap& p(w_{mk}\mid \bA,   \bW_{-mk}, \bZ, \sigma^2, \blambda ), \\
&p(\lambda_k \mid \bW, \bZ, \blambda_{-k}, \sigma^2,  \alpha_\lambda, \beta_\lambda),
&\gap& p(z_{kn}\mid \bA,  \bW, \bZ_{-kn}, \sigma^2, \blambda ), \\
\end{aligned}

where $\blambda\in \real_+^K$ is a vector including all $\lambda_k$ values, and $\blambda_{-k}$ denotes all elements of $\blambda$ except $\lambda_k$.
The conditional posteriors for $\sigma^2$, $w_{mk}$'s, and $z_{kn}$'s remain consistent with the GGG model, except now we replace $\lambda_{mk}^W$ and $\lambda_{kn}^Z$ for all $m\in\{1,2,\ldots, M\}$ and $n\in\{1,2,\ldots, N\}$ by $\lambda_k$.
The conditional posterior density of  $\lambda_{k}$ depends on its parents ($\alpha_\lambda$, $\beta_\lambda$), children ($k$-th row $\bw_k$ of $\bW$, $k$-th column $\bz_k$ of $\bZ$, see definition in Equation (<ref>)), and coparents ($\blambda_{-k}$) [See Figure <ref> and Section <ref>.].
Then it follows that,

\begin{equation}\label{equation:posterior-ggga_lambdak}
\begin{aligned}
&\gap p(\lambda_k\mid  \bW, \bZ, \blambda_{-k}, \sigma^2,  \alpha_\lambda, \beta_\lambda)\\
&\propto  p(\bw_k, \bz_k\mid  \lambda_k) \cdot p(\lambda_k) =
\prod_{i=1}^{M}  \normal(w_{ik}\mid 0, (\lambda_{k})^{-1}) \cdot
\prod_{j=1}^{N} \normal(z_{kj}\mid  0,(\lambda_{k})^{-1})
\cdot \gammadist(\lambda_k \mid  \alpha_{\lambda}, \beta_\lambda)\\
\prod_{i=1}^{M} \lambda_k^{1/2} \exp\left\{-\frac{\lambda_k w_{ik}^2}{2}\right\} \cdot
\prod_{j=1}^{N} \lambda_k^{1/2} \exp\left\{-\frac{\lambda_k z_{kj}^2}{2}\right\}
\cdot \frac{\beta_{\lambda}^{\alpha_\lambda}}{\Gamma(\alpha_\lambda)} \lambda_k^{\alpha_\lambda-1} \exp(- \lambda_k \beta_\lambda)\\
&\propto\lambda_k^{\frac{M+N}{2} +\alpha_\lambda -1}
\exp\left\{-\lambda_k \cdot  \left(\frac{1}{2}\sum_{i=1}^{M}w_{ik}^2 +
\frac{1}{2}\sum_{j=1}^{N}z_{kj}^2 +\beta_\lambda\right)\right\}\\
&\propto \gammadist(\lambda_k\mid  \widetilde{\alpha_\lambda}, \widetilde{\beta_\lambda}),
\end{aligned}
\end{equation}

\widetilde{\alpha_\lambda} = \frac{M+N}{2}+\alpha_\lambda,
\qquad
\widetilde{\beta_\lambda}=\frac{1}{2}\sum_{i=1}^{M}w_{ik}^2 +
\frac{1}{2}\sum_{j=1}^{N}z_{kj}^2 +\beta_\lambda.

According to the definition of the Gamma distribution (Definition <ref>, p. definition:gamma-distribution), we obtain the moments of the posterior density for $\lambda_k$:

\Exp[\lambda_k] = \frac{\widetilde{\alpha_\lambda}}{\widetilde{\beta_\lambda}},
\qquad
\Var[\lambda_k] = \frac{\widetilde{\alpha_\lambda}}{\widetilde{\beta_\lambda}^2}.

Therefore, upon observing the posterior parameters, it is evident that  with a larger  shape of the raw matrix $\bA$ (i.e., $M+N$ is larger), there is a preference for  a larger value of $\lambda_k$ (during sampling, since $\widetilde{\alpha_\lambda}$ tends to be larger).
As indicated in Equation (<ref>), this preference imposes a larger and sparser regularization over the model.
This is reasonable since a larger shape  indicates the vector product for each entry of $\bA$ through $\bW$ and $\bZ$ involves more entries to sum up.

Moreover, if the elements of the factored components $\bW$ and $\bZ$ are larger, the value of $\widetilde{\beta_\lambda}$ tends to be larger as well. Consequently, a smaller value of $\lambda_k$ will be drawn. This is reasonable in the sense that we want to explore in a larger space if the factored components have larger elements from previous Gibbs iterations.

Gibbs sampling.
Again we can construct a Gibbs sampler for the GGGA model, as formulated in Algorithm <ref>. By default, uninformative hyper-parameters are $\alpha_\sigma=\beta_\sigma=1$, $\alpha_\lambda=\beta_\lambda=1$.

Gibbs sampler for GGGA model in one iteration (prior on variance $\sigma^2$ here, similarly for the precision $\tau$).
The procedure presented here may not be efficient but is explanatory. A more efficient one can be implemented in a vectorized manner.
By default, uninformative hyper-parameters are $\alpha_\sigma=\beta_\sigma=1$, $\alpha_\lambda=\beta_\lambda=1$.

[1]
Choose initial $\alpha_\sigma, \beta_\sigma, \alpha_\lambda, \beta_\lambda$;
$k=1$ to $K$
$m=1$ to $M$
Sample $w_{mk}$ from $p(w_{mk}\mid \bA,    \bW_{-mk}, \bZ,\sigma^2,  \lambda_k )$; Equation (<ref>)
$n=1$ to $N$
Sample $z_{kn}$ from $p(z_{kn}\mid \bA ,   \bW, \bZ_{-kn}, \sigma^2, \lambda_k)$;
Symmetry of Eq. (<ref>)

Sample $\lambda_k$ from $p(\lambda_k\mid   \bW, \bZ, \sigma^2, \alpha_\lambda, \beta_\lambda)$;
Equation (<ref>)

Sample $\sigma^2$ from $p(\sigma^2\mid \bA, \bW,\bZ, \alpha_\sigma, \beta_\sigma)$;
Equation (<ref>)
Report loss in Equation (<ref>), stop if it converges.

Decomposition: GGGW

§ ALL GAUSSIAN MODEL WITH WISHART HIERARCHICAL PRIOR (GGGW)

The hierarchical prior with Wishart density is proposed by Salakhutdinov and Mnih, 2008 to increase flexibility and calibration based on the GGG model.
Instead of assuming independence of each entry in the factored components $\bW$ and $\bZ$, we now assume each row $\bw_m$ of $\bW$ and each column $\bz_n$ of $\bZ$ comes from a multivariate Gaussian density (Definition <ref>, p. definition:multivariate_gaussian), whose parameters are placed over a further normal-inverse-Wishart (NIW) prior (Definition <ref>, p. definition:normal_inverse_wishart).
Figure <ref> shows the graphical representation of the GGGW model.

Prior and hyperprior.
Same as the GGG model, we consider the Gaussian likelihood over the data matrix $\bA$, and the variance parameter $\sigma^2$ is placed over an inverse-Gamma prior with shape $\alpha_\sigma$ and scale $\beta_\sigma$.
Given the $m$-th row $\bw_m$ of $\bW$ and the $n$-th column $\bz_n$ of $\bZ$, we consider the multivariate Gaussian density and the normal-inverse-Wishart prior as follows:

\begin{align}
\bw_m&\sim \normal(\bw_m\mid  \bmu_w, \bSigma_w), \gap\gap& \bmu_w, \bSigma_w &\sim \niw(\bmu_w, \bSigma_w\mid  \bmm_0, \kappa_0, \nu_0, \bS_0);\label{equa:gggw_prior_hyper1}\\
\bz_n&\sim \normal(\bz_n\mid  \bmu_z, \bSigma_z), \gap\gap& \bmu_z, \bSigma_z &\sim \niw(\bmu_z, \bSigma_z\mid  \bmm_0, \kappa_0, \nu_0, \bS_0),\label{equa:gggw_prior_hyper2}
\end{align}

where $\niw (\bmu, \bSigma\mid  \bmo, \kappa_0, \nu_0, \bso)
= \mathcal{N}(\bmu\mid  \bmo, \frac{1}{\kappa_0}\bSigma) \cdot  \inversewishart(\bSigma\mid  \bso, \nu_0)$ is the density of a normal-inverse-Wishart distribution, and $ \inversewishart(\bSigma\mid  \bso, \nu_0)$ is the inverse-Wishart distribution (Definition <ref>, p. definition:multi_inverse_wishart).
While  we can also place  normal and inverse-Wishart priors over the mean and covariance parameters separately, i.e., a semi-conjugate prior. We do not repeat the details here, see  Section <ref> and <ref> (p. section:sep_mu_niw and  p. section:sep_mu_niw).

Following from the discussion in Section <ref> (p. sec:niw_posterior_conjugacy), the posterior density of $\{\bmu_w, \bSigma_w\}$ also follows a NIW distribution with updated parameters:

\begin{equation}\label{equation:gggw_post_niw}
\bmu_w, \bSigma_w\sim \niw(\bmu_w, \bSigma_w\mid \bmm_M, \kappa_M, \nu_M, \bS_M),
\end{equation}

\begin{align}
\bmm_M &= \frac{\kappa_0\bmo + M\overline{\bw}}{\kappa_M} =
\frac{\kappa_0 }{\kappa_M}\bmo+\frac{M}{\kappa_M}\overline{\bw} , \label{equation:niw_posterior_equation_2_gggw}\\
\kappa_M  &= \kappa_0 + M,  \label{equation:niw_posterior_equation_3_gggw}\\
\nu_M     &=\nu_0 + M,  \label{equation:niw_posterior_equation_4_gggw}\\
\bS_M
&=\bS_0 + \bS_{\overline{w}} + \frac{\kappa_0 M}{\kappa_0 + M}(\overline{\bw} - \bmo)(\overline{\bw} - \bmo)^\top \label{equation:niw_posterior_equation_5_gggw}\\
&=\bS_0 + \sum_{m=1}^M \bw_m \bw_m^\top + \kappa_0 \bmo \bmo^\top - \kappa_M \bmm_M \bmm_M^\top,  \label{equation:niw_posterior_equation_6_gggw}\\
\overline{\bw} &= \frac{1}{M}\sum_{m=1}^{M} \bw_m, \label{equation:niw_posterior_equation_7_gggw}\\
\bS_{\overline{w}} &= \sum_{m=1}^{M} (\bw_m - \overline{\bw})(\bw_m - \overline{\bw})^\top. \label{equation:niw_posterior_equation_8_gggw}
\end{align}

An intuitive interpretation for the parameters in NIW can be obtained from the updated parameters above.
The parameter $\nu_0$ is the prior number of samples to observe the covariance matrix, and $\nu_M =\nu_0 + M$ is the posterior number of samples.
The posterior mean $\bmm_M$ of the model mean $\bmu_w$ is a weighted average of the prior mean $\bmm_0$ and the sample mean $\overline{\bw}$. The posterior scale matrix $\bS_M$ is the sum of the prior scale matrix $\bS_0$, empirical covariance matrix $\bS_{\overline{w}}$, and an extra term due to the uncertainty in the mean.
And due to symmetry, a similar form for $\{\bmu_z, \bSigma_z\}$ can be derived.

Gibbs sampling.
Once again, we can construct a Gibbs sampler for the GGGW model as formulated in Algorithm <ref>. By default, uninformative hyper-parameters are $\alpha_\sigma=\beta_\sigma=1$, $\bmm_0=\bzero, \kappa_0=1, \nu_0=K+1, \bS_0=\bI$.

Gibbs sampler for GGGW model in one iteration (prior on variance $\sigma^2$ here, similarly for the precision $\tau$).
By default, uninformative hyper-parameters are $\alpha_\sigma=\beta_\sigma=1$, $\bmm_0=\bzero, \kappa_0=1, \nu_0=K+1, \bS_0=\bI$.

[1]
Choose initial $\alpha_\sigma, \beta_\sigma, \bmm_0, \kappa_0, \nu_0, \bS_0$;

$m=1$ to $M$
Sample $\bw_m$ from $p(\bw_m\mid \bmu_w, \bSigma_w )$; Equation (<ref>)
$n=1$ to $N$
Sample $\bz_n$ from $p(\bz_n\mid\bmu_z, \bSigma_z)$;
Symmetry of Eq. (<ref>)

Sample $\bmu_w, \bSigma_w$ from $p(\bmu_w, \bSigma_w\mid \bmm_M, \kappa_M, \nu_M, \bS_M )$
Equation (<ref>)
Sample $\bmu_z, \bSigma_z$ from $p(\bmu_z, \bSigma_z\mid \bmm_N, \kappa_N, \nu_N, \bS_N )$
Symmetry of Eq. (<ref>)

Sample $\sigma^2$ from $p(\sigma^2\mid \bA, \bW,\bZ, \alpha_\sigma, \beta_\sigma)$;
Equation (<ref>)
Report loss in Equation (<ref>), stop if it converges.

Decomposition: GVG

§ GAUSSIAN LIKELIHOOD WITH VOLUME AND GAUSSIAN PRIORS (GVG)

The Gaussian likelihood with volume and Gaussian prior (GVG) is introduced by Arngren et al., 2011.
While the original paper applies the volume prior to unmix a set of pixels into pure spectral signatures (endmembers) and corresponding fractional abundances in hyperspectral image analysis, in which case the factored components are nonnegative.
However, it can also be applied in the real-valued applications.

In the GVG model, the prior over $\bZ$ is still a Gaussian density as that in the GGG model. Instead of assuming Gaussian prior over $\bW$, Arngren et al., 2011 put a volume prior for the factored component $\bW$ with density $\bW\propto \exp\{-\gamma \det(\bW^\top\bW)\}$ (Figure <ref>).
The prior has a single parameter $\gamma$ that is manually determined.

Graphical representation of GVG model. Green circles denote prior variables, orange circles represent observed and latent variables, and plates represent repeated variables.
The slash “/" in the variable represents “or."

Denote vector $\bw_{m,-k}\in \real^{K-1}$ as the $m$-th row of $\bW$ excluding column $k$; vector $\bw_{-m, k}\in\real^{M-1}$ as the $k$-th column of $\bW$ excluding row $m$; matrix $\bW_{-m,-k}\in \real^{(M-1)\times (K-1)}$ as $\bW$ excluding row $m$ and column $k$; matrix $\bW_{:,-k}\in \real^{M\times (K-1)}$ as $\bW$ excluding column $k$;
scalar value $D_{-k,-k}=\det\big( \bW_{:,-k}^\top \bW_{:,-k} \big)$;
and the matrix adjugate of $\big( \bW_{:,-k}^\top \bW_{:,-k} \big)$ as
$\bA_{-k,-k}=\det\big( \bW_{:,-k}^\top \bW_{:,-k} \big)\big( \bW_{:,-k}^\top \bW_{:,-k} \big)^{-1}\in\real^{(K-1)\times(K-1)}$.
Then the posterior density of $w_{mk}$ can be obtained by

\begin{equation}\label{equation:gvg_post_w}
w_{mk} \sim \normal(w_{mk}\mid  \widetilde{\mu_{mk}}, \widetilde{\sigma_{mk}^{2}}),
\end{equation}

\begin{aligned}
\widetilde{\mu_{mk}}&=
\widetilde{\sigma_{mk}^{2}}
\left\{
\gamma \bw_{m,-k}^\top \bA_{-k,-k}(\bW_{-m,-k}^\top)\bw_{-m,k} +
\frac{1}{\sigma^2} \sum_{j=1}^{N} (a_{mj} - \sum_{i\neq k}^{K} \bw_{m,-k}^\top \bz_{j, -k})z_{kj}
\right\},
\end{aligned}

\begin{aligned}
\widetilde{\sigma_{mk}^{2}}&=
1\bigg/ \left(	\frac{1 }{\sigma^2} \sum_{j=1}^{N} z_{kj}^2 +
\gamma\left(D_{-k,-k} - \bw_{m,-k}^\top \bA_{-k,-k}\bw_{m,-k}\right)
\right).
\end{aligned}

The above result can be obtained by extracting $w_{mk}$ from $\exp\{-\gamma \det(\bW^\top\bW)\}$ and from the fact that
$\det(\bM) = \det(\bD)\det(\bA-\bB\bD^{-1}\bC) = \det(\bA)\det(\bD-\bC\bA^{-1}\bB)$
if the  matrix $\bM$ has the block formulation $\bM=
\begin{bmatrix}
\bA & \bB\\
\bC & \bD
\end{bmatrix}$.

* Use the “MovieLens 100K" data set introduced in Section <ref> (p. section:als_movie_rec) to assess and compare the effectiveness of Bayesian real matrix factorization approaches introduced in this Chapter.

*  Following the derivation in Equation (<ref>), derive the conditional distribution over the user feature $z_{kn}$ for all $k\in\{1,2,\ldots, K\}$ and $n\in\{1,2,\ldots, N\}$ for the GGG model.

Similarly, following the derivation in Equation (<ref>), obtain the conditional distribution over the user feature $\bz_{n}$ for all  $n\in\{1,2,\ldots, N\}$ within the context of the GGGM model.

* Derive the automatic relevance determination related results for the linear regression models based on the discussions in Section <ref> and Section <ref> (p. section:bayes_appetizers).

CHAPTER: BAYESIAN NONNEGATIVE MATRIX FACTORIZATION

§ INTRODUCTION

The nonnegative matrix factorization (NMF) method is employed for analyzing data matrices with nonnegative elements, which are common in data sets derived from texts and images [Berry et al., 2007].
In cases where the entries in $\bA, \bW$, and $\bZ$ are nonnegative, NMF algorithms have frequently improved performance.
Thus, the scope of NMF research has grown rapidly in recent years, particularly in the fields of machine learning [Lee and Seung, 1999, Lee and Seung, 2000].

Early work on nonnegative matrix factorizations was conducted by a Finnish group of researchers in the 1990s under the name positive matrix factorization [Paatero et al., 1991, Paatero and Tapper, 1994, Anttila et al., 1995]. This work is rarely cited by later researchers  partly due to the unfortunate phrasing of positive matrix factorization, which is misleading as Paatero and Tapper, 1994 actually create a nonnegative matrix factorization.
Since its introduction by Lee and Seung, 1999, Lee and Seung, 2000, the NMF problem has received a significant amount of attention, both in published and unpublished work, in various fields such as science, engineering, and medicine. Different authors have also proposed alternative formulations for the NMF problem [Schmidt et al., 2009, Tan and Févotte, 2013, Brouwer and Lio, 2017, Lu and Ye, 2022].

The NMF problem can be stated as $\bA=\bW\bZ+\bE$, where $\bA\in \real^{M\times N}$ is approximately factorized into an $M\times K$ nonnegative matrix $\bW\in \real_+^{M\times K}$ and a  $K\times N$ nonnegative matrix $\bZ\in \real_+^{K\times N}$. The data set $\bA$ needs not be complete such that the indices of observed entries can be represented by the mask matrix $\bM\in \real^{M\times N}$, where an entry of one indicates the element is observed, and zero indicates an unobserved element.
This nonnegativity makes the resulting matrices easier to inspect and more intuitive to interpret, especially in the image analysis context.

To simplify the problem, let us assume that there are no missing entries firstly.
Treating   missing entries in the Bayesian NMF context is just the same as that in the Bayesian RMF cases (Section <ref>, p. section:markov-blanket).
Projecting data vectors $\ba_n$ to a smaller
dimension $\bz_n \in \real^K$  with $K<M$,
such that the reconstruction error measured by Frobenius norm is minimized (assume $K$ is known):

\begin{equation}\label{equation:als-per-example-loss_bnmf}
\mathop{\min}_{\bW,\bZ} L(\bW,\bZ) =  \mathop{\min}_{\bW,\bZ}\sum_{n=1}^N \sum_{m=1}^{M} \left(a_{mn} - \bw_m^\top\bz_n\right)^2,
\end{equation}

where $\bW=[\bw_1^\top; \bw_2^\top; \ldots; \bw_M^\top]\in \real^{M\times K}$ and $\bZ=[\bz_1, \bz_2, \ldots, \bz_N] \in \real^{K\times N}$, containing $\bw_m$'s and $\bz_n$'s as rows and columns, respectively.

Continuing the terminology established in Bayesian RMF, the Bayesian NMF models are denoted by their density functions, arranged in the order of likelihood, priors, and hyperpriors  (Section <ref>, p. section:bmf_real_intro).
Table <ref> summarizes the Bayesian models for nonnegative matrix factorization in this chapter.

Name     Likelihood    Prior $\bW$     Prior $\bZ$      Hierarchical prior

GEE     $\normal(a_{mn}|\bw_m^\top\bz_n, \sigma^2)$
$\exponential(w_{mk}|\lambda_{mk}^W)$     $\exponential(z_{kn}|\lambda_{kn}^Z)$
GEEA     $\normal(a_{mn}|\bw_m^\top\bz_n, \sigma^2)$
$\exponential(w_{mk}|\lambda_k)$     $\exponential(z_{kn}|\lambda_k)$     $\gammadist(\lambda_k|\alpha_\lambda, \beta_\lambda)$
GTT    $\normal(a_{mn}|\bw_m^\top\bz_n, \sigma^2)$
$\truncatednormal(w_{mk}|\mu_{mk}^W, \frac{1}{\tau_{mk}^W})$     $\truncatednormal(z_{kn}|\mu_{kn}^Z, \frac{1}{\tau_{kn}^Z})$
GTTN    $\normal(a_{mn}|\bw_m^\top\bz_n, \sigma^2)$
$\truncatednormal(w_{mk}|\mu_{mk}^W, \frac{1}{\tau_{mk}^W})$     $\truncatednormal(z_{kn}|\mu_{kn}^Z, \frac{1}{\tau_{kn}^Z})$
$\{\mu_{mk}^W, \tau_{mk}^W\}$, $\{\mu_{kn}^Z, \tau_{kn}^Z\}\sim$

$\tnsng(\mu_\mu, \tau_\mu, a, b)$
GRR     $\normal(a_{mn}|\bw_m^\top\bz_n, \sigma^2)$
$\rectifieddist(\cdot|\mu_{mk}^W, \frac{1}{\tau_{mk}^W}, \lambda_{mk}^W)$     $\rectifieddist(\cdot|\mu_{kn}^Z, \frac{1}{\tau_{kn}^Z}, \lambda_{kn}^Z)$
GRRN     $\normal(a_{mn}|\bw_m^\top\bz_n, \sigma^2)$
$\rectifieddist(\cdot|\mu_{mk}^W, \frac{1}{\tau_{mk}^W}, \lambda_{mk}^W)$     $\rectifieddist(\cdot|\mu_{kn}^Z, \frac{1}{\tau_{kn}^Z}, \lambda_{kn}^Z)$
$\{\mu_{mk}^W, \tau_{mk}^W, \lambda_{mk}^W\}$,

$\{\mu_{kn}^Z, \tau_{kn}^Z, \lambda_{kn}^Z\}\sim$

$\rnsng(\mu_\mu, \tau_\mu, a, b, \alpha_\lambda, \beta_\lambda)$

$\normal(a_{mn}|\bw_m^\top\bz_n, \sigma^2)$

$\{ \frac{-\lambda_k^W}{2}\sum_{m}(\sum_{k}w_{mk})^2 \}$

$\{ \frac{-\lambda_k^Z}{2}\sum_{n}(\sum_{k}z_{kn})^2 \}$

$\normal(a_{mn}|\bw_m^\top\bz_n, \sigma^2)$

$\{ \frac{-\lambda_k^W}{2}\sum_{m}(\sum_{k}w_{mk}^2) \}$

$\{ \frac{-\lambda_k^Z}{2}\sum_{n}(\sum_{k}z_{kn}^2) \}$

$\normal(a_{mn}|\bw_m^\top\bz_n, \sigma^2)$

$\bW\sim\exp\{ -$

$\lambda_k^W\sum_{m}(\max_k\abs{w_{mk}}) \}$

$\bZ\sim\exp\{ -$

$\lambda_k^Z\sum_{n}(\max_k\abs{z_{kn}}) \}$

GEG     $\normal(a_{mn}|\bw_m^\top\bz_n, \sigma^2)$
$\exponential(w_{mk}|\lambda_{mk}^W)$     $\normal(z_{kn}|0, (\lambda_{kn}^Z)^{-1})$
GnVG     $\normal(a_{mn}|\bw_m^\top\bz_n, \sigma^2)$

$\exp\{-\gamma \bW^\top\bW\}u(\bW)$
$\normal(z_{kn}|0, (\lambda_{kn}^Z)^{-1})$

Overview of Bayesian nonnegative and semi-nonnegative matrix factorization models.

Graphical model representation of GEE and GEEA models. Green circles denote prior variables, orange circles represent observed and latent variables, and plates represent repeated variables.
The slash “/" in the variable represents “or."

Decomposition: GEE

§ GAUSSIAN LIKELIHOOD WITH EXPONENTIAL PRIORS (GEE)

The Gaussian likelihood with exponential priors (GEE) model stands out as one of the simplest Bayesian NMF models, employing a Gaussian likelihood for  data elements and exponential priors  over factor matrices [Schmidt et al., 2009].

Again, we view the data $\bA$ as being produced according to the probabilistic generative process shown in Figure <ref>.
We assume the residuals, $e_{mn}$, are i.i.d. drawn from a  zero-mean Gaussian distribution with variance $\sigma^2$. This is equivalent to assuming the observed $(m,n)$-th data entry $a_{mn}$ of matrix $\bA$ is modeled using a Gaussian likelihood with variance $\sigma^2$ and a mean given by the latent decomposition $\bw_m^\top\bz_n$ (Equation (<ref>)).
Consequently, this gives rise to the following  likelihood function:

\begin{equation}\label{equation:gee_likelihood}
\begin{aligned}
p(\bA\mid \btheta) &= \prod_{m,n=1}^{M,N} \normal \left(a_{mn}\mid (\bW\bZ)_{mn}, \sigma^2 \right)\\
&= \prod_{m,n=1}^{M,N} \normal \left(a_{mn}\mid (\bW\bZ)_{mn}, \tau^{-1} \right),
\end{aligned}
\end{equation}

where $\btheta=\{\bW,\bZ,\sigma^2\}$ denotes all parameters in the model, $\sigma^2$ is the variance, $\tau^{-1}=\sigma^2$ is the precision, and  $\normal(x\mid \mu,\sigma^2)$ is the normal density function.

We treat the latent variables $w_{mk}$'s (and $z_{kn}$'s) as random variables.
And we need prior densities over these latent variables to express beliefs about their values, e.g., nonnegativity in this context, though there are many other possible constraints (semi-nonnegativity in Ding et al., 2008, discreteness in Gopalan et al., 2014, Gopalan et al., 2015).
Here, we assume $w_{mk}$'s and $z_{kn}$'s are independently drawn from   exponential priors with rate parameters $\lambda_{mk}^W$ and $\lambda_{kn}^Z$, respectively (Definition <ref>, p. definition:exponential_distribution):

\begin{equation}\label{equation:gee_prior_density_exponential}
\begin{aligned}
w_{mk} &\sim \exponential(w_{mk}\mid \lambda_{mk}^W),
\gap
&z_{kn}\sim&  \exponential(z_{kn}\mid \lambda_{kn}^Z);\\
p(\bW) &=\prod_{m,k=1}^{M,K} \exponential(w_{mk}\mid \lambda_{mk}^W),
\gap
&p(\bZ) =&\prod_{k,n=1}^{K,N} \exponential(z_{kn}\mid \lambda_{kn}^Z),
\end{aligned}
\end{equation}

where $\exponential(x\mid \lambda)=\lambda\exp(-\lambda x)u(x)$ is the exponential density with $u(x)$ being the unit step function.
This prior serves to enforce the nonnegativity constraint on the components $\bW$ and $\bZ$.
Same as the GGG model, the prior for the noise variance $\sigma^2$ is chosen as an inverse-Gamma density with shape ${\alpha_\sigma}$ and scale ${\beta_\sigma}$ (Definition <ref>, p. definition:inverse_gamma_distribution):

\begin{equation}\label{equation:geg_sigma_prior}
p(\sigma^2)= \inversegammadist(\sigma^2\mid \alpha_\sigma, \beta_\sigma) = \frac{{\beta_\sigma}^{\alpha_\sigma}}{\Gamma({\alpha_\sigma})} (\sigma^2)^{-\alpha_\sigma-1} \exp\left( -\frac{{\beta_\sigma}}{\sigma^2} \right).
\end{equation}

Again by Bayes' rule (Equation (<ref>), p. equation:posterior_abstract_for_mcmc), the posterior is proportional to the product of likelihood and prior, it can be maximized to yield an estimate of $\bW$ and $\bZ$.

Markov blanket

For NMF, following the Bayes' rule and MCMC, this means we need to be able to draw from distributions (by Markov blanket, Section <ref>, p. section:markov-blanket):

\begin{aligned}
&p(w_{mk}\mid \bA,  \bW_{-mk}, \bZ, \sigma^2,\blambda^W), \\
& p(z_{kn}\mid \bA,  \bW, \bZ_{-kn},  \sigma^2,\blambda^Z ), \\
& p(\sigma^2 \mid \bA, \bW, \bZ, \blambda^W, \blambda^Z ), \\
\end{aligned}

where $\blambda^W$ is an $M\times K$ matrix containing all $\{\lambda_{mk}^W\}$ entries, $\blambda^Z$ is a $K\times N$ matrix including all $\{\lambda_{kn}^Z\}$ values, and $\bW_{-{mk}}$ denotes all elements of $\bW$ except $w_{mk}$.
Using Bayes' theorem,
the conditional density of $w_{mk}$ depends on its parents ($\lambda_{mk}^W$), children ($a_{mn}$), and coparents ($\tau$ or $\sigma^2$, $\bW_{-mk}, \bZ$) [See Figure <ref> and Section <ref> (p. section:markov-blanket).]. And it can be obtained by

\begin{equation}\label{equation:gee_poster_wmk1}
\begin{aligned}
&\gap p(w_{mk} \mid \bA ,   \bW_{-mk}, \bZ, \sigma^2, \lambda_{mk}^W ) \\
&\propto p(\bA\mid \bW, \bZ, \sigma^2) \times p(w_{mk}\mid \lambda_{mk}^W)
=\prod_{i,j=1}^{M,N} \normal \left(a_{ij}\mid \bw_i^\top\bz_j, \sigma^2 \right)\times \exponential(w_{mk}\mid \lambda_{mk}^W) \\
&\propto \exp\left\{   -\frac{1}{2\sigma^2}  \sum_{i,j=1}^{M,N}(a_{ij} - \bw_i^\top\bz_j  )^2\right\}  \times \cancel{\lambda_{mk}^W }\exp(-\lambda_{mk}^W \cdot w_{mk})u(w_{mk})\\
&\propto \exp\left\{   -\frac{1}{2\sigma^2}  \sum_{j=1}^{N}(a_{mj} - \bw_m^\top\bz_j  )^2\right\}  \cdot  \exp(-\lambda_{mk}^W\cdot w_{mk})u(w_{mk})\\
&\propto \exp\left\{   -\frac{1}{2\sigma^2}  \sum_{j=1}^{N}
\left( w_{mk}^2z_{kj}^2 + 2w_{mk} z_{kj}\bigg(\sum_{i\neq k}^{K}w_{mi}z_{ij} - a_{mj}\bigg)  \right)
\right\}  \cdot \exp(-\lambda_{mk}^W\cdot w_{mk})u(w_{mk})\\
&\propto \exp\left\{
-\underbrace{\left(\frac{\sum_{j=1}^{N} z_{kj}^2 }{2\sigma^2}  \right) }_{\textcolor{mylightbluetext}{1/(2\widetilde{\sigma_{mk}^{2}})}}
w_{mk}\underbrace{\left( -\lambda_{mk}^W+ \frac{1}{\sigma^2} \sum_{j=1}^{N} z_{kj}\bigg( a_{mj} - \sum_{i\neq k}^{K}w_{mi}z_{ij}\bigg)  \right)}_{\textcolor{mylightbluetext}{\widetilde{\sigma_{mk}^{2}}^{-1} \widetilde{\mu_{mk}}}}
\right\}  \cdot u(w_{mk})\\
&\propto   \normal(w_{mk} \mid \widetilde{\mu_{mk}}, \widetilde{\sigma_{mk}^{2}})\cdot u(w_{mk})
= \truncatednormal(w_{mk} \mid \widetilde{\mu_{mk}}, \widetilde{\sigma_{mk}^{2}}),
\end{aligned}
\end{equation}

where $u(x)$ is the unit function with value 1 if $x\geq 0$ and value 0 if $x<0$,

\begin{equation}\label{equation:gee_posterior_variance}
\widetilde{\sigma_{mk}^{2}}= \frac{\sigma^2}{\sum_{j=1}^{N} z_{kj}^2}
\end{equation}

is the posterior “parent" variance of the normal distribution with “parent" mean $\widetilde{\mu_{mk}}$,

\begin{equation}\label{equation:gee_posterior_mean}
\widetilde{\mu_{mk}} = \left( -\lambda_{mk}^W+ \frac{1}{\sigma^2} \sum_{j=1}^{N} z_{kj}\bigg( a_{mj} - \sum_{i\neq k}^{K}w_{mi}z_{ij}\bigg)  \right)\cdot \widetilde{\sigma_{mk}^{2}}
\end{equation}

and $\truncatednormal(x \mid \mu, \sigma^2)$ is the truncated-normal (TN) density with “parent" mean $\mu$ and “parent" variance $\sigma^2$ (Definition <ref>, p. definition:truncated_normal).

$\ell_1$ norm
$\ell_1$ constraint
[hidealllines=,backgroundcolor=,frametitle=Interpretation of the Posterior: Sparsity Constraint]
The exponential prior serves as an equivalent to imposing an $\ell_1$ norm, introducing a sparsity constraint in the GEE model.
The sparsity arises from the negative term $-\lambda_{mk}^W$ in Equation (<ref>).
When $\lambda_{mk}^W$ becomes larger, the posterior “parent" mean becomes smaller, and the TN distribution will have a  larger probability for smaller values since the draws of $\truncatednormal(w_{mk} \mid \widetilde{\mu_{mk}}, \widetilde{\sigma_{mk}^{2}})$ will be around zero, thus imposing sparsity (see Figure <ref>, p. fig:dists_truncatednorml_mean).

Or after rearrangement, the posterior density of $w_{mk}$ can be equivalently described by

\small
\begin{aligned}
&\gap p(w_{mk}\mid \bA,  \bW_{-mk}, \bZ, \sigma^2, \blambda^W )\\
&\propto \exp\left\{   (\frac{-1 }{2\sigma^2}  \sum_{j=1}^{N} z_{kj}^2)  w_{mk}^2 +
w_{mk}\underbrace{\left( \frac{1}{\sigma^2} \sum_{j=1}^{N} z_{kj}\bigg( a_{mj} - \sum_{i\neq k}^{K}w_{mi}z_{ij}\bigg)  \right)}_{\textcolor{mylightbluetext}{\widehat{\sigma_{mk}^{2}}^{-1} \widehat{\mu_{mk}}}}
\right\}   \exp(-\lambda_{mk}^W w_{mk}) u(w_{mk})\\
&\propto \normal(w_{mk} \mid \widehat{\mu_{mk}}, \widehat{\sigma_{mk}^{2}})\cdot \exponential(w_{mk}\mid\lambda_{mk}^W) = \rectifieddist(w_{mk} \mid   \widehat{\mu_{mk}}, \widehat{\sigma_{mk}^{2}}, \lambda_{mk}^W ),
\end{aligned}

where $\widehat{\sigma^2_{mk}}=\widetilde{\sigma_{mk}^{2}}= \frac{\sigma^2}{\sum_{j=1}^{N} z_{kj}^2}$ is the posterior “parent" variance of the normal distribution with “parent" mean $\widehat{\mu_{mk}}$,

\widehat{\mu_{mk}} =  \frac{1}{\sum_{j=1}^{N} z_{kj}^2} \cdot \sum_{j=1}^{N} z_{kj}\bigg( a_{mj} - \sum_{i\neq k}^{K}w_{mi}z_{ij}\bigg),

and $\rectifieddist(x\mid \mu, \sigma^2, \lambda)\propto \normal(x\mid \mu, \sigma^2)
\exponential(x\mid \lambda)$ is the rectified-normal (RN) density (Definition <ref>, p. definition:reftified_normal_distribution). And due to symmetry, a similar expression for $z_{kn}$ can be easily derived. The conditional density of $\sigma^2$ depends on its parents ($\alpha_\sigma$, $\beta_\sigma$), children ($\bA$), and coparents ($\bW$, $\bZ$). And it is an inverse-Gamma distribution (by conjugacy in Equation (<ref>), p. equation:inverse_gamma_conjugacy_general):

\begin{equation}\label{equation:gee_posterior_sigma2}
\begin{aligned}
p(\sigma^2 \mid \bA, \bW,\bZ, \alpha_\sigma, \beta_\sigma) = \inversegammadist (\sigma^2\mid \widetilde{\alpha_{\sigma}}, \widetilde{\beta_{\sigma}}), \qquad \\
& \widetilde{\alpha_{\sigma}} = \frac{MN}{2} +{\alpha_\sigma}, \qquad
\widetilde{\beta_{\sigma}}  =  \frac{1}{2} \sum_{m,n=1}^{M,N} (\bA-\bW\bZ)_{mn}^2 + {\beta_\sigma}.
\end{aligned}
\end{equation}

Gibbs sampler for GEE model in one iteration (prior on variance $\sigma^2$ here, similarly for the precision $\tau$). The procedure presented here may not be efficient but is explanatory. A more  efficient one can be implemented in a vectorized manner. By default,
uninformative hyper-parameters are $\alpha_\sigma=\beta_\sigma=1$, $\{\lambda_{mk}^W\}= \{\lambda_{kn}^Z\}=0.1$.

[1]
Choose initial $\alpha_\sigma, \beta_\sigma, \lambda_{mk}^W, \lambda_{kn}^Z$;
$k=1$ to $K$
$m=1$ to $M$
Sample $w_{mk}$ from $p(w_{mk} \mid \bA ,   \bW_{-mk}, \bZ, \sigma^2, \lambda_{mk}^W)$;
Equation (<ref>)
$n=1$ to $N$
Sample $z_{kn}$ from $p(z_{kn} \mid \bA ,   \bW, \bZ_{-kn},\sigma^2 \lambda_{kn}^Z )$;
Symmetry of Eq. (<ref>)
Sample $\sigma^2$ from $p(\sigma^2 \mid  \bA, \bW,\bZ,\alpha_\sigma,\beta_\sigma)$;
Equation (<ref>)
Report loss in Equation (<ref>), stop if it converges.

Inverse-Gamma distribution
Rectified-normal distribution

Gibbs sampling.
By this Gibbs sampling method introduced in Section <ref> (p. section:gibbs-sampler), we can construct a Gibbs sampler for the GEE model as formulated in Algorithm <ref>. And also in practice, all the parameters of the exponential distribution are set to be the same value $\lambda=\{\lambda_{mk}^W\}'s = \{\lambda_{nk}^Z\}'s$ for all $m,k,n$.
By default,
uninformative hyper-parameters are $\alpha_\sigma=\beta_\sigma=1$, $\{\lambda_{mk}^W\}= \{\lambda_{kn}^Z\}=0.1$.

Decomposition: GEEA
Automatic relevance determination
Markov blanket

§ GAUSSIAN LIKELIHOOD WITH EXPONENTIAL PRIORS AND ARD HIERARCHICAL PRIOR (GEEA)

The Gaussian likelihood with exponential priors and hierarchical prior  (GEEA) model is proposed by Tan and Févotte, 2013 based on the GEE model.
The difference lies in that the GEEA model applies a hyperprior on the exponential prior. Moreover, GEEA favors an automatic relevance determination (ARD) that helps perform automatic model selection.
This ARD works by replacing the individual scale parameter of exponential prior for the factored components $\bW$ and $\bZ$ with one that is shared by all entries in the same column of $\bW$ and the same row of $\bZ$. In other words, the parameters are shared for each factor.

Gamma distribution

For each prior density in Equation (<ref>), we assume a Gamma distribution on the parameters of exponential distributions in Equation (<ref>):

w_{mk}\sim \exponential(w_{mk}\mid \lambda_{k}), \gap
z_{kn}\sim \exponential(z_{kn}\mid \lambda_{k}), \gap
\lambda_{k} \sim \gammadist(\lambda_{k} \mid \alpha_\lambda, \beta_\lambda),

where $\lambda_k$ is shared by all entries in the same column of $\bW$ and the same row of $\bZ$. The entire factor $k$ is then either activated if $\lambda_k$ has a low value or “turned off" if $\lambda_k$ has a high value (see Figure <ref>, p. fig:dists_exponential).
Therefore, we can establish an upper bound on the number of hidden factors, $K$, instead of choosing the precise value of $K$. The graphical model is shown in Figure <ref>.

For NMF, following  Bayes' rule and MCMC, this means we need to be able to draw from distributions (again by Markov blanket, Section <ref>, p. section:markov-blanket):

\begin{aligned}
&p(\sigma^2 \mid \bA, \bW, \bZ, \blambda  ), &\gap&
p(w_{mk}\mid \bA,  \bW_{-mk}, \bZ,\sigma^2,  \blambda ), \\
&p(\lambda_k \mid \bW, \bZ, \blambda_{-k},    \alpha_\lambda, \beta_\lambda),
&\gap& p(z_{kn}\mid\bA , \bW,\bZ_{-kn},  \sigma^2, \blambda ), \\
\end{aligned}

where $\blambda\in \real_+^K$ is a vector including all $\lambda_k$ values, and $\blambda_{-k}$ denotes all elements of $\blambda$ except $\lambda_k$. The posteriors for $w_{mk}$'s and $z_{kn}$'s are the same as those in the GEE model, except now we replace $\lambda_{mk}^W$ and $\lambda_{kn}^Z$ by $\lambda_k$. The posteriors for $\lambda_k$ can be obtained using Bayes' theorem.
The conditional density of $\lambda_k$ depends on its parents ($\alpha_\lambda, \beta_\lambda$), children ($k$-th column $\widehat{\bw}_k$ or $\bW$, $k$-th row $\widehat{\bz}_k$ of $\bZ$, note we define $\bw_m$ as the $m$-th row of $\bW$ and $\bz_n$ as the $n$-th column of $\bZ$ in Equation (<ref>)), and coparents (none) [See Figure <ref> and Section <ref> (p. section:markov-blanket).]. Then, it follows that,

\begin{equation}\label{equation:posterior-geea_lambdak}
\begin{aligned}
&\gap p(\lambda_k \mid  \bW, \bZ, \alpha_\lambda, \beta_\lambda)\\
&\propto  p(\widehat{\bw}_k, \widehat{\bz}_k \mid \lambda_k) \times p(\lambda_k) =
\prod_{i=1}^{M}  \exponential(w_{ik}\mid \lambda_{k}) \cdot
\prod_{j=1}^{N} \exponential(z_{kj}\mid \lambda_{k})
\times \gammadist(\lambda_k \mid \alpha_{\lambda}, \beta_\lambda)\\
\prod_{i=1}^{M}  \lambda_k \exp(-\lambda_k w_{ik}) \cdot
\prod_{j=1}^{N} \lambda_k \exp(-\lambda_k z_{kj})
\times \frac{\beta_{\lambda}^{\alpha_\lambda}}{\Gamma(\alpha_\lambda)} \lambda_k^{\alpha_\lambda-1} \exp(- \lambda_k \beta_\lambda)\\
&\propto\lambda_k^{M+N +\alpha_\lambda -1} \exp\left\{-\lambda_k \cdot  \left(\sum_{k=1}^{K}(w_{mk}+ z_{kn}) +\beta_\lambda\right)\right\}\\
&\propto \gammadist(\lambda_k \mid \widetilde{\alpha_\lambda}, \widetilde{\beta_\lambda}),
\end{aligned}
\end{equation}

\widetilde{\alpha_\lambda} = M+N+\alpha_\lambda,
\qquad
\widetilde{\beta_\lambda}=\sum_{k=1}^{K}(w_{mk}+ z_{kn}) +\beta_\lambda.

From this posterior form, the prior parameter $\alpha_\lambda$ can be interpreted as the number of prior observations, and $\beta_\lambda$ as the sum of the prior observations. Therefore, week prior parameter can be chosen as $\alpha_\lambda=\beta_\lambda=1$.

Gibbs sampling.
Again we can construct a Gibbs sampler for the GEEA model as formulated in Algorithm <ref>.
By default,
uninformative hyper-parameters are $\alpha_\sigma=\beta_\sigma=1$, $\alpha_\lambda=\beta_\lambda=1$.

Gibbs sampler for GEEA model in one iteration (prior on variance $\sigma^2$ here, similarly for the precision $\tau$).
The procedure presented here may not be efficient but is explanatory. A more efficient one can be implemented in a vectorized manner.
By default,
uninformative hyper-parameters are $\alpha_\sigma=\beta_\sigma=1$, $\alpha_\lambda=\beta_\lambda=1$.

[1]
Choose initial $\alpha_\sigma, \beta_\sigma, \alpha_\lambda, \beta_\lambda$;
$k=1$ to $K$
$m=1$ to $M$
Sample $w_{mk}$ from $p(w_{mk} \mid \bA,   \bW_{-mk}, \bZ, \sigma^2, \lambda_{k} )$; Equation (<ref>)
$n=1$ to $N$
Sample $z_{kn}$ from $p(z_{kn} \mid \bA,   \bW, \bZ_{-kn},\sigma^2, \lambda_{k} )$;
Symmetry of Eq. (<ref>)

Sample $\lambda_k$ from $p(\lambda_k \mid  \bW, \bZ, \alpha_\lambda, \beta_\lambda)$;
Equation (<ref>)

Sample $\sigma^2$ from $p(\sigma^2 \mid  \bA, \bW,\bZ,\alpha_\sigma,\beta_\sigma)$;
Equation (<ref>)
Report loss in Equation (<ref>), stop if it converges.

Decomposition: GTT
Truncated-normal distribution

§ GAUSSIAN LIKELIHOOD WITH TRUNCATED-NORMAL PRIORS (GTT)

The Gaussian likelihood with truncated-normal priors (GTT) model is discussed in Brouwer and Lio, 2017, where truncated-normal (TN) priors
are used over factored matrices (Figure <ref>).
The truncated-normal distribution, a variant of the normal distribution, excludes values smaller than zero (Definition <ref>, p. definition:truncated_normal), allowing it to impose nonnegativity in Bayesian models.
The chosen likelihood function is identical to that of the GEE model (Equation (<ref>)).

We assume $\bW$ and $\bZ$ are independently truncated-normal distributed with mean and precision given by $\{\bmu^W ,\btau^W\}$ and  $\{\bmu^Z,\btau^Z\}$, respectively:

\begin{equation}\label{equation:gtt_prior_wmk}
w_{mk} \sim \truncatednormal(w_{mk} \mid \mu_{mk}^{W}, (\tau_{mk}^W)^{-1} ), \gap
z_{kn} \sim \truncatednormal(z_{kn} \mid \mu_{kn}^Z, (\tau_{kn}^Z)^{-1} ),
\end{equation}

where $\bmu^W$ is an $M\times K$ matrix containing all $\{\mu_{mk}^W\}$ entries, $\bmu^Z$ is a $K\times N$ matrix including all $\{\mu_{kn}^Z\}$ values, $\btau^W$ is an $M\times K$ matrix containing all $\{\tau_{mk}^W\}$ entries, and $\btau^Z$ is a $K\times N$ matrix including all $\{\tau_{kn}^Z\}$ values.

Graphical model representation of GTT and GTTN models. Green circles denote prior variables, orange circles represent observed and latent variables, and plates represent repeated variables.
The slash “/" in the variable represents “or," and the comma “," in the variable represents “and."

Markov blanket

Again, following  Bayes' rule and MCMC, this means we need to be able to draw from distributions (by Markov blanket, Section <ref>, p. section:markov-blanket):

\begin{aligned}
&p(w_{mk}\mid\bA, \bW_{-mk}, \bZ, \sigma^2, \mu_{mk}^W,\tau_{mk}^W), \\
& p(z_{kn}\mid\bA, \bW, \bZ_{-kn}, \sigma^2,  \mu_{kn}^Z, \tau_{kn}^Z), \\
& p(\sigma^2 \mid \bA, \bW, \bZ, \alpha_\sigma, \beta_\sigma ), \\
\end{aligned}

where $\bW_{-{mk}}$ denotes all elements of $\bW$ except $w_{mk}$, and $\bZ_{-kn}$ denotes all elements of $\bZ$ except $z_{kn}$.
Using Bayes' theorem,
the conditional density of $w_{mk}$ depends on its parents ($\mu_{mk}^W$, $\tau_{mk}^W$), children ($a_{mn}$), and coparents ($\tau$ or $\sigma^2$, $\bW_{-mk}, \bZ$) [See Figure <ref> and Section <ref> (p. section:markov-blanket).]. And it can be obtained by (similar to computing the conditional density of $w_{mk}$ in the GEE model, Equation (<ref>))

\begin{equation}\label{equation:gtt_posterior_wmk1}
\small
\begin{aligned}
&\gap p(w_{mk} \mid \sigma^2,   \bW_{-mk}, \bZ, \mu_{mk}^W, \tau_{mk}^W, \bA)
\propto p(\bA\mid \bW, \bZ, \sigma^2) \cdot p(w_{mk}\mid  \mu_{mk}^W, (\tau_{mk}^W)^{-1})\\
&=\prod_{i,j=1}^{M,N} \normal \left(a_{ij}\mid \bw_i^\top\bz_j, \sigma^2 \right)\times
\truncatednormal(w_{mk} \mid \mu_{mk}^W, (\tau_{mk}^W)^{-1} ) \\
%&\propto \exp\left\{   -\frac{1}{2\sigma^2}  \sum_{i,j=1}^{M,N}(a_{ij} - \bw_i^\top\bz_j  )^2\right\}  \times
%\exp\left(- \frac{\tau_{mk}^W}{2} (w_{mk}- \mu_{mk}^W)^2 \right)u(w_{mk})\\
%&\propto \exp\left\{   -\frac{1}{2\sigma^2}  \sum_{j=1}^{N}(a_{mj} - \bw_m^\top\bz_j  )^2\right\}  \cdot  \exp(-\lambda_{mk}^W\cdot w_{mk})u(w_{mk})\\
%&\propto \exp\left\{   -\frac{1}{2\sigma^2}  \sum_{j=1}^{N}
%\left( w_{mk}^2z_{kj}^2 + 2w_{mk} z_{kj}\bigg(\sum_{i\neq k}^{K}w_{mi}z_{ij} - a_{mj}\bigg)  \right)
%\right\}  \cdot \exp(-\lambda_{mk}^W\cdot w_{mk})u(w_{mk})\\
&\propto \exp\bigg\{ - (\frac{ \sum_{j=1}^{N} z_{kj}^2 }{2\sigma^2} +\textcolor{black}{\frac{\tau_{mk}^W}{2}})  w_{mk}^2 +
\big\{ \frac{1}{\sigma^2} \sum_{j=1}^{N} z_{kj}( a_{mj}- \sum_{i\neq k}^{K}w_{mi}z_{ij})  + \textcolor{black}{\tau_{mk}^W \mu_{mk}^W} \big\}
%_{\textcolor{black}{\widetilde{\sigma_{mk}^{2}}^{-1} \widetilde{\mu_{mk}}}}
\bigg\}  u(w_{mk})\\
&\propto   \normal(w_{mk} \mid \widetilde{\mu_{mk}}, \widetilde{\sigma_{mk}^{2}})\cdot u(w_{mk})
= \truncatednormal(w_{mk} \mid \widetilde{\mu_{mk}}, \widetilde{\sigma_{mk}^{2}}),
\end{aligned}
\end{equation}

where $\widetilde{\sigma_{mk}^{2}}= \frac{\sigma^2}{\sum_{j=1}^{N} z_{kj}^2 + \tau_{mk}^W\cdot \sigma^2}$ is the posterior “parent" variance of the normal distribution with “parent" mean $\widetilde{\mu_{mk}}$,

\widetilde{\mu_{mk}} = \left\{ \frac{1}{\sigma^2} \sum_{j=1}^{N} z_{kj}\bigg( a_{mj} - \sum_{i\neq k}^{K}w_{mi}z_{ij}  \bigg)  + \textcolor{black}{\tau_{mk}^W \mu_{mk}^W}\right\}\cdot \widetilde{\sigma_{mk}^{2}}.

Again, due to symmetry, a similar expression for $z_{kn}$ can be immediately derived. Finally, the conditional density of $\sigma^2$ is the same as that in the  GEE model (Equation (<ref>)).

Gibbs sampler for GTT model in one iteration (prior on variance $\sigma^2$ here, similarly for the precision $\tau$).
The procedure presented here may not be efficient but is explanatory. A more   efficient one can be implemented in a vectorized manner.
By default,
uninformative hyper-parameters are $\alpha_\sigma=\beta_\sigma=1$, $\{\mu_{mk}^W\}=\{\mu_{kn}^Z\}=0$, $\{\tau_{mk}^W\}=\{\tau_{kn}^Z\}=0.1$.

[1]
Choose initial $\alpha_\sigma, \beta_\sigma, \mu_{mk}^W, \tau_{mk}^W,\mu_{kn}^Z, \tau_{kn}^Z$;
$k=1$ to $K$
$m=1$ to $M$
Sample $w_{mk}$ from $p(w_{mk} \mid \bA, \bW_{-mk},\bZ,\sigma^2,\mu_{mk}^W,\tau_{mk}^W)$; Equation (<ref>)
$n=1$ to $N$
Sample $z_{kn}$ from $p(z_{kn} \mid \bA, \bW,\bZ_{-kn},\sigma^2,\mu_{kn}^Z,\tau_{kn}^Z)$;
Symmetry of Eq. (<ref>)
Sample $\sigma^2$ from $p(\sigma^2 \mid  \bA, \bW,\bZ,\alpha_\sigma,\beta_\sigma)$;  Equation (<ref>)
Report loss in Equation (<ref>), stop if it converges.

Gibbs sampling.
We can again construct a Gibbs sampler for the GTT model as formulated in Algorithm <ref>. And also in practice, all the parameters of the truncated-normal priors are set to be the same value $\mu^W=\{\mu_{mk}^W\}'s, \mu^Z=\{\mu_{nk}^Z\}'s$, $\tau^W=\{\tau_{mk}^W\}'s,  \tau^Z=\{\tau_{nk}^Z\}'s$ for all $m,k,n$.
By default,
uninformative hyper-parameters are $\alpha_\sigma=\beta_\sigma=1$, $\{\mu_{mk}^W\}=\{\mu_{kn}^Z\}=0$, $\{\tau_{mk}^W\}=\{\tau_{kn}^Z\}=0.1$.

Decomposition: GTTN
Hierarchical prior
Rectified-normal distribution
TN-scaled-normal-Gamma (TNSNG) prior

§ GAUSSIAN LIKELIHOOD WITH TRUNCATED-NORMAL AND HIERARCHICAL PRIORS (GTTN)

This hierarchical prior is proposed in Schmidt and Mohamed, 2009 over a rectified-normal distribution originally, and further discussed in Brouwer and Lio, 2017 based on the GTT model, where the difference lies in that the GTTN model puts a hyperprior on the two parameters of the truncated-normal distribution (Figure <ref>).

We have shown in Equation (<ref>) (p. equation:conjugate_truncated_nonnegative_mean) that the truncated-normal density is a conjugate prior over the nonnegative mean parameter of a Gaussian distribution that is favored in the GTT model. Moreover, if the prior over $\{w_{mk}\}$'s and $\{z_{kn}\}$'s had been a Gaussian, appropriate conjugate priors for the mean and variance would be a normal-inverse-Gamma or a normal-inverse-Chi-square distribution (Equation (<ref>), p. equation:conjugate_nigamma_general; Equation (<ref>), p. equation:nix-posterior). However, these priors are not conjugate to the truncated-normal density. And instead, a convenient prior called TN-scaled-normal-Gamma (TNSNG) distribution is used (or a TN-scaled-normal-inverse-Gamma prior for “parent" mean and “parent" variance parameters) [The original hyperprior in  Schmidt and Mohamed, 2009 is a rectified-normal (RN) scaled one. Here we scale it in the context of truncated-normal density.]:

\begin{aligned}
\mu_{mk}^W, \tau_{mk}^W \mid \mu_\mu, \tau_\mu, a,  b
&\sim \tnsng(\mu_{mk}^W, \tau_{mk}^W \mid \mu_\mu, \tau_\mu, a,  b)\\
&\propto  \frac{1}{\sqrt{\tau_{mk}^W}} \left(1 - \Phi\big(-\mu_{mk}^W\sqrt{\tau_{mk}^W} \big) \right) \cdot
\normal(\mu_{mk}^W\mid \mu_\mu, (\tau_\mu)^{-1})\cdot
\gammadist(\tau_{mk}^W \mid a, b);\\
\mu_{kn}^Z, \tau_{kn}^Z \mid \mu_\mu, \tau_\mu, a,  b
&\sim \tnsng(\mu_{kn}^Z, \tau_{kn}^Z \mid \mu_\mu, \tau_\mu, a,  b).
\end{aligned}

Here we use the same hyper-parameters $\{ \mu_\mu, \tau_\mu, a,  b\}$ over different entries  $\{u_{mk}^W, \tau_{mk}^W \}$ and $\{u_{kn}^Z, \tau_{kn}^Z\}$. However, in rare cases, one may favor different behaviors over $\bW$ and $\bZ$, e.g., small values in $\bW$ and large values in $\bZ$; the hyper-parameters can be chosen as different values accordingly (see the comparison in Figure <ref>). Note that the scaled-normal-Gamma distribution is not simply a product of a normal and a Gamma distribution. It is not easy to sample from this distribution. However, we will see the posteriors have simple forms due to this prior; that's the reason why we add scaled terms in the prior density. The prior can decouple parameters $\mu_{mk}^W$ and $\tau_{mk}^W$, and the posterior conditional densities of them are normal and Gamma respectively due to this convenient scale.

[GTTN with same hyper-parameters. Same as Figure <ref>.]

[GTTN with different hyper-parameters.]

Graphical model representation of GTTN with same and different hyperparameters. Green circles denote prior variables, orange circles represent observed and latent variables, and plates represent repeated variables. The slash “/" in the variable represents “or," and the comma “," in the variable represents “and."

The posteriors for $\{w_{mk}\}$'s, $\{z_{kn}\}$'s, and $\sigma^2$ are the same as those in the GTT model. The posteriors for $\{u_{mk}^W, \tau_{mk}^W\}$ can be obtained using Bayes' rule, where the conditional density of $\{\mu_{mk}^W, \tau_{mk}^W\}$ depend on their parents ($\mu_\mu, \tau_\mu, a,  b$), children ($w_{mk}$), and coparents (none). Then it follows from the likelihood in Equation (<ref>) that the conditional densities of $\mu_{mk}^W$ is

\begin{equation}\label{equation:posterior_gttn_mu_tau1}
\begin{aligned}
&\gap p(\mu_{mk}^W \mid\textcolor{black}{ \tau_{mk}^W}, w_{mk}, \mu_\mu, \tau_\mu, a,  b)\\
&\propto \truncatednormal(w_{mk} \mid \mu_{mk}^W, (\tau_{mk}^W)^{-1} )
\cdot  \frac{1}{\sqrt{\tau_{mk}^W}} \left(1 - \Phi\big(-\mu_{mk}^W\sqrt{\tau_{mk}^W} \big) \right)
\normal(\mu_{mk}^W\mid \mu_\mu, \tau_\mu)
\gammadist(\tau_{mk}^W \mid a, b)\\
&\propto \exp\left\{ -
\underbrace{\frac{\tau_{mk}^W + \tau_\mu}{2}}_{\textcolor{mylightbluetext}{\widetilde{t}/2}}
(\mu_{mk}^W)^2 + \mu_{mk}^W
\underbrace{(\tau_{mk}^W w_{mk} +\tau_{\mu}\mu_\mu)}_{\textcolor{mylightbluetext}{\widetilde{m}\cdot  \widetilde{t}}} \right\}
\propto \normal(\mu_{mk}^W \mid \widetilde{m}, \widetilde{t}^{-1}),
\end{aligned}
\end{equation}

\widetilde{t}=\tau_{mk}^W + \tau_\mu,
\widetilde{m}=(\tau_{mk}^W w_{mk} +\tau_{\mu}\mu_\mu)/\widetilde{t}.
And the conditional density of $\tau_{mk}^W$ is

\begin{equation}\label{equation:posterior_gttn_mu_tau2}
\begin{aligned}
&\gap p(\tau_{mk}^W \mid\textcolor{black}{\mu_{mk}^W }, w_{mk}, \mu_\mu, \tau_\mu, a,  b)\\
&\propto \truncatednormal(w_{mk} \mid \mu_{mk}^W, (\tau_{mk}^W)^{-1} )
\cdot  \frac{1}{\sqrt{\tau_{mk}^W}} \left(1 - \Phi\big(-\mu_{mk}^W\sqrt{\tau_{mk}^W} \big) \right)
\normal(\mu_{mk}^W\mid \mu_\mu, \tau_\mu)
\gammadist(\tau_{mk}^W \mid a, b)\\
&\propto (\tau_{mk}^W)^{a-1} \exp\left\{-\left( b+ \frac{(w_{mk}-\mu_{mk}^W)^2}{2} \right) \tau_{mk}^W \right\}
\propto \gammadist(\tau_{mk}^W \mid \widetilde{a}, \widetilde{b}),
\end{aligned}
\end{equation}

where $\widetilde{a} = a, \widetilde{b}= b+ \frac{(w_{mk}-\mu_{mk}^W)^2}{2}$.
And again due to symmetry, the expressions for $\mu_{kn}^Z$ and $\tau_{kn}^Z$ can be derived accordingly. The Gibbs sampler for the GTTN model is then formulated in Algorithm <ref>. By default,
uninformative hyper-parameters are $\alpha_\sigma=\beta_\sigma=1$, $\mu_\mu=0, \tau_\mu=0.1$,  $a=b=1$.

GibbssSampler for GTTN model in one iteration (prior on variance $\sigma^2$ here, similarly for the precision $\tau$).
The procedure presented here may not be efficient but is explanatory. A more efficient one can be implemented in a vectorized manner.
By default,
uninformative hyper-parameters are $\alpha_\sigma=\beta_\sigma=1$, $\mu_\mu=0, \tau_\mu=0.1$,  $a=b=1$.

[1]
Choose initial $\alpha_\sigma, \beta_\sigma, \mu_\mu, \tau_\mu, a,  b$;
$k=1$ to $K$
$m=1$ to $M$
Sample $w_{mk}$ from $p(w_{mk} \mid\bA ,   \bW_{-mk}, \bZ, \sigma^2, \mu_{mk}^W, \tau_{mk}^W)$; Equation (<ref>)
Sample $\mu_{mk}^W$ from $p(\mu_{mk}^W \mid \tau_{mk}^W, w_{mk}, \mu_\mu, \tau_\mu, a,  b)$;
Equation (<ref>)
Sample $\tau_{mk}^W$ from $p(\tau_{mk}^W \mid\mu_{mk}^W, w_{mk} , \mu_\mu, \tau_\mu, a,  b)$;
Equation (<ref>)
$n=1$ to $N$
Sample $z_{kn}$ from $p(z_{kn} \mid \bA,   \bW, \bZ_{-kn}, \sigma^2,\mu_{kn}^Z, \tau_{kn}^Z )$;
Symmetry of Eq. (<ref>)
Sample $\mu_{kn}^Z$ from $p(\mu_{kn}^Z \mid \tau_{kn}^Z, z_{kn}, \mu_\mu, \tau_\mu, a,  b)$;
Symmetry of Eq. (<ref>)
Sample $\tau_{kn}^Z$ from $p(\tau_{kn}^Z \mid\mu_{kn}^Z, z_{kn} , \mu_\mu, \tau_\mu, a,  b)$;
Symmetry of Eq. (<ref>)

Sample $\sigma^2$ from $p(\sigma^2 \mid  \bA, \bW,\bZ,\alpha_\sigma,\beta_\sigma)$;  Equation (<ref>)
Report loss in Equation (<ref>), stop if it converges.

Decomposition: GRR
Decomposition: GRRN

§ GAUSSIAN LIKELIHOOD WITH RECTIFIED-NORMAL PRIORS (GRR) AND HIERARCHICAL PRIOR (GRRN)

Graphical representation of GRR and GRRN models. Green circles denote prior variables, orange circles represent observed and latent variables, and plates represent repeated variables. The slash “/" in the variable represents “or," and the comma “," in the variable represents “and."

Going further, the Gaussian likelihood with rectified-normal and hierarchical priors (GRR and GRRN) models are proposed in Lu and Ye, 2022 to  enhance flexibility building upon the GTT and GTTN models.
In this context, again, we view the data $\bA$ as being produced according to the probabilistic generative process shown in Figure <ref>. The observed $(m,n)$-th data  entry $a_{mn}$ of matrix $\bA$ is modeled using a Gaussian likelihood function with variance $\sigma^2$ and a mean given by the latent decomposition $\bw_m^\top\bz_n$ (Equation (<ref>)). The chosen likelihood is consistent with that employed in the GEE model (Equation (<ref>)).

Rectified-normal distribution
Exponentially rectified-normal distribution
Hierarchical prior
RN-scaled-normal-Gamma prior

We treat the latent variables $w_{mk}$'s (and $z_{kn}$'s) as random variables. And we need prior densities over these latent variables to express beliefs for their values, e.g., nonnegativity in this context.
Here, we assume further that the latent variables $w_{mk}$'s and $z_{kn}$'s are independently drawn from a rectified-normal (RN) priors (a.k.a., an exponentially rectified-normal distribution, Definition <ref>, p. definition:reftified_normal_distribution),

\begin{equation}\label{equation:rn_prior_grrn}
\begin{aligned}
p(w_{mk} \mid \cdot ) &= \rectifieddist(w_{mk} \mid \mu_{mk}^W, (\tau_{mk}^W)^{-1}, \lambda_{mk}^W);\\
p(z_{kn} \mid \cdot ) &= \rectifieddist(z_{kn} \mid \mu_{kn}^Z, (\tau_{kn}^Z)^{-1}, \lambda_{kn}^Z).
\end{aligned}
\end{equation}

This prior serves to enforce the nonnegativity constraint on the components $\bW$ and $\bZ$, and it is conjugate to the Gaussian likelihood (Equation (<ref>), p. equation:conjugate_rectified_nonnegative_mean).
In some scenarios, the two sets of latent variables can be drawn from two different rectified-normal priors, e.g., enforcing sparsity in $\bW$ while non-sparsity in $\bZ$. And we shall not consider this case for our later examples as it is not the main interest of this book.
The posterior density is a truncated-normal distribution, which is a special rectified-normal distribution.
The model is then called a Gaussian likelihood with rectified-normal priors (GRR) model. Since the RN distribution is a special TN distribution, the GRR model is the same as the GTT model with a careful choice of prior parameters. What makes the RN prior important is from the hierarchical model that provides flexibility and guidance on the prior parameter choices; see discussion below.

Hierarchical prior.
To further favor flexibility, we choose a convenient joint hyperprior density over the parameters $\{\mu_{mk}^W, \tau_{mk}^W, \lambda_{mk}^W\}$ of RN prior in Equation (<ref>), namely, the RN-scaled-normal-Gamma (RNSNG) prior,

\begin{equation}
\begin{aligned}
&\gap p(\mu_{mk}^W, \tau_{mk}^W, \lambda_{mk}^W \mid\cdot)
= \rnsng(\mu_{mk}^W, \tau_{mk}^W, \lambda_{mk}^W\mid \mu_\mu, \tau_\mu, a, b, \alpha_\lambda, \beta_\lambda)\\
&=	C(\mu_{mk}^W, \tau_{mk}^W, \lambda_{mk}^W)
\cdot
\normal(\mu_{mk}^W\mid \mu_\mu, (\tau_\mu)^{-1})
\cdot \gammadist(\tau_{mk}^W \mid a, b)
\cdot \gammadist(\lambda_{mk}^W \mid \alpha_\lambda, \beta_\lambda),
\end{aligned}
\end{equation}

where $C(\mu_{mk}^W, \tau_{mk}^W, \lambda_{mk}^W)$ is a constant in terms of $\{\mu_{mk}^W, \tau_{mk}^W, \lambda_{mk}^W\}$.
This prior can decouple parameters $\mu_{mk}^W, \tau_{mk}^W$, and $\lambda_{mk}^W$, and their posterior conditional densities are Gaussian, Gamma, and Gamma respectively due to this convenient scale. A similar RNSNG prior is given over $\{\mu_{kn}^Z, \tau_{kn}^Z, \lambda_{kn}^Z \}$.

Markov blanket

Again, following  Bayes' rule and MCMC, this means we need to be able to draw from distributions (by Markov blanket, Section <ref>, p. section:markov-blanket):

\begin{aligned}
&p(w_{mk}\mid\bA, \bW_{-mk}, \bZ, \sigma^2, \mu_{mk}^W,\tau_{mk}^W, \lambda_{mk}^W), \\
& p(z_{kn}\mid\bA, \bW, \bZ_{-kn}, \sigma^2,  \mu_{kn}^Z, \tau_{kn}^Z, \lambda_{kn}^Z), \\
& p(\sigma^2 \mid \bA, \bW, \bZ, \alpha_\sigma, \beta_\sigma ), \\
\end{aligned}

where $\bW_{-{mk}}$ denotes all elements of $\bW$ except $w_{mk}$, and $\bZ_{-kn}$ denotes all elements of $\bZ$ except $z_{kn}$.
Using Bayes' theorem,
the conditional density of $w_{mk}$ depends on its parents ($\mu_{mk}^W$, $\tau_{mk}^W$, $\lambda_{mk}^W$), children ($a_{mn}$), and coparents ($\tau$ or $\sigma^2$, $\bW_{-mk}, \bZ$) [See Figure <ref> and Section <ref> (p. section:markov-blanket).].
The conditional density of $w_{mk}$ follows a truncated-normal density.
And it can be obtained by (similar to computing the conditional density of $w_{mk}$ in the GEE model, Equation (<ref>))

\begin{equation}\label{equation:posterior_grrn_wmk122_app}
\small
\begin{aligned}
&\gap p(w_{mk}\mid \bA, \bW_{-mk}, \bZ,\sigma^2, \mu_{mk}^W,\tau_{mk}^W, \lambda_{mk}^W)
\propto p(\bA\mid\bW, \bZ, \sigma^2) \times p(w_{mk} \mid \mu_{mk}, \tau_{mk}, \lambda_{mk})\\
&\propto \prod_{i,j=1}^{M,N} \normal(a_{ij}\mid \bw_i^\top\bz_j, \sigma^2)
\times \rectifieddist(\mu_{mk}, (\tau_{mk})^{-1}, \lambda_{mk})\\
&\stackrel{\star}{\propto} \prod_{i,j=1}^{M,N} \normal(a_{ij}\mid \bw_i^\top\bz_j, \sigma^2)
\times \truncatednormal\left( \underbrace{\frac{\tau_{mk}^W\mu_{mk}^W - \lambda_{mk}^W}{\tau_{mk}^W}}_{\textcolor{mylightbluetext}{:=\mu^\prime}}, (\tau_{mk})^{-1} \right)\\
\exp\left\{ -\left(\frac{\sum_{j=1}^{N}z_{kj}^2 }{2\sigma^2} + \frac{\tau_{mk}^W}{2}\right)w_{mk}^2 + w_{mk} \left(
\frac{1}{\sigma^2} \sum_{j=1}^{N}z_{kj}(a_{mj}- \sum_{i\neq k}^{K}w_{mk}z_{ij}) + \tau_{mk}^W \mu^\prime
\right)  \right\} u(w_{mk})\\
& \propto \normal(w_{mk}\mid \widetilde{\mu_{mk}} , \widetilde{\sigma_{mk}^2})u(w_{mk})
= \truncatednormal(w_{mk}\mid \widetilde{\mu_{mk}} , \widetilde{\sigma_{mk}^2}),
\end{aligned}
\end{equation}

where the equality $(\star)$ is from the equivalence between the RN and TN distributions (Definition <ref>, p. definition:reftified_normal_distribution), $\widetilde{\sigma_{mk}^2} = \frac{\sigma^2}{ \sum_{j=1}^{N} z_{kj}^2 + \tau_{mk}^W \cdot \sigma^2 }$ is the posterior “parent" variance of the normal distribution with posterior “parent" mean

\widetilde{\mu_{mk}} =
\left(
\frac{1}{\sigma^2} \sum_{j=1}^{N}z_{kj}(a_{mj}- \sum_{i\neq k}^{K}w_{mk}z_{ij}) + \tau_{mk}^W \mu^\prime
\right)\cdot \widetilde{\sigma_{mk}^2}.

The  $\mu^\prime = \frac{\tau_{mk}^W\mu_{mk}^W - \lambda_{mk}^W}{\tau_{mk}^W}$ is the “parent" mean of the truncated-normal density.
Due to symmetry, the conditional posterior for $z_{kn}$ can be  derived similarly.

Extra update for GRRN.
Following   the graphical representation of the GRRN model in Figure <ref>,
we need to draw samples iteratively from

\begin{aligned}
&p(\mu_{mk}^W \mid\textcolor{black}{ \tau_{mk}^W, \lambda_{mk}^W}, \mu_\mu, \tau_\mu, a,  b,\alpha_\lambda, \beta_\lambda, w_{mk}),\\
&p(\tau_{mk}^W \mid\textcolor{black}{ \mu_{mk}^W, \lambda_{mk}^W}, \mu_\mu, \tau_\mu, a,  b,\alpha_\lambda, \beta_\lambda, w_{mk}),\\
&p(\lambda_{mk}^W \mid\textcolor{black}{ \mu_{mk}^W,\tau_{mk}^W} , \mu_\mu, \tau_\mu, a,  b,\alpha_\lambda, \beta_\lambda, w_{mk}).
\end{aligned}

The conditional density for $\mu_{mk}^W$ is a truncated-normal (a special rectified-normal):

\begin{equation}\label{equation:posterior_grrn_mu_tau1222_app}
\begin{aligned}
&\gap p(\mu_{mk}^W \mid\textcolor{black}{ \tau_{mk}^W, \lambda_{mk}^W}, \mu_\mu, \tau_\mu, a,  b,\alpha_\lambda, \beta_\lambda, w_{mk})\\
&\propto \rectifieddist(w_{mk} \mid \mu_{mk}^W, (\tau_{mk}^W)^{-1}, \lambda_{mk}^W) \cdot
\rnsng(\mu_{mk}^W, \tau_{mk}^W, \lambda_{mk}^W\mid \mu_\mu, \tau_\mu, a, b, \alpha_\lambda, \beta_\lambda)\\
&\propto \rectifieddist(w_{mk} \mid \mu_{mk}^W, (\tau_{mk}^W)^{-1}, \lambda_{mk}^W)
% \times C(\mu_{mk}^W, \tau_{mk}^W, \lambda_{mk}^W)
\cdot  \normal(\mu_{mk}^W\mid \mu_\mu, (\tau_\mu)^{-1})
\cdot \gammadist(\tau_{mk}^W \mid a, b)
\cdot  \gammadist(\lambda_{mk}^W \mid \alpha_\lambda, \beta_\lambda)\\
&= %\frac{1}{C^{RN}(\mu_{mk}^W, \tau_{mk}^W, \lambda_{mk}^W)}
\normal(w_{mk}| \mu_{mk}^W, (\tau_{mk}^W)^{-1})\cdot
\cancel{\exponential(w_{mk}| \lambda_{mk}^W)}
\cdot  \normal(\mu_{mk}^W| \mu_\mu, (\tau_\mu)^{-1})
\cdot \cancel{\gammadist(\tau_{mk}^W | a, b)}
\cdot  \cancel{\gammadist(\lambda_{mk}^W | \alpha_\lambda, \beta_\lambda)}\\
&\propto \normal(w_{mk}\mid \mu_{mk}^W, (\tau_{mk}^W)^{-1})\normal(\mu_{mk}^W\mid \mu_\mu, (\tau_\mu)^{-1})
%&\propto \exp\left( -\frac{\tau_{mk}^W + \tau_\mu}{2} (\mu_{mk}^W)^2 + \mu_{mk}^W
%\underbrace{(\tau_{mk}^W w_{mk} +\tau_{\mu}\mu_\mu)}_{\textcolor{mylightbluetext}{\widetilde{m}\cdot  \widetilde{t}}} \right)
\propto \normal(\mu_{mk}^W \mid \widetilde{m}, \widetilde{t}^{-1}),
\end{aligned}
\end{equation}

\widetilde{t}=\tau_{mk}^W + \tau_\mu$ and $
\widetilde{m}=(\tau_{mk}^W w_{mk} +\tau_{\mu}\mu_\mu)/\widetilde{t}
are the posterior precision and mean, respectively. The samples $w_{mk}$'s are nonnegative due to the rectification in the distribution (by exponential distribution with the density). However, this “parent" mean parameter $\mu_{mk}^W$ is not limited to be nonnegative.

The conditional density for $\tau_{mk}^W$ is a Gamma distribution:

\begin{equation}\label{equation:posterior_grrn_tau_tau1222_app}
\begin{aligned}
&\gap p(\tau_{mk}^W \mid\textcolor{black}{ \mu_{mk}^W, \lambda_{mk}^W}, \mu_\mu, \tau_\mu, a,  b,\alpha_\lambda, \beta_\lambda, w_{mk})\\
&\propto \rectifieddist(w_{mk} \mid \mu_{mk}^W, (\tau_{mk}^W)^{-1}, \lambda_{mk}^W)
\cdot \rnsng(\mu_{mk}^W, \tau_{mk}^W, \lambda_{mk}^W\mid \mu_\mu, \tau_\mu, a, b, \alpha_\lambda, \beta_\lambda) \\
&\propto \rectifieddist(w_{mk} \mid \mu_{mk}^W, (\tau_{mk}^W)^{-1}, \lambda_{mk}^W)
\cdot  \normal(\mu_{mk}^W\mid \mu_\mu, (\tau_\mu)^{-1})
\cdot \gammadist(\tau_{mk}^W \mid a, b)
\cdot  \gammadist(\lambda_{mk}^W \mid \alpha_\lambda, \beta_\lambda)\\
\normal(w_{mk}| \mu_{mk}^W, (\tau_{mk}^W)^{-1})\cdot
\cancel{\exponential(w_{mk}| \lambda_{mk}^W)}
\cdot  \cancel{\normal(\mu_{mk}^W| \mu_\mu, (\tau_\mu)^{-1})}
\cdot {\gammadist(\tau_{mk}^W | a, b)}
\cdot  \cancel{\gammadist(\lambda_{mk}^W | \alpha_\lambda, \beta_\lambda)}\\
&\propto \normal(w_{mk}\mid \mu_{mk}^W, (\tau_{mk}^W)^{-1})\gammadist(\tau_{mk}^W \mid a, b)\\
&\propto (\tau_{mk}^W)^{a+\frac{1}{2}-1} \exp\left\{-\left( b+ \frac{(w_{mk}-\mu_{mk}^W)^2}{2} \right) \tau_{mk}^W \right\}
\propto \gammadist(\tau_{mk}^W \mid \widetilde{a}, \widetilde{b}),
\end{aligned}
\end{equation}

where $\widetilde{a} = a+\frac{1}{2}$ and $\widetilde{b}= b+ \frac{(w_{mk}-\mu_{mk}^W)^2}{2}$ are the posterior shape and rate parameters, respectively.

Furthermore, the conditional density for $\lambda_{mk}^W$ follows also a Gamma distribution:

\begin{equation}\label{equation:posterior_grrn_lambda122_app}
\begin{aligned}
&\gap p(\lambda_{mk}^W \mid\textcolor{black}{ \mu_{mk}^W,\tau_{mk}^W} , \mu_\mu, \tau_\mu, a,  b,\alpha_\lambda, \beta_\lambda, w_{mk})\\
&\propto \rectifieddist(w_{mk} \mid \mu_{mk}^W, (\tau_{mk}^W)^{-1}, \lambda_{mk}^W)\cdot
\rnsng(\mu_{mk}^W, \tau_{mk}^W, \lambda_{mk}^W\mid \mu_\mu, \tau_\mu, a, b, \alpha_\lambda, \beta_\lambda) \\
&\propto \rectifieddist(w_{mk} \mid \mu_{mk}^W, (\tau_{mk}^W)^{-1}, \lambda_{mk}^W)
\cdot  \normal(\mu_{mk}^W\mid \mu_\mu, (\tau_\mu)^{-1})
\cdot \gammadist(\tau_{mk}^W \mid a, b)
\cdot  \gammadist(\lambda_{mk}^W \mid \alpha_\lambda, \beta_\lambda)\\
\cancel{\normal(w_{mk}| \mu_{mk}^W, (\tau_{mk}^W)^{-1})}\cdot
{\exponential(w_{mk}| \lambda_{mk}^W)}
\cdot  \cancel{\normal(\mu_{mk}^W| \mu_\mu, (\tau_\mu)^{-1})}
\cdot \cancel{\gammadist(\tau_{mk}^W | a, b)}