Computing entropy of tensor networks containing classical probability distributions

[apiedrafTNO]

Hi all,

we are looking into using tensor networks to encode classical probability distributions. The goal is to build an inference engine that will replace the back-end of a machine diagnostics library.
One crucial step for diagnostics is computing the entropy of joint distributions. The naïve approach that we've so far used is to take a network containing (implicitly) our joint distribution, contracting all the inner indices, and then computing the entropy of the resulting tensor array, but we're starting to need better speed.
There are some entropy-related methods available in quimb, but they seem to only apply to structured networks like MPS or matrix-like networks, with left and right indices, rather than taking the network as a multilinear operator.

My question is: is there any clever way of computing the largest k eigenvalues of a general non-negative tn?
What comes to mind is that this non-negative tn with a bunch of open indices is kinda-sorta like an unnormalized pure state, and what I want to do is compute a truncated Schmidt decomposition.
I have a background in quantum information, but not in condensed matter or many-body physics. I know that DMRG is used to compute eigenvalues of a Hamiltonian, but I doubt it would serve here since the TNs we get aren't usually even planar, let alone 1D.

Does anybody have any ideas?

[jcmgray]

Hi @apiedrafTNO, I am not very familiar with the language of inference and machine diagnostics, so maybe you could give some more specific details to narrow it down:

1. Do you a have tensor network description of the (un-normalized?) prob dist directly? (i.e. one output index per variable)? In this case which are the left and right degrees of freedom for the Schmidt decomposition, or do you simply need $\sum_{\{x\}} p_x \log p_x$ ?
2. Or do you have some "density matrix" where the probability distribution is the diagonal? (i.e. a bra and ket output index per variable) so that $Tr \rho \log \rho$ is what you need?

There are possible methods for both which can get past computing the full output tensor, but non-linear quantities are trickier to scale up.

[apiedrafTNO]

Hi @jcmgray,
thanks for such a quick response!

It is definitely the first case.
The link between machine diagnostics and tensor networks is through probabilistic graphical models. https://arxiv.org/pdf/1710.01437.pdf This paper from 2017 makes that connection explicitly. The idea is that the probabilistic graphical model is a graphical representation of a factorizable joint probability distribution that looks like $$p(x_1,\dots,x_n)=\frac{1}{Z}\prod_{c\in\mathcal{C}}\psi_c(x_c)$$ , where $c$ are subsets of variables.
This formula is functionally the same as that of the tensor represented by a tensor network before contraction. the variables in the graphical model are on the edges of the network, and the factors in the factorization are the tensors. Computing marginals is nothing but contracting all but one edge. Our tensor networks end up looking something like the example one in https://quimb.readthedocs.io/en/latest/tensor-basics.html#graph-orientated-tensor-network-creation, albeit with more internal structure.
What I'm interested in is indeed computing $\sum_{{x}}p_x\log p(x)$ where x is a subset of the indices of the network. There are no left and right degrees of freedom, so I'm guessing the Schmidt decomposition is not quite what I want.

The non-linearity is a bit of an issue, indeed. After all, from a functional analysis pov, what I actually want is to compute the expected value of $\log{T}$, where $T$ is a tensor resulting from contracting some edges of a network, but not others. That is why I was thinking that perhaps it would be easier to understand the tensor as a multilinear operator, compute the first $k$ singular values, and then take the entropy of that 1-d array.

I don't know if I'm making any sense. Perhaps it is quite a tall order. Anyway, I hope this is enough information.

[jcmgray]

So for classical tensor networks, I don't know of any really scalable methods, but you can push a bit further by contracting the output chunks lazily and computing $\sum p \log p$ on each of them separately - I added a cotengra example here https://cotengra.readthedocs.io/en/latest/examples/ex_large_output_lazy.html that shows this for a size $2^{36}$ marginal.

Beyond that probably sampling or using a replica trick might be possible, but nothing non-trivial comes to mind!