# Bayesian learning of tree distribution

this is my first post here.

I'm currently trying to compute the posterior predictive likelihood for a tree-structured distribution, following the paper Tractable Bayesian learning of tree belief networks by Marina Meila. Specifically, I'm asking about equation (46) in the paper:

$$P(x|\mathcal{D}) = \frac{w_0(x)\left|Q(\beta Ww(x))\right|}{\left|Q(\beta W)\right|}.$$

where $x$ is a $d$ dimensional discrete random variable, in my application, you can assume it to be binary-valued; $\mathcal{D}$ is the observed (marginal and bivariate marginal) statistics from data. $\beta$ is a $d$-by-$d$ real valued matrix denoting the prior distribution for the structure of the tree. $W$ is another $d$-by-$d$ matrix where each entry is computed from $\mathcal{D}$. $w(x)$ is also a $d$-by-$d$ matrix computed from each dimension of $x$ and bivariate interactions of $x$. The matrix multiplication here is entrywise matrix product. $Q$ is an operator that takes a $d$-by-$d$ matrix, compute its Laplacian, and take the first $d-1$-by-$d-1$ submatrix. Each entry of $W$, $\beta$, $w$ are all positive.

Theoretically, what $\left|Q(\beta W)\right|$ does is to use Matrix-Tree theorem to sum over all possible spanning trees for a weighted adjacency matrix $\beta W$.

In practice, $W$ may have a large range, the maximum entry might be $e^{-2300}$, minimum entry might be $e^{-2600}$, which means the edge weights in the adjacency matrix are close to 0.

When computing the Laplacian for $\beta W$, for summing over the edge weights, because some entries (e.g. $e^{-2300}$ is a lot larger than $e^{-2600}$), the diagonal sum will be (I'm doing computation in log domain) as same as the largest entry. Therefore, when exponentiating this log-entry for determinant calculation, there might be some large-off-diagonal entries making the determinant negative, while in theory the determinant should always be positive. Thus this posterior likelihood cannot be computed numerically stable.

I'm wondering whether someone has implemented equation (46) in Meila2006, and could share some experience.

Currently I have tried two ways of going around this numerical issue:

1. Rescale both $\beta W$ and $\beta W w(x)$, however, because the issue comes from the huge log-domain difference, this rescaling doesn't work.

2. I also tried doing an Eigen-decomposition of the Laplacian of $\beta W$, and found out that there might be some very small negative eigenvalues (like $-e^{-142}$), which makes the determinant negative. I then try to use the approach in Andrew Ng's 2001 NIPS paper ( On Spectral Clustering: Analysis and an algorithm) to do spectral clustering on $\beta W$, and decompose $\beta W$ into different connected components. But the kmeans algorithm gives some random cluster results based on the eigenvalues.

Thanks all in advance for any suggestions.