##  1 Estimating information: Theory

###  1.1 Overview

Using the minimal probabilistic model of an optical system described in
Section LABEL:main_paper_intro_to_model, the object-relevant information in
the noisy measurements can be quantified. Here, $\mathbf{X}$ is a random 2D
array of pixels with distribution $p(\mathbf{x})$, where each element
quantifies the expected amount of energy (in photons) at a pixel of a
noiseless image formed on a sensor. $\mathbf{Y}$ is a random 2D array of the
same shape with distribution $p(\mathbf{y})$, where each element quantifies
the amount of energy at that pixel in the noisy measurement.

The objective is to estimate $I(\mathbf{X};\mathbf{Y})$, the mutual
information between the noiseless image and the noisy measurement. This
quantity has units of bits, which has the operational significance that for
each bit collected, $2^{I(\mathbf{X};\mathbf{Y})}$ different objects can be
perfectly discriminated in the presence of noise (i.e. with no overlap in
their distributions of noisy measurements). In practice, however, different
objects cannot be unambiguously identified in the presence of noise, so each
bit quantifies multiple states that can be partially discriminated and how
well.

Optical systems often exhibit varying performance across their field of view,
both in resolution and in information content. However, it is often
mathematically and conceptually convenient to ignore their field-varying
optical performance and instead model them as shift-invariant, thereby
enabling a single global quantification of their resolution
[goodmanIntroductionFourierOptics2005]. Similarly, it is a convenient
simplification to express the variation in information across the field of
view with a single scalar quantity. To do so, we model the random variables
$\mathbf{X}$ and $\mathbf{Y}$ as stationary stochastic processes, meaning
their probability distributions are invariant to translations across the field
of view. The information in such a model is quantified by the mutual
information rate, which is the mutual information per pixel.

###  1.2 Mutual information estimator

Mutual information estimation is a well-studied problem in many fields,
including machine learning, neuroscience, and computational biology, and a
number of different approaches have been proposed, some of which attempt to
estimate mutual information directly [kraskovEstimatingMutualInformation2004a,
darbellayEstimationInformationAdaptive1999,
rossMutualInformationDiscrete2014], and others which try to infer its value by
computing upper or lower bounds [mcallesterFormalLimitationsMeasurement2018,
pooleVariationalBoundsMutual2019]. It is a challenging statistical problem to
solve in general, and many estimators and bounds are known to suffer from high
bias and/or variance, particularly on high-dimensional problems
[mcallesterFormalLimitationsMeasurement2018,
pooleVariationalBoundsMutual2019].

Many mutual information estimators rely on first estimating entropy, which is
a measure of the uncertainty of a random variable (Fig.
LABEL:fig2_informaiton_estimatora). The outcome of random variables with
higher entropy are more uncertain and harder to predict than the outcome of
those with lower entropy, because they are more random. Equivalently, the
higher a distribution’s entropy is, the more spread out it is, and the harder
it would be to digitally compress samples from it. Mathematically, entropy is
defined as the expected value of the negative log of the probability
density/mass function:

$\displaystyle\mathrm{H}\left(\mathbf{Y}\right)=-\mathbb{E}\left[\log
p(\mathbf{Y})\right]$ (1)

Like mutual information estimation, there is a large body of literature on
entropy estimation and many different approaches have been proposed
[beirlantNonparametricEntropyEstimation1997], including those that form an
estimate of the probability density function $\hat{p}(\mathbf{y})$ and plug it
into the definition of entropy, and those that estimate entropy based on the
similarity of samples from $p(\mathbf{y})$
[kozachenkoSampleEstimateEntropy1987,
lombardiNonparametricKnearestNeighbour2016,
jiaoNearestNeighborInformation2018]. Both approaches, however, face
difficulties in high dimensions.

There are multiple ways of decomposing mutual information into a difference of
entropies. Our approach is based upon the decomposition (Fig.
LABEL:fig2_informaiton_estimatorb):

$\displaystyle
I(\mathbf{X};\mathbf{Y})=\mathrm{H}\left(\mathbf{Y}\right)-\mathrm{H}\left(\mathbf{Y}\mid\mathbf{X}\right)$
(2)

Here, $\mathrm{H}\left(\mathbf{Y}\right)$ is the entropy of the noisy image
distribution $p(\mathbf{y})$ (Fig. LABEL:fig2_informaiton_estimatorc). Both
variations in the object and measurement noise contribute to the randomness of
measurements. We are interested in the information about the object, thus to
quantify how much of these variations the measurements contain, we must
subtract the entropy contributed by the noise,
$\mathrm{H}\left(\mathbf{Y}\mid\mathbf{X}\right)$ (Fig.
LABEL:fig2_informaiton_estimatord).

This decomposition is particularly useful for the problem of estimating mutual
information in optical imaging, because the conditional entropy
$\mathrm{H}\left(\mathbf{Y}\mid\mathbf{X}\right)$ can be greatly simplified
due to the fact that commonly used analytic noise models in optical imaging
can generally be computed for each pixel independently.

(A technical note: the outcomes of the random variables $\mathbf{X}$ and
$\mathbf{Y}$ in our model are discrete, because they are 2D arrays of pixels,
where each pixel takes on an intensity value that is an integer number of
photons. However, for computational simplicity, we will approximate them as
continuous random variables and use differential entropy instead of discrete
entropy [coverElementsInformationTheory2005,
mackayInformationTheoryInference2003]. These approximations break down at very
low photon counts, less than $\sim 20$ photons, so in this paper, we use only
data with photon counts greater than this. A possible direction of future work
is to extend our approach to work with discrete random variables, which would
allow it to be applied to data with lower photon counts.)

We begin by describing the simpler of the two terms to estimate, the
conditional entropy $\mathrm{H}\left(\mathbf{Y}\mid\mathbf{X}\right)$.

###  1.3 Estimating Conditional Entropy of Measurement Noise

The conditional entropy of noisy measurements given noiseless images can be
written out in terms of expectations over the logarithm conditional
probability of noisy measurements given a noiseless image
$p(\mathbf{y}\mid\mathbf{x})$:

$\displaystyle\mathrm{H}\left(\mathbf{Y}\mid\mathbf{X}\right)$
$\displaystyle=-\mathbb{E}_{\mathbf{X}}\left[\mathbb{E}_{\mathbf{Y}}\left[\log
p(\mathbf{Y}\mid\mathbf{X})\right]\right]$ (3)

$p(\mathbf{y}\mid\mathbf{x})$ embodies the various sources of noise in the
detection process, including photon shot noise, detector read noise, etc.
Here, we utilize established analytic models of detection noise in optical
imaging [goodmanStatisticalOptics2015]. Empirical results suggest that the
true noise in systems deviates from these models in low-light conditions
[monakhovaDancingStarsVideo2022], but since our experiments are conducted in
the high-light regime, we will assume that these models are accurate for the
purposes of this paper. A possible direction for future work is to learn more
accurate noise models from the data, as was done in
[krullProbabilisticNoise2VoidUnsupervised2019].

Assuming we have access to a dataset of $N$ samples from the distribution of
noiseless images
$\\{\mathbf{x}^{\left(1\right)},\mathbf{x}^{\left(2\right)},\dots\mathbf{x}^{\left(N\right)}\\}$,
the outer expectation can be estimated through Monte Carlo approximation:

$\displaystyle-\mathbb{E}_{\mathbf{X}}\left[\mathbb{E}_{\mathbf{Y}}\left[\log
p(\mathbf{Y}\mid\mathbf{X})\right]\right]$
$\displaystyle\approx-\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}_{\mathbf{Y}}\left[\log
p(\mathbf{Y}\mid\mathbf{x}^{\left(i\right)})\right]$ (4)
$\displaystyle=-\frac{1}{N}\sum_{i=1}^{N}\mathrm{H}\left(\mathbf{Y}\mid\mathbf{x}^{\left(i\right)}\right)$
(5)

Here, $\mathrm{H}\left(\mathbf{Y}\mid\mathbf{x}^{\left(i\right)}\right)$ is
the conditional entropy of the distribution of noisy images given the $i$th
noiseless image.

##### Conditionally independent noise at each pixel

$\mathrm{H}\left(\mathbf{Y}\mid\mathbf{x}^{\left(i\right)}\right)$ is a
function of the noise introduced in the detection process, which is modeled by
the probability distribution $p(\mathbf{y}\mid\mathbf{x}^{\left(i\right)})$.
In optical imaging, common analytic noise models like additive Gaussian and
Poisson shot noise typically assume that the noise at each pixel is
conditionally independent of the noise at other pixels, given the true
(noiseless) pixel value. When this is true,
$p(\mathbf{y}\mid\mathbf{x}^{\left(i\right)})$, which is a joint distribution
over all $D$ pixels in the noisy measurement, can be simplified by factoring
it into a product of scalar distributions for each pixel, where $y_{k}$ and
$x_{k}^{\left(i\right)}$ are the intensity values at the $k$th pixel in the
noisy measurement and the $i$th noiseless image, respectively:

$\displaystyle p(\mathbf{y}\mid\mathbf{x})$ $\displaystyle=p(y_{1},y_{2},\dots
y_{D}\mid x_{1}^{\left(i\right)},x_{2}^{\left(i\right)},\dots
x_{D}^{\left(i\right)})$ $\displaystyle=\prod_{k=1}^{D}p(y_{k}\mid
x_{k}^{\left(i\right)})$

This factorization simplifies the calculation of conditional entropy, because
it is much easier to compute $D$ scalar conditional entropies than a single
joint conditional entropy over $D$ variables. Mathematically, this
simplification can be seen by plugging the factorized distribution into the
definition of conditional entropy:

$\displaystyle\mathrm{H}\left(\mathbf{Y}\mid\mathbf{x}^{\left(i\right)}\right)$
$\displaystyle=-\mathbb{E}_{\mathbf{Y}}\left[\log\prod_{k=1}^{D}p(y_{k}\mid
x_{k}^{\left(i\right)})\right]$ (6)
$\displaystyle=-\mathbb{E}_{\mathbf{Y}}\left[\sum_{k=1}^{D}\log p(y_{k}\mid
x_{k}^{\left(i\right)})\right]$ (7)
$\displaystyle=-\sum_{k=1}^{D}\mathbb{E}_{\mathbf{Y}}\left[\log p(y_{k}\mid
x_{k}^{\left(i\right)})\right]$ (8)
$\displaystyle=\sum_{k=1}^{D}\mathrm{H}\left(\mathrm{Y}_{k}\mid
x_{k}^{\left(i\right)}\right)$ (9)

Here, $H(\mathrm{Y}_{k}\mid x_{k}^{\left(i\right)})$ is the conditional
entropy of the $k$th pixel in the noisy image given the intensity of the $k$th
pixel in the noiseless image.

This is a scalar quantity and can be calculated analytically for many common
noise models. We will discuss two such models here: additive Gaussian noise
and Poisson noise.

####  1.3.1 Conditional entropy with additive Gaussian noise

Additive Gaussian noise is a simple noise model often used in optical images,
especially in low-light conditions where the read noise of the detector is the
dominant source of noise. In this model, the noise at each pixel is drawn from
a Gaussian distribution with mean zero and variance $\sigma^{2}$.
Mathematically:

$\displaystyle\mathrm{Y}_{k}=\mathrm{X}_{k}+\mathrm{N}_{k}$
$\displaystyle\mathrm{N}_{k}\sim\mathcal{N}(0,\sigma^{2})$

The entropy of a (scalar) Gaussian distribution $\mathcal{N}(\mu,\sigma)$ is
[coverElementsInformationTheory2005]:

$\displaystyle\mathrm{H}\left(\mathrm{N}_{k}\right)=\frac{1}{2}\log_{2}(2\pi
e\sigma^{2})$

Since the noise is independent of the noiseless image, the conditional entropy
of the noise at each pixel is the same, and the full conditional entropy
(Equation 9) simplifies to:

$\displaystyle\mathrm{H}\left(\mathbf{Y}\mid\mathbf{x}^{\left(i\right)}\right)$
$\displaystyle=\sum_{k=1}^{D}\mathrm{H}\left(\mathrm{Y}_{k}\mid
x_{k}^{\left(i\right)}\right)$ (10)
$\displaystyle=\sum_{k=1}^{D}\mathrm{H}\left(\mathrm{N}_{k}\right)$ (11)
$\displaystyle=\frac{D}{2}\log_{2}(2\pi e\sigma^{2})$ (12)

Plugging this result into equations 5 and 3 yields:

$\displaystyle\mathrm{H}\left(\mathbf{Y}\mid\mathbf{X}\right)=\frac{D}{2}\log_{2}(2\pi
e\sigma^{2})$ (13)

To summarize, the conditional entropy of the noisy image given the noiseless
image is a constant, independent of the intensity values of the noiseless
images, and is equal to the number of pixels in the image times the entropy of
the noise distribution at each pixel.

####  1.3.2 Conditional entropy with shot noise

Images with high photon counts are limited in their signal-to-noise ratio by
the inherently random arrival times of photons, due to the quantum nature of
light. This noise is called shot noise, and it follows a Poisson distribution
with rate parameter equal to the expected number of photons at that pixel
[goodmanStatisticalOptics2015]. As the rate parameter of a Poisson
distribution increases, it becomes more and more similar to a Gaussian
distribution with equal mean and variance. In regimes in which shot noise is
the dominant source of noise, this noise can thus be very accurately
approximated with a Gaussian distribution with equal mean and variance
[goodmanStatisticalOptics2015]. Thus, the conditional entropy of the noise at
pixel $k$ for the $i$th noiseless image can be approximated as:

$\displaystyle\mathrm{H}\left(\mathrm{Y}_{k}\mid
x_{k}^{\left(i\right)}\right)=\frac{1}{2}\log_{2}(2\pi
ex_{k}^{\left(i\right)})$ (14)

Once again making use of the fact that the measurement noise at each pixel is
independent of the noise at other pixels conditional on the intensity of the
noiseless image at that pixel (Section 1.3), we can write the conditional
entropy for noiseless image $\mathbf{x}^{\left(i\right)}$ as:

$\displaystyle\mathrm{H}\left(\mathbf{Y}\mid\mathbf{x}^{\left(i\right)}\right)$
$\displaystyle=\sum_{k=1}^{D}\mathrm{H}\left(\mathrm{Y}_{k}\mid
x_{k}^{\left(i\right)}\right)$ (15)
$\displaystyle=\sum_{k=1}^{D}\frac{1}{2}\log_{2}(2\pi
ex_{k}^{\left(i\right)})$ (16)

The full conditional entropy (Equation 9) simplifies to:

$\displaystyle\mathrm{H}\left(\mathbf{Y}\mid\mathbf{X}\right)\approx\frac{1}{N}\sum_{i=1}^{N}\sum_{k=1}^{D}\frac{1}{2}\log_{2}(2\pi
ex_{k}^{\left(i\right)})$ (17)

To summarize, the conditional entropy under a Poisson noise model can be
approximated as a sum of the log of the intensity values of the noiseless
image at each pixel, averaged over $N$ noiseless images. This approximation is
accurate when the photon counts are high, and breaks down when the photon
counts are low, and is discussed further in Section 2.5.

###  1.4 Estimating Entropy of Noisy Images

The second, more challenging term to estimate in the decomposition of mutual
information (Equation 2) is the entropy of the noisy images,
$\mathrm{H}\left(\mathbf{Y}\right)$. It is challenging because unlike the
conditional entropy estimates in Section 1.3, the relevant joint probability
distribution $p(y)$ cannot be factorized into a product of independent
distributions for each pixel. That is, the pixels in a noisy measurement (like
those of a noiseless image) are not independent of each other.

Our approach for computing this entropy is to compute an upper bound by
fitting a probabilistic model to it, which empirically has been shown to be
more accurate than other types of entropy bounds in high dimensions
[mcallesterFormalLimitationsMeasurement2018].

Specifically, we create a parametric model for the distribution of noisy
images $p_{\mathbf{\theta}}(\mathbf{y})$, where $\theta$ are the parameters of
our model. We estimate these parameters $\theta$ by fitting the model to
samples from the empirical distribution of noisy images ${p}(\mathbf{y})$
using maximum likelihood estimation. Mathematically, the maximum likelihood
fit $\hat{\theta}_{\mathrm{MLE}}$ is found as the solution to the following
optimization problem, which uses the negative log likelihood of the data as
the objective function:

$\displaystyle\hat{\theta}_{\mathrm{MLE}}=\operatorname*{arg\,min}_{\theta}-\mathbb{E}\left[\log
p_{\theta}(\mathbf{Y})\right]$ (18)

This loss function, $-\mathbb{E}\left[\log p_{\theta}(\mathbf{Y})\right]$, is
also known as the cross entropy between the model distribution
$p_{\theta}(\mathbf{y})$ and the empirical distribution $p(\mathbf{y})$.

In practice, it is fit using a dataset of $N$ samples from the empirical
distribution $p(\mathbf{y})$:

$\displaystyle\hat{\theta}=\operatorname*{arg\,min}_{\theta}-\frac{1}{N}\sum_{i=1}^{N}\log
p_{\theta}(\mathbf{y}^{\left(i\right)})$ (19)

When the model distribution $p_{\theta}(\mathbf{y})$ is identical to the
empirical distribution $p(\mathbf{y})$, the data has been fit perfectly, and
the value of the cross entropy loss function is equal to the entropy of the
noisy images, $\mathrm{H}\left(\mathbf{Y}\right)$.

$\displaystyle-\mathbb{E}\left[\log p_{\theta}(\mathbf{Y})\right]$
$\displaystyle=-\mathbb{E}\left[\log p(\mathbf{Y})\right]$ (20)
$\displaystyle=\mathrm{H}\left(\mathbf{Y}\right)$ (21)

In practice, the model will not be able to fit the true distribution exactly,
and the average value of the loss function will be greater than the entropy of
the noisy images. The gap between the entropy that we are interested in
estimating and the cross-entropy loss function is given by the Kullback-
Leibler divergence between the model distribution and the empirical
distribution:

$\displaystyle-\mathbb{E}\left[\log
p_{\theta}(\mathbf{Y})\right]=\mathrm{H}\left(\mathbf{Y}\right)+D_{KL}(p_{\theta}\;\|\;p)$
(22)

The Kullback-Leibler divergence is a measure of the difference between two
probability distributions. It is always non-negative and is zero only when the
two distributions are identical. Thus, the cross-entropy loss function is an
upper bound on the entropy of the noisy images. The better our model fits the
data, the tighter this bound will be. Finding the right model that balances
the accuracy of this bound with the computational cost of fitting the data is
an important choice that is discussed further in Section 2.2.

In practice, the cross entropy loss function is evaluated on separate test set
of samples from the empirical distribution, to avoid overfitting to a subset
of the data and generating a model that is overly optimistic about its
performance.

This process is mathematically equivalent to the problem of data compression
in information theory111Technically, this is only completely true when the
model distribution is a discrete distribution, because continuous data cannot
be losslessly compressed in a discrete setting. For example, it would take an
infinite number of bits on a computer to represent an arbitrary real number
exactly. But the intuition remains the same.. In this problem, the goal is to
compress data by mapping each outcome to a string of bits, such that the total
number of bits needed on average is as small as possible. The optimal
compression scheme maps more probable outcomes to shorter bit strings, and
less probable outcomes to longer bit strings. The efficiency with which this
can be done is determined by how random the true distribution is (i.e. its
entropy).

###  1.5 Probabilistic models

In this section, we describe two different probabilistic models
$p_{\theta}(\mathbf{y})$ that we will use to fit the empirical distribution of
noisy measurements $p(\mathbf{y})$, and thus estimate via an upper bound the
entropy of noisy measurements, as described in Section 1.4. The two models
differ in their complexity, the tightness with which they upper bound the true
entropy, and the computational resources and data needed to fit them.

The first model is based on a multivariate Gaussian distribution. It is
simpler and faster to fit, but it is less flexible in its ability to fit the
data, and thus provides a looser upper bound on the entropy of the noisy
images.

The second model, a PixelCNN, is a type of model called an autoregressive
model, and is based on a neural network. It is more complex and
computationally expensive to fit in terms of time and data, but it is more
flexible in its ability to fit the data and can thus provide a tighter upper
bound on the entropy of the noisy images given sufficient resources.

There are many other possible models that could be used to fit the data using
this general procedure. Any model for which the likelihood function can be
directly evaluated can be used. There are several classes of such models based
on neural networks that have been studied in recent years, including other
types of auto-regressive models like transformers
[vaswaniAttentionAllYou2017], normalizing flows
[kobyzevNormalizingFlowsIntroduction2021], and energy-based models such as
diffusion models [sohl-dicksteinDeepUnsupervisedLearning2015,
songHowTrainYour2021].

Notably, the perceptual quality of samples is not necessarily an indicator
that it fits the data in terms of cross entropy loss
[theisNoteEvaluationGenerative2016], though state-of-the-art models in terms
of perceptual quality of samples can be adapted to perform well for both
[kingmaVariationalDiffusionModels]. We chose the PixelCNN model for its
relative simplicity and computational efficiency and because it has been shown
to produce good fits in terms of cross-entropy of test data
[salimansPixelCNNImproving2017]. However, it seems quite likely that there
already is or will be a better model for this task. One that improves on our
approach in terms of data and compute requirements and tightness of the upper
bound on the entropy of the noisy images.

####  1.5.1 Stationary stochastic processes and entropy rate

In order to be able to quantify the performance over the full field of view of
an imaging system by a single scalar quantity, we make a simplifying
assumption that the noiseless image and the noisy image are stationary
stochastic processes. This is the probabilistic analog of the assumption of a
constant optical resolution across the field of view
[goodmanIntroductionFourierOptics2005].

A stationary stochastic process with distribution $p(\mathbf{y})$ is one in
which the joint distribution of any set of pixels is invariant to translations
across the field of view. Focusing on a 1D stationary process (i.e. a single
row of pixels in an image) for simplicity, this means that the joint
distribution of any set of pixels depends only on their relative locations in
the image, not on their absolute locations. Mathematically, this means that
for any vector of pixels $(y_{1},y_{2},\dots,y_{D})$ with joint distribution
$p(y_{1},y_{2},\dots,y_{D})$, the joint distribution of pixels with the same
relative locations but all offset by a constant amount
$y_{1+k},y_{2+k},\dots,y_{D+k}$ is the same:

$\displaystyle p(y_{1},y_{2},\dots,y_{D})=p(y_{1+k},y_{2+k},\dots,y_{D+k})$
(23)

Practically, this means that we can view cropped patches from different parts
of the field of view as samples from the same distribution. This simplifies
the model and reduces the number of parameters that need to be fit.

Stationary stochastic processes have the useful property that they have a
constant entropy per pixel, called the entropy rate
[coverElementsInformationTheory2005] (explained in Figure 8 of
[pinkardVisualIntroductionInformation2022]). This entropy rate can be
equivalently viewed as the limit of the joint entropy
$H(y_{1},y_{2},\dots,y_{D})$ as the number of pixels $D$ goes to infinity
divided by $D$, or as the conditional entropy of the $D$th pixel given the
previous $D-1$ pixels as $D$ goes to infinity. Mathematically:

$\displaystyle\lim_{D\to\infty}\frac{1}{D}H(y_{1},y_{2},\dots,y_{D})=\lim_{D\to\infty}\mathrm{H}\left(y_{D}\mid
y_{1},y_{2},\dots,y_{D-1}\right)$ (24)

This can be understood in practical terms by considering the conditional
entropy of the $D$th pixel given the previous $D-1$ pixels, which quantifies
how random the $D$th pixel is given the previous $D-1$ pixels. As $D$
increases, we have seen more pixels of the image, and thus predicting the next
pixel becomes easier, and the conditional entropy decreases. Essentially, it
can be much easier to predict the value of a pixel in an image if you know the
values of neighboring pixels.

Since in practice we will not have images with an infinite number of pixels,
the models that we fit to the data will be finite length approximations of the
true underlying stochastic process. In general, modeling more pixels at once
will be able to fit more complex dependencies in the joint distribution, and
thus provide a tighter upper bound on the entropy of the noisy images.

####  1.5.2 Stationary Gaussian Process

In the stationary Gaussian model, the true distribution of noisy images is
approximated with a multivariate Gaussian distribution. A multivariate
Gaussian distribution is a generalization of the univariate Gaussian
distribution to multiple dimensions. It is fully specified by a mean vector
$\boldsymbol{\mu}$ and a covariance matrix $\mathbf{\Sigma}$. Its probability
density function is given by:

$\displaystyle
p_{\theta}(\mathbf{y})=\frac{1}{\sqrt{(2\pi)^{D}\det(\mathbf{\Sigma})}}\exp\left(-\frac{1}{2}(\mathbf{y}-\boldsymbol{\mu})^{T}\mathbf{\Sigma}^{-1}(\mathbf{y}-\boldsymbol{\mu})\right)$
(25)

Compared to an arbitrary multivariate Gaussian process, the stationary
Gaussian process has additional constraints on its mean vector and covariance
matrix, which simplify the model by reducing the number of parameters that
need to be fit. Mean vectors for non-stationary Gaussian processes can have
arbitrary values, but for stationary Gaussian processes they must be constant
vectors, reducing the number of parameters from the number of pixels $D$ to
just 1.

Covariance matrices for non-stationary Gaussian processes can be any
symmetric, positive definite matrix. For stationary Gaussian processes, the
covariance between any two pixels is only a function of their relative
positions in the image, not their absolute positions. For a 1-Dimensional
stationary Gaussian process (i.e. a single row of pixels in the image), the
matrix must be a Toeplitz matrix: a matrix in which each descending diagonal
from left to right is constant. For a 2-Dimensional stationary Gaussian
process (i.e. the full image), the matrix is a doubly Toeplitz matrix: It is
made up of submatrix blocks of Toeplitz matrices, and these blocks are
themselves arranged in a Toeplitz pattern (Fig. 1). These symmetric
Toeplitz/doubly Toeplitz matrices have only $D$ parameters, since the full
matrix can be constructed by repeating the values in the first column. This is
less than the $\frac{D(D+1)}{2}$ parameters needed to specify an arbitrary
symmetric, positive definite matrix.

Figure 1: The doubly Toeplitz $\mathbf{D\times D}$ covariance matrix of a
stationary 2D Gaussian process.

####  1.5.3 PixelCNN

The PixelCNN [oordPixelRecurrentNeural2016,
vandenoordConditionalImageGeneration2016, salimansPixelCNNImproving2017] is a
type of autoregressive probabilistic model that uses neural networks to model
complex dependencies between pixels.

Autoregressive models in general are a class of probabilistic models that
factorize a joint distribution into a product of conditional distributions,
and then model each of the conditional distributions. Mathematically:

$\displaystyle p(\mathbf{y})$ $\displaystyle=\prod_{k=1}^{D}p(y_{k}\mid
y_{1},y_{2},\dots,y_{k-1})$ (26) $\displaystyle=p(y_{1})p(y_{2}\mid
y_{1})p(y_{3}\mid y_{1},y_{2})\dots p(y_{M}\mid y_{1},y_{2},\dots,y_{D-1})$
(27)

This factorization does not require specific assumptions about the joint
distribution, so in theory it can be applied to any distribution.

Creating the full model thus requires creating $D$ models for each of the
conditional distributions. Each conditional distribution model takes between
$1$ and $D-1$ arguments, which are the values of the previous pixel of a
particular image, and outputs a 1D probability distribution over the possible
values of the next pixel.

While it is possible to implement this approach with $D$ different models for
each conditional distribution, the PixelCNN instead uses a specific type of
convolutional neural network that enables all the conditionals to be modeled
simultaneously, thus dramatically reducing computational complexity.

####  1.5.4 Sampling and likelihood evaluation

The primary purpose of fitting probabilistic models to data is to be able to
evaluate the likelihood of data under the model and compute an upper bound of
the data distribution entropy, as described in Section 1.4. However, it is
also useful to generate samples from the model, which can be used to visualize
the model and qualitatively assess its fit to the data. The visual quality of
samples and the likelihood of samples from the data distribution are
correlated for simple models like the stationary Gaussian process, and thus
inspecting quality of the samples can be a useful way to assess the quality of
the model fit. However, this is not necessarily true for more complex models
like the PixelCNN [theisNoteEvaluationGenerative2016]. Nonetheless, in these
cases it can still be useful to visualize samples to better understand its fit
to the data.

Both the PixelCNN and the stationary Gaussian process models are trained on
image patches of a fixed size $\sqrt{D}\times\sqrt{D}$. These models are thus
incapable of capturing dependencies between pixels that are larger than this
patch size. Increasing the patch size will give the models the ability to
capture longer range dependencies, and thus more accurately fit the data.

However, despite the fixed patch size, since the models are both stationary,
they can be used to evaluate likelihood and generate samples for images of
arbitrary size, by iteratively sliding their fixed extent over a large area of
an image (though these procedures are iterative and substantially slower than
likelihood computation/sampling that is less than or equal to the model patch
size). This is because both models can be autoregressively factored into a
product of conditional distributions. Using a 1-dimensional model with $D$
pixels for simplicity, this means that the joint distribution can be written
as:

$\displaystyle p(\mathbf{y})$ $\displaystyle=\prod_{k=1}^{D}p(y_{k}\mid
y_{1},y_{2},\dots,y_{k-1})$ (28) $\displaystyle=p(y_{1})p(y_{2}\mid
y_{1})p(y_{3}\mid y_{1},y_{2})\dots p(y_{D}\mid y_{1},y_{2},\dots,y_{D-1})$
(29)

Taking the log likelihood turns this product into a sum:

$\displaystyle\log p(\mathbf{y})$ $\displaystyle=\sum_{k=1}^{D}\log
p(y_{k}\mid y_{1},y_{2},\dots,y_{k-1})$ (30) $\displaystyle=\log p(y_{1})+\log
p(y_{2}\mid y_{1})+\log p(y_{3}\mid y_{1},y_{2})+\dots+\log p(y_{D}\mid
y_{1},y_{2},\dots,y_{D-1})$ (31)

Computing the likelihood over an image larger than the length $D$ model the
pixel was trained on can be accomplished by adding additional terms to this
sum. For example, computing the log likelihood of an $D+1$ length image would
require adding an additional term to the sum of the log likelihood of the the
$D+1$ pixel conditioned on the previous $D-1$ pixels (the maximum extent of
the model) that preceded it:

$\displaystyle\log p(\mathbf{y})$ $\displaystyle=\left(\sum_{k=1}^{D}\log
p(y_{k}\mid y_{1},y_{2},\dots,y_{k-1})\right)+\log p(y_{D+1}\mid
y_{2},y_{3},\dots,y_{D})$ (32)

Details of the likelihood computations for PixelCNNs are described in
[oordPixelRecurrentNeural2016]. For stationary Gaussian processes, there is a
closed form solution for finding the mean $\mu_{D}$ and variance
$\sigma^{2}_{D}$ of a 1-dimensional Gaussian distribution conditioned on $D-1$
previous values. This involves decomposing the covariance matrix into a top
left $(D-1)\times(D-1)$ block $\mathbf{\Sigma}_{1,1}$, a top right
$(D-1)\times 1$ column vector $\mathbf{\Sigma}_{1,2}$, and a bottom left
$1\times(D-1)$ row vector $\mathbf{\Sigma}_{2,1}$:

$\displaystyle\mathbf{\Sigma}=\begin{bmatrix}\mathbf{\Sigma}_{1,1}&\mathbf{\Sigma}_{1,2}\\\
\mathbf{\Sigma}_{2,1}&\sigma^{2}_{M}\\\ \end{bmatrix}$ (33)

Given $D-1$ previous values $\mathbf{y}_{D-1}$, the mean and variance of the
$D$th value can be computed as:

$\displaystyle\mu_{D}$
$\displaystyle=\mu+\mathbf{\Sigma}_{2,1}\mathbf{\Sigma}_{1,1}^{-1}(\mathbf{y}_{D-1}-\mu)$
(34) $\displaystyle\sigma^{2}_{D}$
$\displaystyle=\sigma^{2}-\mathbf{\Sigma}_{2,1}\mathbf{\Sigma}_{1,1}^{-1}\mathbf{\Sigma}_{1,2}$
(35)

The likelihood the the $D$th pixel can then be evaluated using the probability
density function of a 1D Gaussian distribution with mean $\mu_{D}$ and
variance $\sigma^{2}_{D}$.

Similarly, sampling images larger than the patch size on which the models were
trained can be accomplished by iteratively sampling each pixel conditioned on
the previous $D-1$ pixels.

###  1.6 Accounting for encoder uncertainty

As described in Section LABEL:sec:encoder_uncertainty_and_calibration,
measurements often contain information about both the object and the system
itself. To avoid overestimating object information, the strategy for
estimating mutual information must be modified.

Specifically, instead of estimating the full entropy of noisy measurements
$H(\mathbf{Y})$, the conditional entropy $H(\mathbf{Y}|\mathbf{E})$ should be
estimated. The conditional entropy estimate $H(\mathbf{Y}|\mathbf{X})$ is
unaffected. But $H(\mathbf{Y})$ now includes randomness from the object,
system variations, and noise. Conditioning on the encoder state $\mathbf{E}$
eliminates the contribution of system variations, leaving only object and
noise entropy.

In practice, conditioning can be accomplished by providing additional encoder
state information like field position or calibration images to the entropy
model during training. Flexible models like PixelCNN can readily accept such
conditional data [vandenoordConditionalImageGeneration2016]. This reduces the
entropy estimate by eliminating variations predictable from the system state.
We leave developing such conditional entropy estimation to future work.

##  2 Estimating information–Experiments

Having described the theoretical background of information estimation, we now
turn to experiments on simulated and experimental data to validate its
performance. For both methods of entropy estimation described below, we treat
noiseless images and noisy measurements as stationary stochastic processes, as
described in Section 1.5.1. The training data used to optimize the models
described below are thus cropped patches from random locations in the images,
since patches from any location should be statistically identical under such a
model. When simulating data, we can ensure that images are in fact stationary
by generating them from shift invariant probability distributions.
Experimental data does not have this property in reality, so the estimates
produced should be regarded as an average over spatially varying information
rates.

###  2.1 Probabilistic models for entropy estimation

####  2.1.1 Fitting stationary Gaussian processes

Fitting a stationary Gaussian process to a dataset of $N$
$\sqrt{D}\times\sqrt{D}$ image patches consists of using the image patches to
estimate the mean vector $\boldsymbol{\mu}$ and covariance matrix
$\mathbf{\Sigma}$ of the model. The mean vector can be easily estimated by
taking the mean of all pixels in the dataset across all patches.

The covariance matrix of a stationary Gaussian process has a special
structure: it is symmetric and positive definite (like every covariance
matrix), as well as being doubly Toeplitz (Section 1.5.2). This means that it
is composed of blocks of submatrices with constant diagonals, and the blocks
themselves are arranged in a repeated pattern along the diagonals.

Our initial estimator for the covariance matrix consisted of taking the set of
image patches, vectorizing them, computing the covariance matrix of the
sample, and then averaging along the diagonals and blocks to enforce this
special structure. However, this procedure is not guaranteed to produce a
covariance matrix that is positive definite, and we found empirically that the
averaging operation led the resulting estimate to have negative eigenvalues,
making an invalid covariance matrix.

To address this, we then regularized the estimates by enforcing an eigenvalue
floor, a minimum small, positive value for the eigenvalues. This was
insufficient for two reasons. 1) Modifying the eigenvalues in practice tended
to then break the doubly Toeplitz structure of the matrix, making it slightly
violate the requirements of a stationary Gaussian process. 2) More
importantly, the purpose of fitting this model was to be able to evaluate the
likelihood of data under the model to upper bound entropy, and this likelihood
is highly sensitive to the small eigenvalues of the covariance matrix. Thus,
the heuristic choice of eigenvalue floor has a significant impact on the
estimated entropy upper bound, and that effect varied from dataset to dataset.

To address these limitations, we developed an iterative optimization procedure
to improve the fit of the covariance matrix to the data. Using a dataset of
$N$ $\sqrt{D}\times\sqrt{D}$ image patches vectorized into $N$ $D$ length
vectors
$\mathbf{y}^{\left(1\right)},\mathbf{y}^{\left(2\right)},\dots,\mathbf{y}^{\left(N\right)}$,
we solve the following optimization problem to minimize the negative log
likelihood of a multivariate Gaussian distribution
$\mathcal{N}(\cdot;\boldsymbol{\hat{\mu}},\mathbf{\Sigma})$:

$\displaystyle\hat{\mathbf{\Sigma}}=\operatorname*{arg\,min}_{\mathbf{\Sigma}}-\sum_{i=1}^{N}\log\mathcal{N}(\mathbf{y}^{\left(i\right)};\boldsymbol{\hat{\mu}},\mathbf{\Sigma})$
(36)

Where $\boldsymbol{\hat{\mu}}$ is the fixed mean vector estimated from the
data as described above.

We solved this optimization problem using proximal gradient descent with
momentum. At each step, we computed the eigendecomposition of the covariance
matrix, fixed the eigenvectors, and computed and applied the gradient to the
eigenvalues. After each gradient step, a proximal operator was applied that
averaged along diagonals and blocks as described above, and then enforced a
minimum eigenvalue floor. The procedure was regularized using early stopping,
in which the optimization was stopped when the loss had not decreased for a
fixed number of iterations, and the parameters from the iteration with the
lowest loss were used as the final estimate. The parameters were initialized
with the estimate directly computed from the data as described above.

Empirically, we found that the optimization procedure was unstable and prone
to divergence. This seemed to be because the likelihood was very sensitive to
small eigenvalues of the covariance, which tended to produce extremely large
gradients for these eigenvalues. To account for this, we implemented gradient
clipping to limit the magnitude of the gradients for the eigenvalues.

The likelihood of a probability distribution can be very sensitive to small
changes in parameters because the Euclidean distance between parameter vectors
does not always accurately reflect the dissimilarity of the resulting
distributions (as described in section 2.3 of
[hoffmanStochasticVariationalInference2013]). In theory this can be corrected
by multiplying the gradient by the inverse of the Fisher information matrix to
reorient it from the steepest direction in Euclidean space to the steepest
direction in Riemannian space, which accounts for the natural geometry of the
parameters. However, in practice we found this to be expensive to compute
without obvious performance benefits over gradient clipping.

Figure 2 shows the results when fitting stationary Gaussian processes to
images from the BSCCM dataset [pinkardBerkeleySingleCell2024]. The direct
estimate contained several negative eigenvalues, which were then set to an
arbitrary eigenvalue floor. The optimization procedure was able to correct
these negative eigenvalues, and produced robust results for a variety of
settings of the eigenvalue floor parameter (provided it was smaller than the
true minimum eigenvalue) (Fig. 2b). As a consequence of the small, incorrect
eigenvalues in the direct estimate, samples produced from the model that were
larger than the patch size on which they were trained (Sec. 1.5.4) exhibited
numerical instability that created large oscillation in the samples. In
contrast, samples from the optimized model were stable and did not exhibit
this oscillation (Fig. 2b). The optimization quickly converged in a few
iterations (Fig. 2c), taking $\sim$ 3 seconds to complete on an NVIDIA GeForce
RTX 3090 GPU. Comparing estimates of a non-stationary Gaussian process, the
direct estimate stationary Gaussian process, and the optimized stationary
Gaussian process showed that the optimized stationary Gaussian process
produced more robust results in terms of its eigenvalue distribution using
smaller datasets (Fig. 2d).

Figure 2: Fitting Gaussian Process models. a) Eigenvalues of the covariance
matrix before and after optimization, indicating improved conditioning post-
optimization. The conditioning of the covariance matrix of non-optimized fits
is highly dependent on the choice of the eigenvalue floor hyperparameter. b)
Comparison of numerical stability during sampling pre- and post-optimization,
demonstrating that optimization prevents the instability that manifests in the
initial fit. c) The optimization process takes only a few iterations to
converge and shows an extremely large improvement in the negative log
likelihood of the data for poor initial fits. d) Efficacy of the optimization
across varying training set sizes, illustrating that optimized fitting
requires fewer images to achieve a comparable or better fit than the non-
optimized approach.

####  2.1.2 Fitting PixelCNNs

PixelCNNs were constructed and optimized as described previously
[oordPixelRecurrentNeural2016, vandenoordConditionalImageGeneration2016].
Briefly, the architecture consists of a series of masked convolutions, to
ensure an autoregressive ordering to pixels, with a final output layer that
outputs a conditional probability distribution at each pixel. This
distributions can be used both to iteratively generate samples, and to
evaluate the likelihood of data under the model.

We made one important modification to the standard PixelCNN architecture.
Rather than using a softmax layer as the final output for the conditional
distribution at each pixel, we use a mixture of Gaussians parameterized by the
output of the network (also known as a mixture density network
[bishopMixtureDensityNetworks1994]). This strategy is capable in theory of
approximating any conditional probability distribution, in the same way that a
neural network is capable of approximating any function
[hornikMultilayerFeedforwardNetworks1989]. This modification is essential for
estimating entropy in our framework, because it means the output distributions
will be a continuous probability density, instead of discrete probability mass
functions. In order to subtract the conditional entropy
$H(\mathbf{Y}\mid\mathbf{X})$ from the entropy estimate of $H(\mathbf{Y})$
made by the PixelCNN, it is essential that these are either both discrete
entropies (i.e. based on a probability mass function) or continuous entropies
(based on a probability density function). Since the optical noise models
$p(\mathbf{y}\mid\mathbf{x})$ we employ are continuous (e.g. Gaussian
distributions), the entropy estimates for $H(\mathbf{Y})$ must also be
continuous.

###  2.2 Comparing Stationary Gaussian Process and PixelCNN estimates

The fit of PixelCNN and Gaussian process models was compared by evaluating the
negative log likelihood on held out test data, and by visualizing samples from
the models. An important consideration for fitting either model to data is the
size of image patches used in the training set. Larger patches allow the model
to capture longer range dependencies and thus fit the data better in terms of
per pixel negative log likelihood. However, this comes at the cost of greater
data requirements and computational resources to train. Thus we sought to find
the largest patch size that produced a sufficient fit to the data to form an
accurate estimate of entropy.

Figure 3 shows samples produced when fitting models to patches of different
sizes. As expected, larger patch sizes were able to capture more complex
statistical relationships between pixels (Fig. 3a), decrease the per pixel
loss (Fig. 3b), and more tightly upper bound the entropy of the noisy images
(Fig. 3c). The gains in performance were minimal beyond a patch size of
$35\times 35$ pixels for the Gaussian model, while the PixelCNN was able to
continue to produce a more accurate estimate as patch size increased further.
This is presumably due to its much greater flexibility in modeling complex
dependencies between pixels. However, the magnitude of the gains was small
relative to the differences between the estimates of mutual information for
different contrast modalities (Fig. 4) for these experimental images.

Figure 3: The effects of patch size on model fit. a) Samples generated from a
Stationary Gaussian Process and PixelCNN, alongside samples from the the true
distribution. As the patch size increases, the models are able to capture more
long range dependencies between pixels, and b) achieve lower negative log
likelihood per pixel and thus c) a tighter upper bound on mutual information.
Figure 4: Mutual information and photon count. As the average number of
photons per pixel increases, the signal-to-noise ratio and mutual information
both increase. (Top) Example images of three different contrast modalities at
varying photon counts. (Bottom) The mutual information per pixel for each
contrast modality as a function of photon count.

When comparing the samples produced by the two models when fit to images with
different illumination patterns on the LED array, both the Gaussian process
and PixelCNN models were able to capture statistical patterns to each type of
contrast modality (Fig. 5a). Samples from the trained models appear to show
that the Gaussian model was able to learn the statistical patterns of the
texture of the images of objects under different illumination, while the
PixelCNN model additionally learned higher order structures like the edges of
cells.

Figure 5: Model samples for different contrast modalities. Samples from the
stationary Gaussian process and PixelCNN models for Brightfield, Differential
phase contrast, and Single LED illumination.

###  2.3 Failures of stationary Gaussian estimates on highly non-Gaussian data

The similar performance between the stationary Gaussian process and PixelCNN
models on the BSCCM dataset of single cells taken with coded-angle
illumination is likely due to the fact that the data is relatively close to a
multivariate Gaussian distribution, which is not something that can be assumed
to be true for all types of data. An extreme example of this can be seen by
comparing mutual information estimates using the MNIST handwritten digits
dataset as the object [lecunGradientbasedLearningApplied1998].

To demonstrate this, we take the MNIST dataset and simulate a minimal optical
encoder (i.e. a single lens) by convolving the images with a Gaussian point
spread function. We then add Poisson noise to simulate noisy measurements, and
fit both the stationary Gaussian process and PixelCNN models to the data.
Sampling from the fitted models shows that the PixelCNN model can produce
images that appear similar to true measurements, unlike the Gaussian process
model (Fig. 6a).

This occurs because the best approximating Gaussian fit is very dissimilar to
the true distribution of the measurements. This can be seen by looking at the
histograms of image pixels of the true data and samples from the models,
pooled across all many images (Fig. 6b). The true data and PixelCNN samples
have similar histograms, while the histogram from the Gaussian process samples
has a very different shape. As a result of this poor fit, the upper bound
given by the mutual information estimator is quite loose for the stationary
Gaussian model compared to the PixelCNN model (Fig. 6c).

These findings demonstrate that a stationary Gaussian process model is
inadequate to capture certain image distributions. Since the tightness of the
upper bound in our information estimation procedure depends on accurately
fitting the distribution of noisy measurements, this suggests that achieving a
tight upper bound/accurate estimate in many cases may require more flexible
probabilistic models, like the PixelCNN.

Figure 6: PixelCNN can fit non-Gaussian data better than stationary Gaussian
processes. a) Samples from the stationary Gaussian process and PixelCNN models
for MNIST digits. The Gaussian process samples do not resemble the true
distribution, indicating a poor fit. b) Histograms of pixel values for the
true distributions and samples from the two models indicate that the marginal
distribution of pixels is non-Gaussian, which the PixelCNN, unlike the
Gaussian model, is able to fit. c) As a result of the poor fit, the upper
bound on mutual information given by the Gaussian process estimator is much
looser than the PixelCNN estimator.

###  2.4 Consistency of mutual information estimates

An important and desirable property of an estimator is its consistency:
whether it converges to the true value of the parameter being estimated given
enough data. This is also an important practical consideration for using
estimators, because it determines how much data is needed to achieve a desired
level of accuracy.

In order to evaluate the accuracy of the estimators on high-dimensional data,
we must first generate samples from an image distribution with known mutual
information. This can be accomplished using a multivariate Gaussian
distribution with independent additive Gaussian noise at each pixel 13. In
this scenario, $H(\mathbf{Y}\mid\mathbf{X})$ is simply a constant, and
$H(\mathbf{Y})$ can be computed analytically by taking the covariance matrix
of the noiseless data and adding the variance of the noise to its diagonal,
and analytically computing the entropy of the resulting multivariate Gaussian
distribution.

We started by estimating a stationary Gaussian process from real data filtered
with a 3x3 median filter to use as the ground truth distribution of noiseless
images. We then estimated the mutual information of noisy samples from this
distribution using the PixelCNN, stationary Gaussian estimator, and a non-
stationary multi-variate Gaussian estimate (found by directly using the
covariance matrix of the samples) (Fig. 7a). All three estimators were able to
accurately estimate the mutual information of the data given enough samples.
The stationary Gaussian process estimator was the most efficient at this,
showing the lowest variance in its estimate of mutual information for a given
number of samples. The full Gaussian process estimator was the least
efficient, due to its sensitivity to the eigenvalue floor hyperparameter
described in (36).

These results are unsurprising, given that the stationary Gaussian process is
the least flexible class of models and the data itself was generated from a
stationary Gaussian process. To evaluate a more realistic scenario where data
is not perfectly Gaussian, we repeated the same experiment on cell images
under single LED illumination (Fig. 7b). In this scenario, the stationary
Gaussian estimator showed lower variance, more accurate estimates for smaller
amounts of data, but was less accurate than the PixelCNN estimator for larger
amounts of data. Despite not knowing the ground truth value of mutual
information in this scenario, since our estimator is an upper bound, we know
that lower estimates are more accurate. The full Gaussian estimator was the
least accurate in this scenario. Once again these results are unsurprising:
real data is not Gaussian, and flexible, neural network models like the
PixelCNN are able to produce more accurate results given enough data and
computational resources.

Figure 7: Estimator consistency. a) Estimating mutual information on simulated
data consisting of samples from a stationary Gaussian process. With an
increasing number of samples, all three estimators converge to the true value
of mutual information, with the stationary Gaussian estimator producing the
most reliable estimates for a given amount of data. b) When estimating on non-
Gaussian data, the stationary Gaussian process has less variation in its
estimates for a given amount of data, but cannot perfectly fit the non-
Gaussian data, as seen from the fact that it converges to a looser upper bound
than the PixelCNN estimator.

###  2.5 Conditional entropy estimates on noisy data

The accuracy of the mutual information estimate depends on the accuracy of the
separate estimates of its two components: $H(\mathbf{Y})$ and
$H(\mathbf{Y}\mid\mathbf{X})$. To better understand which estimate limits the
accuracy of the mutual information estimator, we examined the consistency of
the conditional entropy estimator $H(\mathbf{Y}\mid\mathbf{X})$ on simulated
data. When the noise model is additive Gaussian noise, the conditional entropy
is independent of the noiseless images, and thus can be computed analytically
exactly (Sec. 1.3.1). We thus focus on the more challenging scenario of
signal-dependent Poisson noise (Sec. 1.3.2).

Once again using images of cells from the BSCCM dataset, we filtered the
images with a 3$\times$3 median filter to produce simulated noiseless images
and then added simulated shot noise to the noiseless images to produce noisy
images for which the ground truth conditional entropy could be computed
analytically. We then compared our estimates of the conditional entropy to
this ground truth value as a function of the number of images used to estimate
it (Fig. 8a). Our conditional entropy rapidly converged to the true value,
with variations in its estimate on the order of $10^{-2}$ differential entropy
per pixel, even when the estimate was made from a single image patch. These
variations were $1-2$ orders of magnitude smaller than the variations seen in
the estimation of $H(\mathbf{Y})$ and the differences between different
encoders, leading us to conclude that the conditional entropy estimation was
an insignificant source of error in our information estimates.

Our approach for estimating conditional entropy relies on a Monte Carlo
approximation using samples from the distribution of clean images (Equation
5). However, when applying the mutual information estimator to experimental
data, noiseless images are not available. Thus we must resort to an
approximation in which we estimate this quantity from samples from the
distribution of noisy measurements instead. This may introduce some error into
the estimator, and it is important to understand how much and what effect this
error will have on the final estimate of mutual information.

In the case of additive Gaussian noise (Sec. 1.3.1), the conditional entropy
is independent of the noiseless images, and thus there is no additional error.
However, in the case of shot noise (Sec. 14), the conditional entropy is
dependent on the noiseless images, and thus the approximation will introduce
some error. Mathematically, our approach to estimating conditional entropy is
to replace the noiseless pixel values in Equation 17 with noisy pixel values:

$\displaystyle\mathrm{H}\left(\mathbf{Y}\mid\mathbf{X}\right)$
$\displaystyle\approx\frac{1}{N}\sum_{i=1}^{N}\sum_{k=1}^{D}\frac{1}{2}\log_{2}(2\pi
ex_{k}^{\left(i\right)})$ (37)
$\displaystyle\approx\frac{1}{N}\sum_{i=1}^{N}\sum_{k=1}^{D}\frac{1}{2}\log_{2}(2\pi
ey_{k}^{\left(i\right)})$ (38)

To test this we generated simulated noiseless images by applying a $3\times 3$
median filter to images from the BSCCM dataset and generated simulated noisy
images by adding Poisson noise to these. We then estimated mutual information
for different contrast modalities and photon counts, either using the true
noiseless images or the simulated noisy images to estimate the conditional
entropy. We found that the error introduced by using the noisy images was
small relative to the differences between the estimates of mutual information
for different contrast modalities (Fig. 8b), with appreciable differences only
occurring at photon counts below $\sim 40$ photons per pixel. Furthermore, the
errors introduced a consistent bias in the estimates, giving overestimates of
mutual information due to underestimates of the conditional entropy. Given the
predictable nature of this bias, it could be possible to explicitly introduce
a correction factor to the estimator to account for it. However, given that
other assumptions of our estimator, like the Gaussian approximation to Poisson
noise (Sec. 1.3.2), are also not valid in this regime, we leave this for
future work and instead focus on the higher photon count regime where our
estimator is more accurate.

Figure 8: Estimating conditional entropy. a) The estimate of conditional
entropy for shot noise quickly converges to the true value as the number of
samples increases, showing much smaller variations than the entropy estimators
considered in Figure 7. b) When noisy experimental images are used to estimate
conditional entropy instead of ground truth noiseless images, the error
introduced is small relative to the differences between the estimates of
mutual information for different contrast modalities, and only deviates from
the true value at low photon counts.

###  2.6 Analytic Gaussian entropy vs. upper bound

In addition to estimating entropy by upper-bounding it using the test set
negative log likelihood, an alternative is to fit a probabilistic model to the
data, and then analytically compute the entropy of the model. Though this is
not possible for all models, it is possible for the stationary Gaussian
process model, since the entropy of a multivariate Gaussian distribution can
be computed analytically. It is unclear what the theoretical relationship
between this analytic estimate and the true entropy is, unlike the test set
negative log likelihood method, which upper bounds the true entropy.

We tested this approach empirically and compared it to the upper bound entropy
estimates produced by the PixelCNN and stationary Gaussian process models on
three different contrast modalities of the BSCCM dataset (Fig. 9). In this
experiment, the values for the analytic stationary Gaussian entropy estimator
were slightly lower than those of the stationary Gaussian upper bound, while
being slightly higher than the PixelCNN upper bound estimates. Since we know
the true entropy is lower than the PixelCNN upper bound, this suggests that
the analytic Gaussian entropy is more accurate than the Gaussian upper bound
estimator. However, these differences were small in comparison to the
differences between the estimates of mutual information for different contrast
modalities and different photon counts. The similarity between the analytic
Gaussian entropy and the upper bound Gaussian entropy is advantageous, as we
employ the analytic Gaussian entropy in our IDEAL (Sec.
LABEL:sec:main_paper_IDEAL) loss function, and its close agreement with the
upper bound estimator supports its validity.

Figure 9: Analytic Gaussian entropy vs. upper bound. The analytic Gaussian
entropy estimator produces slightly different estimates than the upper bound
estimator for the stationary Gaussian process, but these differences are small
relative to the differences between the estimates of mutual information for
different contrast modalities and different photon counts.