Before going into the details of what SLIC is all about, let’s first start by clearly defining the data generation process. In the typical case where the noise is additive, the observed data \(D\) can be expressed as

\begin{equation}
D = M(\theta) + N, \tag{1}
\end{equation}

where \(M\) represents the underlying signal and \(N\) signifies the additive noise. Here, the dependence of \(M(\theta)\) on \(\theta\) stems from the assumption that the signal \(M\) within the observed data is fully describable by a latent set of physical parameters \(\theta\). For the most part, the goal of data analysis is to infer the value of these latent parameters \(\theta\) from the data \(D\). To do this, we need to know the likelihood distribution, which characterizes the probability of observing the data \(D\) given the true signal \(M(\theta)\):
\begin{equation}
P(D \mid M(\theta)).
\label{likelihood}\tag{2}
\end{equation}

This term conveys that even if the underlying true signal \(M(\theta)\) is precisely defined, the data \(D\) remains probabilistic. But why is \(D\) probabilistic? It’s because the data is influenced by random noise, which cannot be precisely predicted or controlled.

In fact, Equation \ref{likelihood} implies that the noise \(N\) is equivalent to the residuals between the data \(D\) and the signal \(M(\theta)\):

\begin{equation}
N = D \, – M(\theta). \tag{3}
\end{equation}

Therefore, if we look at the probability of observing a specific noise realization \(P(N)\), we find that this is equivalent to the probability of the residuals \(P(D – M(\theta))\). Assuming \(M(\theta)\) is deterministic (meaning that \(M\) is always the same for a given value of \(\theta\)), this also means that the probability \(P(D – M(\theta))\) is equivalent to the probability of \(D\) for a given \(M(\theta)\), which is the definition of the likelihood \(P(D \mid M(\theta))\). Hence, this leads us to the conclusion that

\begin{equation}
P(N) = P(D \, – M(\theta)) = P(D \mid M(\theta)). \tag{4}
\end{equation}

In other words, when the noise is additive and \(M(\theta)\) is deterministic, the likelihood distribution IS the probability distribution of the noise \(P(N)\), which is consequently also the probability of the residuals \(P(D – M(\theta))\). Therefore, the task of figuring out the correct likelihood comes down to figuring out the correct noise distribution \(P(N)\).

When the Central Limit Theorem (CLT) holds, \(P(N)\) would be a Gaussian distribution. For example, if we have \(K\) dimensional uncorrelated Gaussian noise, the probability distribution \(P(N)\), or equivalently \(P(D – M(\theta))\), can be written as:

\begin{equation}
P(N) = \prod_{k=1}^K \frac{1}{\sqrt{2\pi\sigma_k^2}} \exp\left(-\frac{1}{2} \frac{(D_k – M_k(\theta))^2}{\sigma_k^2}\right). \tag{5}
\end{equation}

Alternatively, if we look at the logarithm of \(P(N)\), it can be written as

\begin{equation}
\log P(N) = – \frac{1}{2} \sum_{k=1}^K \left[\log(2\pi \sigma_k^2) + \frac{(D_k – M_k(\theta))^2}{\sigma_k^2} \right]. \tag{6}
\end{equation}

Ignoring the first term as it is a constant, the second term in \(\log P(N)\) specifically \(\sum_k^K \frac{(D_k – M_k(\theta))^2}{\sigma_k^2}\), is none other than the well-known chi-squared statistic \(\chi^2\), which we routinely use to measure how well our models of \(M(\theta)\) fit the observational data \(D\). This is to say that the reliance on the chi-squared statistic \(\chi^2\) to assess the goodness of fit inherently assumes Gaussian-distributed noise, a presumption rooted deeply in the Central Limit Theorem.

Fun facts aside, as a reminder, we are interested in the cases where the Central Limit Theorem does not hold and where the noise in the data can significantly deviate from a Gaussian distribution. Then the question is: how do we model \(P(N)\)? This is where machine learning comes to the rescue: if we have real samples of this complex noise, we can train a machine learning model on this dataset to learn \(P(N)\). This could be done, for example, by collecting samples from parts of the data that contain only noise, without any signal. Then, with a model for the noise distribution \(P(N)\), we can evaluate the probability of different noise realizations, or equivalently, evaluate the probability of the residuals \(D – M(\theta)\). The latter application is particularly useful when trying to find values of \(\theta\) in which the residual \(D – M(\theta)\) is consistent with noise, indicating an accurate model fit of the data.

Theory aside, in practice, how feasible is it to model \(P(N)\)? For instance, when dealing with noise in images from the Hubble Space Telescope, the dimensionality of \(N\) can be extremely large, approximately 16.8 million \((4096 \times 4096)\) pixels for a single full-resolution image. For these cases where the dimensionality of \(N\) is large, modeling \(P(N)\) suffers from the curse of dimensionality, in the sense that the normalization of a probability distribution becomes intractable in high dimensions. Consequently, only in the simplest cases (such as Gaussian noise) could we find a tractable model for \(P(N)\).

If \(P(N)\) is intractable, then what is the next best thing that we can learn about \(P(N)\)? Perhaps instead of \(P(N)\), we can learn the score of \(P(N)\), defined as \(\nabla_{N} \log P(N)\). In contrast, \(\nabla_{N} \log P(N)\) doesn’t require \(P(N)\) to be normalized and only requires local information to be determined.

The question now is: how do we go about learning \(\nabla_{N} \log P(N)\)? One possibility is to train score-based generative models to learn to predict \(\nabla_{N} \log P(N)\) from samples of noise-only data. We will skip through the details of learning \(\nabla_{N} \log P(N)\), but for more information you can take a look at some of our recent work that does exactly this (for example, here and here).

Fast-forward a couple days of training some score-generative models, assuming we’ve correctly learned \(\nabla_{N} \log P(N)\), what can we do with it? As a reminder, the goal of inference is to infer the values of the latent parameters \(\theta\) that describe the background signal given the data \(D\). One possibility is to use gradient-based Markov Chain Monte Carlo methods (MCMC) to sample the space of \(\theta\) values given \(D\). However, if we wish to explore the space of \(\theta\) values, we need gradients with respect to \(\theta\), not \(N\); that is, we need \(\nabla_{\theta} \log P(N) = \nabla_{\theta} \log P(D – M(\theta))\), not \(\nabla_{N} \log P(N)\).

Fortunately for us, there exists a powerful mathematical tool, often revered in ancient scripture: the chain rule of calculus. Using the chain rule, \(\nabla_{\theta} \log P(N)\) and \(\nabla_{N} \log P(N)\) are related via

\begin{equation}
\nabla_{\theta} \log P(N) = \nabla_{N} \log P(N) \cdot \nabla_{\theta} N. \label{slic_equation}\tag{7}
\end{equation}

Given the definition of the data generation process \(N = D – M(\theta)\), we have

\begin{equation}
\nabla_{\theta} \log P(N) = \nabla_{N} \log P(N) \cdot \nabla_{\theta} (D\, – M(\theta)) = -\nabla_{N} \log P(N) \cdot \nabla_{\theta} M(\theta).\tag{8}
\end{equation}

And there we have it! We can obtain the gradient of the likelihood with respect to \(\theta\) from the gradient of the noise distribution \(\nabla_{N} \log P(N)\). The only additional quantity needed is \(\nabla_{\theta} M(\theta)\), i.e., the Jacobian of the physical model. \(\nabla_{\theta} M(\theta)\) can easily be computed if we have access to a differentiable forward simulator for \(M(\theta)\), which these days is becoming increasingly common due to the growing use of auto-differentiable coding packages such as PyTorch and JAX.

To summarize, the framework we’ve introduced, referred to as Score-based Likelihood Characterization (SLIC), enables data analysis amid non-Gaussian noise. The key procedures involved in SLIC are:

1. Learning the Noise Score: We estimate the score of the noise distribution, denoted as \(\nabla_N \log P(N)\), using data samples that do not contain the signal of interest.
2. Calculating the Jacobian: We compute the Jacobian of the forward simulator model \(\nabla_{\theta} M(\theta)\).
3. Gradient Computation of Likelihood: Utilizing Equation \ref{slic_equation}, we calculate the gradient of the likelihood with respect to the parameters \(\theta\). This calculation forms the basis for employing gradient-based sampling algorithms to infer \(\theta\) from the observed data \(D\).

In a previous work of ours, we tested SLIC using data from both the Hubble Space Telescope and the James Webb Space Telescope (JWST). Specifically, we used low-level data that underwent minimal preprocessing such that significant non-Gaussian noise, including cosmic rays were present. After learning the score of this non-Gaussian noise distribution, we applied it to infer the true model parameters of simulated background signals from mock observations containing real HST and JWST noise. Using SLIC, we found that we were able to accurately perform inference of the underlying signal, while using a Gaussian likelihood clearly showed significant bias in the inferred model parameters. For example, in Figure 2 we show one such result, where the noise is taken from real HST data and the background signal is a simulation of a strong gravitational lensing system. As demonstrated in Figure 2, learning the correct noise distribution through SLIC allowed for unbiased parameter inference.

In conclusion, the Score-based Likelihood Characterization (SLIC) framework marks a significant leap forward in data analysis, particularly in tackling the complexities of non-Gaussian noise. By revolutionizing our approach to modeling noise distribution through score-based generative models, SLIC transcends the limitations of traditional Gaussian assumptions. This paradigm shift not only improves the accuracy of our inferences but consequently also paves the way for deeper exploration and understanding of the universe. The capability to precisely model complex, high-dimensional noise from sources like the Hubble Space Telescope enhances our ability to interpret astrophysical phenomena by extracting every bit of information from the data. As we continue to push the frontiers of observational astrophysics, frameworks like SLIC will be indispensable, driving new discoveries and fostering a profound connection between our data and the cosmos. The future of data analysis is bright, and with tools like SLIC, we are better equipped than ever to unravel the mysteries of the universe, one noisy dataset at a time.