How to measure the information loss due to coarse-graining of a physical system into a graphical representation? Let's consider a system of bead-spring with $N+1$ beads connected with $N$ springs:

The Hamiltonian of such a chain is:
$$
\mathcal{H} = \frac{1}{2} k \sum_{i=1} ^N (\mathbf{r}_{i+1} - \mathbf{r}_{i})^2 + \frac{1}{2} m\mathbf{\dot{r}}_{i} 
$$
where $k$ is the spring constant, $m$ is the mass, $\mathbf{r}_{i} $ is the displacement, and $\mathbf{\dot{r}}_{i}$ is the velocity of the $i$th particle.
Assume now that now I create a graph $G=(V,E)$ (an ordered pair of nodes and edges) representation of the system as follows:
Let $\mathbf{D}$ be an $N+1 \times N+1$ matrix where $\mathbf{D}_{ij} = |\mathbf{r}_{i}- \mathbf{r}_{j}|$. Now we let matrix $\mathbf{A}$ be the adjacency matrix such that $\mathbf{A}_{ij}=1$ if $\mathbf{D}_{ij} < d_{0}$ or $0$ otherwise, where $d_{0}$ is some constant threshold.
Now I would like to quantify how much information was lost during the encoding of the physical system into the graph. I would love to hear some suggestions how to go around this. I assume that some mutual information related derivation should be considered here. I just would to mention that this question has a general implication in quantifying any lossy encoding algorithms including autoencoders which are very abundant today in machine learning.

Edit:
Might be related too: Fisher Information.
Maybe something with field theory, see Condensed Matter Field Theory.
 A: In order to have a thermodynamic entropy that makes sense, we need to use canonical coordinates. The Lagrangian for this system is
$$\mathcal{L} = \frac{1}{2}\sum_i m\dot{\mathbf{r}}_i^2 - k(\mathbf{r}_{i + 1} - \mathbf{r}_i)^2  $$
We can choose generalized coordinates $\mathbf{q_i} = \mathbf{r}_{i + 1} - \mathbf{r}_i$ so that
$$\mathbf{r}_i = \mathbf{r}_0 + \sum_{j < i} q_j$$
This will make the later steps easier, since $A$ will depend only on the $q_i$. We can also fix $r_0 = 0$ without changing anything important, I think (but this deserves further thought). In terms of these, the Lagrangian is
$$\mathcal{L} = \frac{1}{2}\sum_i  m(\sum_{j < i} \dot{q}_j)^2 - kq_i^2$$
and we can simplify the kinetic term as
$$ \frac{m}{2} \sum_i \sum_j w_{ij}\dot{q}_i\dot{q}_j $$
where the weights are some symmetric combinatorial coefficients which turn out to not matter. We can stack all the $\mathbf{q}_i$ into one big vector to get this as
$\frac{m}{2}\mathbf{\dot{q}}^T\mathbf{w}\mathbf{\dot{q}}$
The canonical momenta are
$$ \mathbf{p} = \nabla_{\mathbf{\dot{q}}} \mathcal{L} = m\mathbf{w}\mathbf{\dot{q}}  $$
so we can invert $\mathbf{w}$ to get the velocities from the momenta. The Hamiltonian is
$$\mathcal{H} = \frac{1}{2}kq_i^2 + \frac{m}{2}\mathbf{\dot{q}}^T\mathbf{w}\mathbf{\dot{q}}$$
or, in terms of the momenta,
$$\mathcal{H} = \frac{1}{2}k\mathbf{q}^T\mathbf{q} + \frac{m}{2}\mathbf{p}^T\mathbf{w}^{-1}\mathbf{p}$$
For protein folding, you probably consider the chain in thermal contact with some environment at fixed temperature $T$. The canonical ensemble is
$$ \rho = \frac{1}{Z}\exp\left(-\frac{1}{T}\mathcal{H}\right)$$
which is a multivariate Gaussian distribution in $\mathbf{x} = (\mathbf{q},\mathbf{p})$. The positions of the beads are thus a discretization of the (appropriately scaled) Wiener process, and in the large-$N$ limit this approaches the Wiener process itselt.
So far, this has all been setup to argue that the thermal distribution is what you'd expect. We are now ready to talk about information. We have a specific, known distribution $\rho$, and we wish to measure some property $A$. The goal is to maximize the information we get by measuring $A$. $A$ is fully determined by $\mathbf{x}$, and so $H(A|x) = 0$. By this Venn diagram, you can see that it thus suffices to maximize
$$ H(A) = -\sum p(A)\log p(A) $$
where the sum is over all possible $A$. The remaining problem is that of evaluating $H(A)$ (or its derivative with respect to $d_0$). Three possible strategies are:

*

*Find the distribution of $D$ analytically. For each possible configuration $A$, integrate $p(D)$ to find $p(A)$. These are $2^{O(N^2)}$-dimensional integrals, and there are $2^{O(N^2)}$ of them.

*Simulate many random walks numerically, and compute $A$ for each. Estimate $H(A)$ from the resulting distribution (using one of several known numerical entropy estimators). We need enough draws that some values of $A$ occur many times, so this requires at worst $2^{O(N^2)})$ samples. I think concentration of measure lets us get away with only $2^{O(N)}$, since only a certain diagonal band of $D$ will be likely to have values close to $d_0$. Symmetry arguments suggest that we can't do any better than that.

*Rewrite $H(A) = -\langle \log p(A) \rangle $. Simulate many random walks numerically, then compute $\log p(A)$ for each and average them. Because we get central-limit-theorem convergence, we need only $O(1)$ random walks. However, computing $p(A)$ still requires a $2^{O(N^2)})$ dimensional integral for each.

None of these look easy. However, if we can find an efficient way to do the integrals then the third option becomes tractable. Fix $d_0$ and let $S(A)$ be the subset of $D$-space consistent with a particular $A$. One option is to generate samples from some distribution $q(D)$ supported over only $S(A)$. Then
$$ \int_{S(A)} p(D)  = \left\langle \frac{p(D)}{q(D)}\right\rangle_{D \sim q}$$
(This is a standard Monte Carlo method). We need a $q$ which is easy to sample from, easy to evaluate, and we need to know $p(D)$. To reduce the variance, we also want $q(D)$ close to $p(D)$.
I'll start by finding $p(D)$. The marginal distributions of Gaussians are easy to work with. In particular, the marginal distribution over $\mathbf{q}_i$ is
$$ \frac{1}{Z}\exp\left(-\frac{k}{2T}\mathbf{q}^T\mathbf{q}\right)$$
The marginal distribution over $\sqrt{\mathbf{q}_i^2}$ for any single $i$ will be something like a Maxwell-Boltzmann distribution, depending on how many dimensions we are working in. The entries of $D$ are not independent, however, so we need to figure out the joint PDF of the $D_{ij}$.
We can pick some traversal order through the matrix elements, so that we label with a single index as $D_n$. We can then build the joint pdf in terms of conditional pdfs as
$$p(d_1...d_N) = p(d_1)p(d_2|d_1)p(d_3|d_1,d_2)...p(d_N|d_1...d_{N-1}) $$
If we choose a good traversal order, we can make sure the conditional pdfs are all known. For the diagonal elements, we have
$$ p(d_{ii}) = \delta(d_{ii})$$
so we can start with those. For the elements just below the diagonal, $ p(d_{i,i+1}) $ is a dimension-dependent generalization of the Maxwell-Boltzmann distribution. If the dimension is $m$, it's of the form
$$ p(d_{i,i+1}) = \frac{1}{Z} d_{i,i+1}^me^{-\frac{1}{2\sigma^2}d_{i,i+1}^2} $$
Call this probability density function $f(x)$.
For concreteness, let's think about how $D_{14}$ is determined for $N = 4$. Let $\theta$ be the angle between the line from $r_1$ to $r_3$ and the line from $r_3$ to $r_4$. It is distributed according to some $p_\theta$ which again depends only on dimension $m$. Using the law of cosines, we can write
$$ D_{14}^2 = D_{13}^2 + D_{34}^2 - 2D_{13}D_{34}\cos\theta $$
From this we can get the conditional distribution $ P(D_{14}|D_{13},D_{34} )$. $D_{14}$ is independent of all the other $D_{ij},i,j<4$ once we condition on $D_{13},D_{34}$, and so actually
$$ P(D_{14}|D_{13},D_{34}) = P(D_{14}|\{D_{ij},i,j < 4\})$$
Call this function $g(x|y,z)$. It will depend only on dimension.
Now we are ready to traverse the matrix. We start with the main diagonal and then move by diagonal stripes. We find
$$
p(\{d_{ij}\}) = \left(\prod_{i=1}^N \delta(d_{ii})\right)\left(\prod_{i=1}^{N-1} f(d_{i,i+1})\right)\left(\prod_{i=1}^{N-2} g(d_{i,i+2}|d_{i,i+1}, d_{i+1,i+2})\right)\left(\prod_{i=1}^{N-3} g(d_{i,i+3}|d_{i,i+2}, d_{i+2,i+3})\right) ... \left(\prod_{i=1}^{N-k} g(d_{i,i+k}|d_{i,i+k-1}, d_{i+k-1,i+k})\right)... g(d_{1,N}|d_{1,N-1},d_{N-1,N})
$$
where each parenthesized term accounts for one diagonal stripe. Now we know the full joint distribution of $D$ (assuming we've chosen a dimension and figured out what $f$ and $g$ are). This product has $O(N^2)$ terms in it, so it's relatively efficient to evaluate.
Next, we need to find some good distribution $q$. It's going to be easiest to have an independent distribution for each element of $D$. We know that the marginal distribution of $D_{ij}$ is a generalized Maxwell-Boltzmann distribution with mean proportional to $\sqrt{i-j}$. We can take the left-truncated form if $A_{ij} = 0$ and the right-truncated form if $A_{ij} = 1$ to obtain a distribution with support $S(A)$. The normalization requires access to the CDF associated with $f$, and sampling can be done either by rejection sampling or by using the inverse CDF. (I can come back and explain this part in more detail later if needed).
Summary
For a system in thermal equilibrium, the canonical ensemble provides a well-defined distribution with a meaningful entropy. We wish to choose $d_0$ to maximize the information we learn by measuring $A$. I haven't found an analytical formula for $H$ in terms of $d_0$, but there are some numerical methods that can be used to compute it. One of these numerical methods looks reasonably efficient.
