11institutetext:  Both authors contributed equally to this work.
Faculty of Mathematics, University of Vienna, Oskar-Morgenstern-Platz 1,
A-1090 Vienna
Department of Astrophysics, University of Vienna, Türkenschanzstraße 17,
A-1180 Vienna, Austria
Max-Planck Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg,
Germany
Instituto de Astrofísica de Canarias, C/ Vía Láctea s/n, E-38205 La Laguna,
Spain
Facultad de Ingeniería y Arquitectura, Universidad Central de Chile, Av.
Francisco de Aguirre 0405, La Serena, Coquimbo, Chile
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218,
USA
Johann Radon Institute for Computational and Applied Mathematics (RICAM),
A-4040 Linz
Christian Doppler Laboratory for Mathematical Modeling and Simulation of Next
Generations of Ultrasound Devices (MaMSi), Oskar-Morgenstern-Platz 1, A-1090
Vienna
11email<EMAIL_ADDRESS>
# Uncertainty-Aware Blob Detection with an Application to Integrated-Light
Stellar Population Recoveries

Fabian Parzer **11 Prashin Jethwa **22 Alina Boecker 3344 Mayte Alfaro-Cuello
5566 Otmar Scherzer 117788
Glenn van de Ven 22**1122334455667788

###### Abstract

Context. Blob detection is a common problem in astronomy. One example is in
stellar population modelling, where the distribution of stellar ages and
metallicities in a galaxy is inferred from observations. In this context,
blobs may correspond to stars born in-situ versus those accreted from
satellites, and the task of blob detection is to disentangle these components.
A difficulty arises when the distributions come with significant
uncertainties, as is the case for stellar population recoveries inferred from
modelling spectra of unresolved stellar systems. There is currently no
satisfactory method for blob detection with uncertainties.

Aims. We introduce a method for uncertainty-aware blob detection developed in
the context of stellar population modelling of integrated-light spectra of
stellar systems.

Methods. We develop theory and computational tools for an uncertainty-aware
version of the classic Laplacian-of-Gaussians method for blob detection, which
we call ULoG. This identifies significant blobs considering a variety of
scales. As a prerequisite to apply ULoG to stellar population modelling, we
introduce a method for efficient computation of uncertainties for spectral
modelling. This method is based on the truncated Singular Value Decomposition
and Markov Chain Monte Carlo sampling (SVD-MCMC).

Results. We apply the methods to data of the star cluster M54. We show that
the SVD-MCMC inferences match those from standard MCMC, but are a factor 5-10
faster to compute. We apply ULoG to the inferred M54 age/metallicity
distributions, identifying between 2 or 3 significant, distinct populations
amongst its stars.

###### Key Words.:

[Methods: statistical] – [Galaxies: stellar content] – [Galaxies: star
clusters: individual: M54]

##  1 Introduction

Identifying distinct structures in a distribution is a common problem in
astronomy. When the distribution under consideration is an observed image,
this problem is called source detection (see Masias et al., 2012, for a
review). In many cases, however, the distribution is not an observed image,
but instead is itself the result of some earlier inference problem. Consider
the following examples:

  1. 1.

Gravitational lensing, where a mass distribution is inferred from an image
which displays lensing effects and structures such as dark matter halos are
identified in the mass distribution (e.g. Despali et al., 2022).

  2. 2.

Orbit-superposition modelling of galaxies, where an orbit distribution is
inferred from stellar kinematic data and structures, i.e. dynamical components
of the galaxy, are identified in the orbit distribution (e.g. Zhu et al.,
2022).

  3. 3.

Stellar population modelling of unresolved stellar systems, where a
distribution over stellar ages and metallicities is inferred from spectral
modelling and structures, i.e. distinct populations, are identified there
within (e.g. Boecker et al., 2020).

A theme connecting all three examples - and the important distinction with
respect to observed images - is that the distribution may come with
significant uncertainties as a result of the initial inference.

Let’s consider stellar population modelling in more detail. Non-parametric
modelling of integrated-light spectra is a flexible way to recover stellar
population distributions. In this approach, libraries of several hundred
template spectra of simple stellar populations (SSPs) are superposed to model
observed spectra (e.g. Cappellari & Emsellem, 2004; Cid Fernandes et al.,
2005; Ocvirk et al., 2006; Tojeiro et al., 2007; Koleva et al., 2008). This
approach offers the flexibility to reconstruct complex star formation
histories (SFHs) where metallicity and alpha abundances (Vazdekis et al.,
2015) can vary freely with stellar age. In cases where stellar systems host
several distinct populations, nonparametric models can recover this complexity
(Boecker et al., 2020). A difficulty with nonparametric models is that the
high-dimensionality has hindered robust uncertainty quantification. An
alternative is to use parametric models, where the SFH is characterised with
at most a few dozen parameters (e.g. Gladders et al., 2013; Simha et al.,
2014; Carnall et al., 2018). Whilst uncertainty quantification is easier here,
parametric models impose strong assumptions which, if invalid, may bias
inferences (Lower et al., 2020) and lead to underestimated uncertainties (Leja
et al., 2019).

In this work we describe a method for uncertainty quantification for
nonparametric stellar population modelling. Furthermore, we introduce a
general method for _blob detection_ which can identify distinct blobs in
reconstructed distributions with uncertainties. In the context of stellar
populations, distinct blobs will correspond to groups of stars with different
origins – e.g. those formed in-situ within a galaxy versus those accreted from
merged, destroyed satellites (White & Rees, 1978). Due to the extended star-
formation histories of galaxies and their satellites, blobs may partially
overlap in age/metallicity space, therefore blob detection in this context is
distinct from the task of image segmentation.

Blob detection has traditionally been studied in the case where the underlying
image is known. When the image of interest has to be reconstructed from noisy
measurements, the result of the blob detection itself becomes uncertain.
Currently, the dominating approach for uncertainty quantification in imaging
is Bayesian statistics (Hanson, 1993; Kaipio & Somersalo, 2007; Dashti &
Stuart, 2017). While current research mostly focuses on the development of
computationally efficient methods for generic uncertainty quantification, in
particular MCMC (Bardsley, 2012; Martin et al., 2012; Parno et al., 2021), the
specific task of combining uncertainty quantification with blob detection has
not been explicitly studied. Repetti et al. (2018) proposed a method for
significance tests of image structures called ”Bayesian uncertainty
quantification by convex optimisation” (BUQO). Price et al. (2020) proposed a
similar test that detects a structure in an uncertain image by checking
whether a suitable surrogate image is contained in the posterior credible
region. See also the related works by Cai et al. (2018a 2018b) on uncertainty
quantification for radio interferometric imaging, and by Price et al. (2020,
2021) on sparse Bayesian mass mapping. The cited methods are computationally
very efficient, but do not perform automatic scale-selection, assuming instead
that the structure of interest is well-localised in a user-specified image
region. This makes them also less suitable for the case where structures might
overlap, as is the case with stellar populations.

We have based our approach for blob detection on the classic Laplacian-of-
Gaussians method (Lindeberg, 1998). This method is able to determine the
location and scale of blob-like structures in an arbitrary image. Although the
underlying ideas of the Laplacian-of-Gaussians have been modified and extended
in subsequent work (see e.g. the reviews by Mikolajczyk et al., 2005;
Lindeberg, 2013; Li et al., 2015), we have based our approach on the original
formulation due to its conceptual simplicity and straightforward mathematical
definition.

###  1.1 Overview

An overview of our technique is shown in figure 1. Given an observed spectrum
and a choice of prior, a set of reference blobs is determined from the
maximum-a-posteriori (MAP) estimate. Using posterior sampling we then
determine those components of the image that are distinguishable with
confidence. These _significant blobs_ are then combined with the reference
blobs using a matching procedure, yielding a final visualisation that shows
which blobs are significant, whether blobs are separable with confidence, and
which blobs are explainable by noise alone. An example of this is shown in
figure 3.

Figure 1: High-level overview of our method for uncertainty-aware blob
detection.

In more detail: First, we obtain the MAP estimate through optimisation (see
equation (4) and the surrounding discussion). We then apply the Laplacian-of-
Gaussians (LoG) method (described in section 3) to determine the position and
scale of blobs in the MAP image. An example of this is visualised in figure 2.

When accounting for uncertainties, however, not all of these blobs are
significant. To quantify the uncertainty in our estimate, we generate
posterior samples using MCMC. Standard MCMC algorithms require several hours
for our problem, so we developed a new technique for efficient sampling based
on dimensionality reduction, called SVD-MCMC (described in section 4).

To assess the significance of the image blobs from the posterior samples we
use a tube-based approach: First we compute so-called filtered credible
intervals (see sections 5.2 and 5.3) that, for a given scale, consist of an
upper and lower bound in which the true image lies with a high probability. We
can then determine significant structures by ”spanning a blanket”, i.e.
computing the image with the least amount of blobs that still lies within the
credible intervals (see section 5.4). Since the uncertainty of fine image
features is usually different than the uncertainty of large-scale structures,
this has to be repeated for multiple scales. Applying the LoG method to the
resulting stack of blankets then allows us to determine which image blobs are
significant. Hence, we call this technique ”Uncertainty-aware Laplacian-of-
Gaussians” (ULoG).

To visualise the results of ULoG as in figure 3, we use the MAP image as a
reference and highlight which blobs detected by the LoG method - visualised in
figure 2 \- are significant and which are not. The details of this procedure
are described in sections 5.4.1 and 5.4.2.

Figure 2: Results of the LoG method (see section 3) applied to the MAP
estimate of the stellar distribution function. The detected blobs are
indicated by red circles. Figure 3: Results of the proposed ULoG method to
detect blobs with 95% credibility. The mock data and MAP estimate are the same
as in figure 2. The significant blobs are visualised by dashed blue circles,
their size corresponding to localisation uncertainty. Blobs in the MAP
estimate that correspond to significant structures are marked in green, blobs
that might be explainable by noise are marked in red (see section 5.4.2).

##  2 Problem setup

Our goal is to robustly identify distinct blobs in stellar population
distribution functions by modelling integrated-light spectra. The distribution
function $f$ is a function of stellar metallicity $z$ and age $t$, and it is
physically constrained to be non-negative, i.e.

$f(z,t)\geq 0.$

This can be related to an observed spectrum via simple stellar population
(SSP) templates $S(\lambda,z,t)$ which encode the spectra of single-age and
single-metallicity stellar populations. The spectrum of a composite stellar
population is a combination of SSPs weighted by $f$. Assuming that all stars
in this composite population share a line-of-sight velocity distribution
(LOSVD) $\mathcal{L}(v)$, then the observed spectrum is given by the composite
rest-wavelength spectrum convolved by $\mathcal{L}(v)$. Putting this all
together, the observation operator $\mathcal{G}$ is given by

$\mathcal{G}(f,\mathcal{L})=\int\frac{\mathcal{L}(v)}{1+v/c}\left[\int\int
S\left(\frac{\lambda}{1+v/c},z,t\right)f(z,t)\;\mathrm{d}t\;\mathrm{d}z\right]\;\mathrm{d}v.$
(1)

The observed data $y$ is given by the addition of noise $\epsilon$ which is
wavelength-dependent,

$y(\lambda)=\mathcal{G}(f,\mathcal{L})+\epsilon(\lambda).$ (2)

Since the goal of this work is to develop methods to robustly identify
structures in $f$, we will fix the LOSVD $\mathcal{L}$ in our recoveries.
Whilst equations (1) and (2) relate continuous functions $y,f$ and
$\mathcal{L}$, in reality observations are made at a finite number of
wavelengths $\lambda$. We therefore instead work with the discretised version
of equation (2), for fixed $\mathcal{L}$, i.e.

$\bm{y}=\bm{G}\bm{f}+\bm{\epsilon}.$ (3)

Details of converting between continuous and discretised equations are in
section 2.1. The noise is typically assumed to be distributed as a
multivariate normal with some covariance ${\bm{\Sigma}^{\bm{y}}}$,

$\bm{\epsilon}\sim\mathcal{N}(0,{\bm{\Sigma}^{\bm{y}}}),$

where we assume that an estimate of ${\bm{\Sigma}^{\bm{y}}}$ is available as
is typically output from data-reduction pipelines. This completes the
definition of our forward model. For the time being, we ignore several other
physical ingredients which go into complete spectral models – e.g. LOSVD
modelling, gas emission, reddening and calibration related continuum
distortions. We discuss these effects further in Section 7.2.

Non-parametric recovery of $f$ from noisy observations $y$ is an ill-
conditioned inverse problem. This means that noise in the observations is
amplified in the recovery, leading to unrealistically spiky solutions.
Solutions are found by minimising the mismatch between data and model as
quantified by $\chi^{2}$,

$\chi^{2}=\left\|\bm{y}-\bm{G}\bm{f}\right\|_{\bm{\Sigma}^{\bm{y}}}^{2}.$

The _maximum likelihood_ (ML) solution is the one which minimises the
$\chi^{2}$ mismatch, i.e.

$\bm{f}_{\mathrm{ML}}=\operatorname*{argmin}_{\bm{f}}\chi^{2}.$

To counter ill-conditioning, regularisation can be used to promote smoothness
in the solution. In practice this means that rather than finding the ML
solution, one instead finds the _maximum a posteriori_ (MAP) solution by
adding a penalty term to the objective function,

$\bm{f}_{\mathrm{MAP}}=\operatorname*{argmin}_{\bm{f}}\left(\chi^{2}+\beta\mathcal{R}\right).$
(4)

The penalty consists of a functional $\mathcal{R}(\bm{f})$ which quantifies
the smoothness of $\bm{f}$ multiplied by a tunable scalar $\beta$ called the
regularisation parameter. The functional is non-negative, and decreases for
increasingly smooth $\bm{f}$. When $\beta=0$, the MAP and the ML solutions
coincide, while in the limit of very large $\beta$, $\bm{f}_{\mathrm{MAP}}$
will be very smooth at the expense of no longer well reproducing the data. By
appropriately tuning $\beta$ to some intermediate value however the MAP
solution will exhibit a good data-fit and a desired degree of smoothness. MAP
recoveries are commonly used in stellar population modelling: Ocvirk et al.
(2006) give a detailed description of the theory behind regularisation in this
context, and it is available in the popular spectral modelling software
Penalised Pixel-Fitting (pPXF, Cappellari, 2017).

We will not address the question of how to tune the regularisation parameter
$\beta$ in this work. Several other suggestions have been made elsewhere. For
example, the pPXF documentation 111https://pypi.org/project/ppxf/ offers the
guideline to first perform an un-regularised fit, then increase $\beta$ till
the data mismatch increases to $\Delta\chi^{2}=\sqrt{2n}$ where $n$ is the
size of the spectrum. Further well-known parameter choice rules are the
L-curve method (Hansen & O’Leary, 1993) and Morozov’s discrepancy principle
(Groetsch, 1983, provides a theoretical discussion). Alternative Bayesian
approaches for choosing the regularisation parameter are discussed in section
7.1. While the choice of regularisation parameter can affect the number of
blobs visible in a recovered $\bm{f}_{\mathrm{MAP}}$, in this work we restrict
ourselves to the question of quantifying the uncertainty in the solution at
fixed $\beta$ and then using these uncertainty estimates to identify distinct
blobs.

Regularisation can also be understood from a Bayesian perspective. Under this
interpretation, equation (4) comes from the combination of a likelihood
function

$p(\bm{y}|\bm{f})\propto e^{-\frac{1}{2}\chi^{2}}$

with a prior probability density function

$p(\bm{f})\propto e^{-\frac{1}{2}\beta\mathcal{R}}.$ (5)

Bayes’ theorem tells us that the posterior probability is proportional to the
product of the likelihood and prior, i.e.

$p(\bm{f}|\bm{y})\propto
e^{-\frac{1}{2}\left(\chi^{2}+\beta\mathcal{R}\right)},$

from which we can see that equation (4) picks out the solution
$\bm{f}_{\mathrm{MAP}}$ which maximises the (log) posterior probability
function. This is the origin of the name _maximum a posteriori_.

In this work, we assume a Gaussian prior on $\bm{f}$ that is truncated to the
nonnegative orthant, since the density $\bm{f}$ cannot take negative values.
That is, we make the choice

$p(\bm{f})\propto\cases{e}^{-\frac{\beta}{2}\left\|\bm{f}\right\|_{\bm{\Sigma}^{\bm{f}}}^{2}},&\text{if
}\bm{f}\geq 0,\\\ 0,\text{otherwise}.{}$ (6)

In the spirit of regularisation theory, we use the prior covariance
${\bm{\Sigma}^{\bm{f}}}$ to encode smoothness assumptions on $\bm{f}$
(Calvetti & Somersalo, 2018). Details on our choice for
${\bm{\Sigma}^{\bm{f}}}$ are given in section 2.1.

Viewed from this Bayesian perspective, we can explicitly see that
regularisation incorporates prior information into the solution, and the
regularisation parameter $\beta$ controls the strength of the prior. Whilst a
MAP solution may be considered Bayesian in the sense that it assumes a prior,
it is nevertheless a single point estimate which on its own provides no
information about uncertainty of the solution. In this work, we will introduce
an efficient method to estimate this uncertainty based on dimensionality-
reduced MCMC sampling.

The use of regularisation complicates another question commonly faced in
stellar population modelling: the choice of whether to reconstruct _mass-
weighted_ or _light-weighted_ quantities. These are two options for the
parametrisation of $\bm{f}$ connected to the normalisation of SSP templates
$S$ in equation (1). If the SSPs are normalised to represent the flux of a
unit _mass_ of stars, $\bm{f}$ is said to be mass-weighted. In this case, the
integrated flux of a stellar population is given by

$\text{Flux}(z,t)=\int S(\lambda,z,t)\;\mathrm{d}\lambda.$

SSP integrated flux varies strongly with age and, to a lesser extent,
metallicity (Figure 4). As an alternative to mass-weighting, the SSP templates
$S$ in equation (1) may be normalised to all have the same integrated flux, in
which case the recovery is said to be light-weighted. Maximum-likelihood mass-
weighted and light-weighted solutions should be self-consistent in the sense
that

$\textnormal{light-
weighted}\;\bm{f}_{\mathrm{ML}}(z,t)=\text{Flux}(z,t)\times\textnormal{mass-
weighted}\;\bm{f}_{\mathrm{ML}}(z,t).$

This consistency does not hold for MAP solutions. That is to say, the decision
to regularise mass-weighted or light-weighted $f$ will affect the result.

Figure 4: Examples of SSP templates. Left panel shows 9 representative models
from the MILES SSP library. These are labelled 1-9, corresponding to the
position on the age-metallicity plane as shown in the right panel. The
colorbar shows the integrated flux of the full grid of all 636 SSP templates
as a function of age and metallicity. Younger stars are hotter and therefore
brighter, explaining why integrated flux is a strong function of SSP age.
Integrated flux depends more weakly on metallicity. For fixed age, we see that
SSPs with greater metallicity have stronger absorption features e.g. compare
spectra 7, 8 and 9 in the left panel.

Is there a preferred option between regularised mass- or light-weighted
recoveries? We think there is no single objective correct answer to this
question. Physically, one may imagine a mass of gas converting to stars at a
rate varying smoothly with time, whence regularised mass-weights might seem to
be the preferred option. Given that integrated SSP flux varies smoothly with
age and metallicity however (see Figure 4), mass-weighted smoothness would
also imply light-weighted smoothness. Whilst both options may be reasonably
justified under physical considerations, we raise this discussion here since
_numerically_ the situation is more complicated. For the MCMC based
uncertainty quantification method we will present in Section 4, we find that
sampling a model with regularised mass-weights is numerically unstable; see
that section for further details. In order to facilitate comparison between
all methods presented in this work, we will therefore only consider
regularised light-weighted recoveries throughout. We also remark that this
choice can be interpreted as a ”linear functional strategy” (see Anderssen,
1986) for reducing the ill-posedness of the reconstruction.

###  2.1 Discretisation

Here we outline the details of converting from continuous variables $f$, $y$,
$S$ etc. in equation (2), to the discretised version for a fixed LOSVD, i.e.
equation (3). We assume that data is available at $n$ wavelengths

$\bm{y}=[y(\lambda_{1}),...,y(\lambda_{n})],$

where typical values of $n$ are around 1000. The domains of $t$ and $z$ can be
discretised into $T$ and $Z$ intervals respectively. In pratice, the age and
metallicity discretisations are set by the choice of SSP library. The MILES
library which we use has $(Z,T)=(12,53)$. We next introduce vectors
$\bm{g}_{ij}$ equal to the SSP template of age $t_{i}$ and metallicity
$z_{i}$, pre-convolved by the (known) LOSVD, and multiplied by the
discretisation interval widths, i.e.

$\displaystyle\bm{g}_{ij}=[g_{ij}(\lambda_{1}),...,g_{ij}(\lambda_{n})],$
$\displaystyle g_{ij}(\lambda)=\Delta t_{i}\Delta
z_{j}\int\frac{L(v)}{1+v/c}S\left(\frac{\lambda}{1+v/c},z_{j},t_{j}\right)\mathrm{d}v.$

Equation (2) now reduces to

$\bm{y}=\sum_{i=1}^{Z}\sum_{j=1}^{T}\bm{g}_{ij}f_{ij}+\bm{\epsilon},$ (7)

i.e. the data is simply a sum of (pre-convolved) templates $\bm{g}_{ij}$
weighted by the 2D discretised distribution function $f_{ij}$, plus noise
$\bm{\epsilon}$.

We correspondingly collect the templates $\bm{g}_{ij}$ into a 3D tensor
$\bm{G}$ of shape $(n,Z,T)$. Equation (7) can then be rewritten using tensor
notation, i.e.

$\bm{y}=\bm{G}\bm{f}+\bm{\epsilon},$

Corresponding to our choice to represent the image $\bm{f}$ as a matrix, the
prior covariance ${\bm{\Sigma}^{\bm{f}}}$ in (6) is given by a 4-tensor, where
the entry ${\bm{\Sigma}^{\bm{f}}}_{ij,kl}$ is equal to the covariance of the
pixel intensities $f_{ij}$ and $f_{kl}$. As we remarked above, different
choices for ${\bm{\Sigma}^{\bm{f}}}$ correspond to different regularisation
strategies. In this paper, we use the Ornstein-Uhlenbeck covariance

$\Sigma^{\bm{f}}_{ij,kl}=e^{-(\left|i-k\right|+\left|j-l\right|)/h}.$ (8)

This choice corresponds to the intuitive assumption that neighbouring pixels
in $\bm{f}$ are positively correlated, and that the amount of correlation
decreases exponentially with distance. The parameter $h>0$ is sometimes called
the “correlation length”, as it controls the rate of this decrease.

###  2.2 Mock data

Throughout Sections 3 \- 5, we will use mock data to illustrate our methods.
That is, we simulate a known ground truth LOSVD and distribution function
$\bm{f}$, apply the observation operator and add artificial noise to generate
mock, noisy data. Here we describe the specific choices taken for these steps.

For the LOSVD we will use a fourth order Gauss-Hermite expansion (van der
Marel & Franx, 1993) which takes parameters

$\bm{\theta}_{v}=(V,\sigma,h_{3},h_{4}),$

which correspond roughly to the mean, standard deviation, skewness and
kurtosis of $L(v)$. We will generate mock-data using plausible galactic
LOSVDs, namely using parameters $\bm{\theta}_{v}=[30,100,-0.05,0.1]$. We fix
$\bm{\theta}_{v}$ to the known ground truth during recoveries.

For the SSPs, we use the MILES222http://miles.iac.es/pages/ssp-models.php
models (Vazdekis et al., 2010) with BaSTI isochrones for base alpha-content
(Pietrinferni et al., 2004) and a Chabrier (2003) initial mass function. These
models span 53 ages in the range $0.03<t/\mathrm{Gyr}<14$ and 12 values of
metallicity $-2.27<z/\mathrm{[M/H]}<0.4$. Both age and metallicity values are
non-uniformly spaced in these ranges. The overall 2D grid of SSP models
consists of 636 template spectra. We use the wavelength range between
$[4700,6500]$ Angstroms with a sampling in velocity space of 10 km s-1.

For the distribution functions we use random Gaussian mixtures. We generate
ground-truth, light-weighted distributions $f$ as mixtures of three Gaussian
blobs which are randomly sized, shaped and positioned in the (non-uniformly
spaced) age and metallicity domain shown e.g. in the right panel of Figure 4.
Specifically, the parameters describing the blobs are sampled in normalised
co-ordinates which increase linearly from $[0,1]$ between the corners of the
domain. Blob means are sampled uniformly between $[0,1]$. The per-dimension
variances are sampled independently from log-uniform distributions between
$[0.01,1]$ and the blob covariance matrix is constructed using these variances
and a random rotation. The blob weights are sampled from log-uniform
distributions in the range $[1,30]$.

Finally, for the noise model we sample independent, normally distributed noise
with a constant standard deviation $\sigma^{2}$, i.e.

${\bm{\Sigma}^{\bm{y}}}=\mathbf{diag}\\{(\sigma^{2},...,\sigma^{2})\\},$ (9)

where $\sigma$ is chosen so that the average signal-to-noise ratio (SNR)
equals some specified value

$\frac{\left\|\bar{\bm{y}}\right\|}{\sigma}=\mathrm{SNR}.$

Unless otherwise specified, we will use $\mathrm{SNR}=100$.

##  3 The Laplacian-of-Gaussians method for blob detection

In this section, we recall the necessary prerequisites on the Laplacian-of-
Gaussians method. It is defined in terms of the Gaussian scale-space
representation, which we introduce first.

###  3.1 The scale-space representation of an image

The _scale-space representation_ of an image is a classic tool from computer
vision, where it is used to study image structures at their intrinsic scale.
For an exhaustive reference on scale-space theory in computer vision, see
Lindeberg (1994).

A scale-space representation of a continuous image
$f:\mathbb{R}^{2}\to\mathbb{R}$ is a one-parameter family $(L_{t})_{t\geq 0}$
of images with $L_{0}=f$, where the continuous parameter $t\geq 0$ corresponds
to physical scale. We will use the notation $L(\bm{x},t)=L_{t}(\bm{x})$ to
refer to the scale-space representation as a function both in ${\bm{x}}$ and
$t$.

One distinguishes linear and nonlinear scale space representations, where a
scale space is called linear if, for all $t\geq 0$, $L_{t}$ is obtained from
$f$ by a linear operation, and else is called nonlinear. The Gaussian scale
space is arguably the most well-studied scale space, since it is the unique
linear scale space that satisfies shift-invariance, scale invariance, the
semigroup property and non-enhancement of local extrema (see Florack et al.
(1992)). It is defined as

$L({\bm{x}},t)=k_{t}*f,$

where $k_{t}$ is a two-dimensional Gaussian of the form

$k_{t}=\frac{1}{2\pi t}e^{-\left\|x\right\|^{2}/(2t)}.$

###  3.2 The Laplacian-of-Gaussians method for blob detection

Given the Gaussian scale-space representation $L$ of a continuous image $f$,
the _Laplacian-of-Gaussians_ (LoG) method (see Lindeberg, 1998) detects blob-
like image structures as the local minimisers and maximisers of the scale-
normalised Laplacian of $L$, which is defined as

$\Delta^{\text{norm}}L({\bm{x}},t)=t(\partial_{x_{1}}^{2}+\partial_{x_{2}}^{2})L({\bm{x}},t).$
(10)

A local minimiser $({\bm{x}},t)\in\mathbb{R}^{2}\times[0,\infty)$ of
$\Delta^{\text{norm}}L$ corresponds to a bright blob of scale $t$ with center
${\bm{x}}$. On the other hand, local maximisers correspond to dark blobs. For
the rest of this article, we will only consider the identification of bright
blobs, since as explained in section 1, we expect them to correspond to
distinct stellar populations. While alternative ways to formally define blobs
in an image have been studied (e.g. in Lindeberg, 1994), those usually lead to
less straightforward methods for detection. For the rest of this paper, we
will always mean by a ”blob” of an image $f$ a local minimiser of
$\Delta^{\text{norm}}L$.

Note that the scale-normalisation in (10) is necessary to achieve scale-
invariance of the detected interest points: If a blob of scale $t$ is detected
in the image $f$, then it should be detected as a blob of scale $s\cdot t$ in
the rescaled image $\tilde{f}({\bm{x}})=f(\frac{{\bm{x}}}{s})$.

To build some intuition on the LoG method, consider the prototypical Gaussian
blob with scale parameter $s$ and center $\bm{m}$:

$f({\bm{x}})=\frac{1}{2\pi
s}e^{-\frac{\left\|{\bm{x}}-\bm{m}\right\|^{2}}{2s}}.$

It can then be shown that the scale-normalised Laplacian
$\Delta^{\text{norm}}L$ has exactly one local minimum at the point
$({\bm{x}},t)=(\bm{m},s)$. Hence, we can recover the scale and center of the
Gaussian blob by determining the local minimum of $\Delta^{\text{norm}}L$ in
the 3-dimensional scale space. Furthermore, this corresponds to the intuition
that blobs of an image correspond to points where the intensity of an image
has highly negative or positive curvature, corresponding to structures that
are bright or dark, respectively.

For a practical implementation, the LoG method has to be formulated for
discrete images $\bm{f}\in\mathbb{R}^{Z\times T}$ and given a discrete set of
scales $0\leq t_{1}<\ldots<t_{K}$. Analogously to the continuous setting, the
scale-space representation of a discrete image $\bm{f}$ is obtained from a
discrete convolution,

$\bm{L}_{t}=\bm{k}_{t}*\bm{f}.$

As convolution kernel $\bm{k}_{t}$ a sampled Gaussian is often used. However,
this choice is not consistent with the scale-space axioms (see e.g. section
1.2.2. in Weickert, 1998), as discussed by Lindeberg (1990). Instead, the
correct analogue of the continuous Gaussian kernel in the sense of scale-space
theory has to be defined in terms of the modified Bessel function, i.e.

$\displaystyle k_{t,ij}$ $\displaystyle=$ $\displaystyle
e^{-t^{2}}I_{i}(t)I_{j}(t),$ (11) $\displaystyle I_{i}(t)$ $\displaystyle=$
$\displaystyle\sum_{\nu=0}^{\infty}\frac{1}{\nu!\Gamma(\nu+i+1)}\left(\frac{t}{2}\right)^{2\nu+i}.$
(12)

With this choice, the scale-space properties from the continuous setting hold
exactly in the discrete setting (Lindeberg, 1991). We have observed in our
numerical tests that the use of this kernel leads to slightly better results
than with the sampled Gaussian.

Given a discrete set of scales $0\leq t_{1}<\ldots<t_{K}$, we obtain a
discrete scale-space representation of the image $\bm{f}$ as three-dimensional
matrix $\bm{L}=(\bm{L}_{t_{q}})_{q=1}^{K}\in\mathbb{R}^{Z\times T\times K}$.
Similar to the continuous setting, blobs of the discrete image $\bm{f}$ are
then detected as the local minimizers of
$\bm{\Delta}^{\text{norm}}\bm{L}=(\bm{\Delta}^{\text{norm}}\bm{L}_{t_{q}})_{q=1}^{K}\in\mathbb{R}^{Z\times
T\times K}$, where $\bm{\Delta}^{\text{norm}}$ is the discrete scale-
normalised Laplacian operator,

$\bm{\Delta}^{\text{norm}}\bm{L}_{ijt}=t\cdot(L_{i-1,j,t}+L_{i+1,j,t}+L_{i,j-1,t}+L_{i,j+1,t}-4L_{ijt}),$

which is well-defined for all $(i,j)\in\\{1,\ldots,Z\\}\times\\{1,\ldots,T\\}$
by reflecting the image at the boundary. The numerical value of
$\left|\bm{\Delta}^{\text{norm}}\bm{L}\right|$ at a detected blob can be used
as surrogate for the blob strength.

It is recommended to post-process the results of the LoG method by removing
blobs that are too weak or overlap too much. First, all blobs for which the
value of $\left|\bm{\Delta}^{\text{norm}}\bm{L}\right|$ is below a certain
threshold $l^{\mathrm{tresh}}$ are removed. In the present work, we use a
relative threshold of the form
$l^{\mathrm{tresh}}=r^{\mathrm{tresh}}\cdot\left|\min\Delta^{\text{norm}}\bm{L}^{\text{MAP}}\right|$,
where $r^{\mathrm{tresh}}\in(0,1)$. Without thresholding, the LoG method would
detect blobs at the slightest fluctuations of the image intensity, which would
often not correspond to meaningful features.

In a second step, the relative overlap of the detected blobs is computed.
Hereby, a blob $\bm{b}=(i,j,q)$ is identified with a circle of radius
$2\sqrt{t_{q}}$, centered at the pixel $(i,j)$. If the relative overlap of two
blobs surpasses a treshold $o^{\mathrm{tresh}}_{1}$, the blob with the smaller
value of $\left|\bm{\Delta}^{\text{norm}}\bm{L}\right|$ (i.e. the weaker blob)
is removed.

We visualise a detected blob $(i,j,t)$ by a circle centered at the pixel
$(i,j)$ and radius proportional to $\sqrt{t}$. A common choice for the radius
is $\sqrt{2t}$, motivated by the fact that such a circle contains 68% of the
total mass of a Gaussian with center $(i,j)$ and scale parameter $t$.

In figure 2, we plot the results of applying the LoG method to a MAP estimate
of the stellar distribution function from mock data. In this particular case
we used the relative threshold $r^{\mathrm{thresh}}=0.02$ for the blob
strength and the threshold $o^{\mathrm{thresh}}=0.5$ for the maximal relative
overlap of two blobs.

Before we move on, we want to point out that while the post-processing step is
important, it should not be interpreted as a statistical procedure since it
does not take the measurement error into account but only depends on the given
point estimate of the distribution function. For example, a blob in the MAP
estimate might be weak (i.e. a low value of
$\left|\bm{\Delta}^{\text{norm}}\bm{L}\right|$) but its presence could be very
certain. Vice versa, a very strong blob
($\left|\bm{\Delta}^{\text{norm}}\bm{L}\right|$ large) in the estimate might
at the same time be very uncertain.

##  4 An SVD-MCMC method for scalable uncertainty quantification

We wish to generalise the LoG method for blob detection to handle uncertainty.
As a prerequiste to this, we first must explain how to quantify uncertainties
for our problem. This poses some distinct challenges due to the combination of
ill-posedness, high-dimensionality and nonnegativity constraints, for which
generic methods for uncertainty quantification are not well-suited.

In this section, we present a specialised MCMC method that is able to overcome
the computational challenges in the reconstruction of stellar distribution
functions. To make MCMC sampling computationally feasible for this problem we
employ the following strategy: reduce the dimensions of the problem, sample
the posterior in the reduced space, then deproject the solution back to the
original space. The specific dimensionality reduction technique we will use is
based on truncating the singular value decomposition (SVD) which is closely
related to principal component analysis (PCA). This basic idea follows Tipping
& Bishop (1999) however our method for deprojection is substantially different
since we sample possible deprojections rather than take a single choice. In
the context of stellar population modelling, previous work includes Falcón-
Barroso & Martig (2021) who also combine SVD and MCMC sampling. Here, we build
on this work by outlining how to deproject samples from the low-dimensional
space in order to perform inference on stellar populations.

###  4.1 Notation

In this section it will be useful to rewrite our model in a vectorised form,
i.e. where we reshape 2D matrices to 1D vectors, 3D tensors to 2D matrices
etc. Rather than introduce new baroque notation, for the duration of Section 4
we ask the reader to reinterpret existing notation in this fashion.
Specifically, the foward model in (7), i.e.

$\bm{y}=\bm{G}\bm{f}+\bm{\epsilon},$ (13)

originally had $\bm{G}$ as a 3D tensor of shape $(n,Z,T)$ and $\bm{f}$ as the
2D image of shape $(Z,T)$. Here we abuse notation and instead use $\bm{f}$ to
refer to the 1D vectorised image of length $Z\cdot T$ and $\bm{G}$ to be
correspondingly reshaped into a 2D matrix of shape $(n,Z\cdot T)$. Similarly
the prior covariance $\Sigma_{\bm{f}}$ from equation (6) should now be
considered a 2D covariance matrix rather than a 4-tensor.

###  4.2 SVD regression

The SVD is a matrix decomposition which generalises eigenvalue decomposition
to rectangular matrices. For example, given a matrix $\bm{G}$ with shape
$(n,p=ZT)$, the SVD uniquely factorises $\bm{G}$ into three matrices,

$\bm{G}=\bm{U}\bm{S}\bm{V}^{T}.$

The matrix $\bm{S}$ is a diagonal matrix whose entries are known as singular
values, which are non-negative and are sorted in decreasing order, i.e.
$\bm{S}[0,0]\geq\bm{S}[1,1]\geq...\geq 0$. These are analogous to eigenvalues
of an eigendecomposition. The columns of $\bm{U}$ and $\bm{V}$ contain
orthonormal basis-vectors of the spaces $\mathbb{R}^{n}$ and $\mathbb{R}^{p}$
respectively, which are analogous to eigenvectors.

The SVD can be used for dimensionality reduction if the singular values decay
sufficiently fast. The _truncated SVD_ retains only the first $q$ components
of the factorisation, i.e.

$\bm{U}_{q}=\bm{U}[:,:q],\;\;\bm{S}_{q}=\bm{S}[:q,:q],\;\;\bm{V}_{q}=\bm{V}[:,:q].$
(14)

The total squared-error incurred by approximating $\bm{G}$ by the truncated
SVD is given by the sum of squares of the remaining $p-q$ singular values,
i.e.

$\left\|\bm{G}-\bm{U}_{q}\bm{S}_{q}\bm{V}_{q}^{\top}\right\|_{2}^{2}=\sum_{i=q+1}^{p}\bm{S}[i,i]^{2}.$

Therefore, if the singular values decay sufficiently fast then the truncated
SVD can provide an accurate approximation of $\bm{G}$ for relatively small
values of $q$. For the MILES SSP models (Vazdekis et al., 2010) we find that
$q=15$ provides a sufficiently good approximation. In general this choice will
depend on the particular SSP models used and the data quality, and in appendix
B we suggest a general strategy to choose $q$.

The truncated SVD transforms our problem to a low-dimensional regression
problem. We introduce the notation

$\bm{Z}=\bm{U}_{q}\bm{S}_{q}$ (15)

and insert into equation (13) to give the transformed equation for the forward
model,

$\displaystyle\bm{y}\approx\bm{Z}\bm{\eta}+\bm{\epsilon}\;\;\text{where}$ (16)
$\displaystyle\bm{\eta}=\bm{V}_{q}^{\top}\bm{f}.$ (17)

The original $p$-dimensional parameter vector $\bm{f}$ in equation (13) has
been reduced to a $q$-dimensional vector $\bm{\eta}$ in equation (16). To do
inference on $\bm{f}$, we first do inference on $\bm{\eta}$ using equation
(16), then deproject back to infer $\bm{f}$. Having reduced the dimensions and
the condition-number of the problem, this problem is now computationally more
amenable to MCMC sampling to estimate the posterior distribution
$p(\bm{\eta}|\bm{y})$.

###  4.3 Deprojection

How should we perform the deprojection? The suggestion from Tipping & Bishop
(1999) is to use the least-square deprojection of equation (17), i.e. for each
sample $\bm{\eta}$ from the MCMC, calculate

$\bm{f}_{*}=\bm{V}_{q}\bm{\eta}.$

This satisfies equation (17) because the SVD ensures that $\bm{V}_{q}$ is
orthonormal, i.e. $\bm{V}_{q}\bm{V}_{q}^{\top}=\bm{I}_{q}$. _The least-square
deprojection is not suitable for our problem, however, since it is not
guaranteed to satisfy the non-negativity constraint on $\bm{f}$_. In numerical
tests, these deprojections did contain many negative entries and are therefore
represent unphysical solutions. A solution to this problem would be to instead
perform a constrained deprojection, where we minimise some function
$\mathcal{R}(\bm{f})$ subject to the desired constraints,

$\bm{f}_{*}=\operatorname*{argmin}_{\bm{f}}\;\mathcal{R}(\bm{f})\;\text{subject
to}\;\bm{V}_{q}^{\top}\bm{f}=\bm{\eta}\;\text{and}\;\bm{f}\geq 0.$ (18)

A natural choice for the function $\mathcal{R}(\bm{f})$ is to use the same
regularising penalty function defined in equation (4). For standard, quadratic
choices of $\mathcal{R}(\bm{f})$, equation (18) is a quadratic programming
problem which can be solved efficiently. In practice, however, we find that
this solution is unacceptable since these deprojections are unphysically
spiky. Although the objective function in equation (18) promotes smoothness in
the deprojections, the equality constraint $\bm{V}_{q}^{\top}\bm{f}=\bm{\eta}$
_overwhelms_ this effect, resulting in unphysical deprojections too similar to
the spiky ML solution.

The final strategy we investigated is to again use MCMC sampling. Whereas
deprojection by calculating a point estimate will neglect uncertainty in the
$p-q$ SVD components which are not informed by the data, by sampling we can
account for this. These components will instead be informed by the prior
$p(\bm{f})$ defined in equation (5), i.e. the Bayesian equivalent of the
regularisation term. Since the strength of this prior depends on a tunable
regularisation parameter $\beta$, samples deprojected in this way can be made
as smooth as desired, in analogy with MAP vs ML recoveries.

###  4.4 SVD vs. PCA

SVD and PCA are closely related techniques. If we consider the matrix $\bm{G}$
to be a data matrix with $n$ samples each of dimension $p$, and then we
column-center this matrix, i.e.

$\displaystyle\hat{\bm{G}}=\bm{G}-[{\bm{\mu}},...,{\bm{\mu}}]\;\;\mathrm{where}$
(19) $\displaystyle{\bm{\mu}}=\frac{1}{p}\sum_{i=1}^{p}\bm{g}_{i}$ (20)

and $\bm{g}_{i}$ are the columns of $\bm{G}$, then the SVD of $\hat{\bm{G}}$,

$\hat{\bm{G}}=\hat{\bm{U}}\hat{\bm{S}}\hat{\bm{V}}^{\top},$

is equivalent to PCA of the data. They are equivalent in the sense that the
matrix $\hat{\bm{V}}$ contains the principal axes of the data and the squares
of the singular values in $\hat{\bm{S}}$ acquire the interpretation of being
the variances explained by each principal component. In practice we find that
using PCA works better than performing an SVD directly on $\bm{G}$, so we
adopt this choice throughout this work. We retain the name _SVD-MCMC_ since
the SVD interpretation helps to clarify the underlying ideas, and to prevent
confusion with the more common usage of PCA as applied to observed spectra
rather than model templates (e.g. Yip et al., 2004).

Centering $\bm{G}$ as in equation (19) transforms the forward model in (13) as
follows,

$\bm{y}=m{\bm{\mu}}+\hat{\bm{G}}\bm{f}+\bm{\epsilon},$

where $m$ is a scalar constant subject to the constraint

$m=\sum_{i=1}^{p}f_{i}.$

We define the truncated SVD of $\hat{\bm{G}}$ equivalently to equation (14) to
give $\hat{\bm{U}_{q}},\hat{\bm{S}}_{q}$ and $\hat{\bm{V}_{q}}$, then define
$\hat{\bm{Z}}$ equivalently to equation (15) to give the reduced forward model
after PCA,

$\displaystyle\bm{y}\approx
m{\bm{\mu}}+\hat{\bm{Z}}\bm{\eta}+\bm{\epsilon}\;\;\text{where}$
$\displaystyle\bm{\eta}=\hat{\bm{V}}^{\top}\bm{f}.$

###  4.5 Probabilistic model and MCMC implementation

We can now define our probabilistic models and describe how we sample their
posterior distributions. A probabilistic model is defined as a series of
sampling statements and transformations of variables required to generate the
observed data. Combined with Bayes’ theorem, the model implicitly defines an
expression for the posterior probability function of the model unknowns given
the observed data. An MCMC algorithm evaluates this expression to draw samples
from the posterior.

We will compare two models, the _Full MCMC_ model with no dimensionality
reduction, and the _SVD-MCMC_ model using the truncated SVD. Both models start
by sampling the prior distribution on $\bm{f}$, i.e.

$p(\bm{f})\propto\cases{}{N}\left(\bm{0},\frac{1}{\beta}{\bm{\Sigma}^{\bm{f}}}\right),&\text{if
}\bm{f}\geq 0,\\\ 0,\text{otherwise}.{}$ (21)

The two models are:

  1. 1.

Full MCMC:

$\displaystyle\bm{f}$ $\displaystyle\sim p(\bm{f})$
$\displaystyle\bm{y}\leavevmode\nobreak\ |\leavevmode\nobreak\ \bm{f}$
$\displaystyle\sim\mathcal{N}(\bm{G}\bm{f},{\bm{\Sigma}^{\bm{y}}}).$

  2. 2.

SVD-MCMC (for some choice of $q$):

$\displaystyle\bm{f}$ $\displaystyle\sim p(\bm{f})$ $\displaystyle\bm{\eta}$
$\displaystyle=\hat{\bm{V}}_{q}^{\top}\bm{f},$ $\displaystyle m$
$\displaystyle=\sum_{i=1}^{p}f_{i},$ $\displaystyle\bm{y}\leavevmode\nobreak\
|\leavevmode\nobreak\ \bm{\eta},m$
$\displaystyle\sim\mathcal{N}(m{\bm{\mu}}+\hat{\bm{Z}}\bm{\eta},{\bm{\Sigma}^{\bm{y}}}).$

In both cases, the final sampling statement draws data $\bm{y}$ from a
multivariate normal with noise covariance ${\bm{\Sigma}^{\bm{y}}}$ as defined
in equation (9). Each model implicitly defines an expression for the posterior
density $p(\bm{f}|\bm{y})$, which we sample via MCMC.

For MCMC sampling, we use the Hamiltonian Monte Carlo algorithm (HMC, Neal,
2011) as implemented in the numpyro probabilistic programming language (Phan
et al., 2019). Unless stated otherwise, for all experiments shown we run a
single chain for 5000 warmup steps and 10,000 sampling steps. The warmup
parameters are left at the numpyro defaults. We check convergence with three
criteria: (i) all sampled parameters have an $\hat{R}$ statistic (Gelman &
Rubin, 1992) in the range $[0.95,1.05]$, (ii) all sampled parameters have
effective sample size (Geyer, 2011) of at least 100, and (iii) no divergent
transitions occur during sampling (Betancourt, 2016).

Note that for the truncated SVD model, we directly sample the posterior
$p(\bm{f}|\bm{y})$, i.e. the _deprojection_ occurs internally during sampling.
We did also try the alternative implementation of first sampling
$p(\bm{\eta},m|\bm{y})$ and then later deprojecting as a subsequent step. The
inferred posteriors and total computation times of the two implementations are
very similar, however the first step of the alternative is orders of magnitude
faster than the deprojection step. Modelling with SVD-truncated spectral
templates is therefore an appealing option for purely kinematic models, i.e.
where stellar populations can be considered as a nuisance effect and need not
be deprojected. This is the approach taken by Falcón-Barroso & Martig (2021).
If we do also wish to also perform inference over stellar populations,
however, we see no clear benefit in separating the two sampling steps. Indeed,
one possible downside to their separation is that sampling
$p(\bm{\eta},m|\bm{y})$ requires defining some arbitrary prior
$p(\bm{\eta},m)$. It is not possible for any such low-dimensional prior to be
consistent with (21), or any physical prior, since the non-negativity
constraint can only be satisfied in $p$-dimensions. Nevertheless, we note that
both implementations of SVD-MCMC enjoys a significant speed-up compared to the
full model.

##  5 Uncertainty-aware blob detection

How can we incorporate uncertainties into the LoG method for blob detection?
Our motivation for doing so is demonsrated in figure 2: here we see that the
LoG method identifies multiple blobs in the MAP estimate of the distribution
function. However, the MAP estimate depends on the measurement noise, and so
it is uncertain whether structures that are identified in this way also exist
in the true, unknown stellar distribution function. Since this question needs
to be addressed if one wants to make credible statements about the structure
of the stellar distribution function, one has to find a way to incorporate the
a-posteriori uncertainty of the stellar distribution function into the blob
detection method.

In this section, we show how this can be achieved by combining the ideas from
scale-space theory with the concept of multi-dimensional Bayesian credible
intervals, a central tool for understanding and visualising uncertainty in
imaging. This results in an uncertainty-aware variant of the Laplacian-of-
Gaussians (LoG) method, which we correspondingly refer to as _ULoG_.

In section 5.1, we will review Bayesian posterior credible intervals for
images and related notions. In section 5.2, we introduce a new type of
credible intervals that incorporates the notion of scale into the
quantification of posterior uncertainty. We will refer to these as _filtered
credible intervals_. Section 5.3 discusses how they can be estimated from MCMC
samples. We then describe how filtered credible intervals can be used to
define an uncertainty-aware variant of the Laplacian-of-Gaussians method,
ULoG. The details of this method are given in section 5.4.

###  5.1 Credible intervals for imaging

A central tool for Bayesian uncertainty quantification are _posterior credible
intervals_. In general, they do not coincide with frequentist _confidence
intervals_ and should be interpreted differently. Since credible intervals are
canonically only defined for scalar parameters, the definition has to be
adapted when the unknown quantity of interest is multi-dimensional, as in the
case of a discrete image.

Before we state the definition, it will be useful to introduce the following
notation for multi-dimensional intervals: Let
$\bm{\ell},\bm{u}\in\mathbb{R}^{Z\times T}$ be two images such that
$\ell_{ij}\leq u_{ij}$ for every pixel $(i,j)$. Then, we will refer to the
hyper-box

$[\bm{\ell},\bm{u}]:=\left\\{\,\bm{f}\in\mathbb{R}^{Z\times
T}:\,\bm{\ell}\leq\bm{f}\leq\bm{u}\,\right\\}\subset\mathbb{R}^{Z\times T}$
(22)

as a _multi-dimensional interval_ in $\mathbb{R}^{Z\times T}$. Given a
posterior density function $p(\bm{f}|{\bm{y}})$ and a small credibility
parameter $\alpha\in(0,1)$ (typical values are $\alpha=0.05$ and
$\alpha=0.01$, respectively representing 95%- and 99%-credibility), we then
call a set
$[\bm{f}^{\mathrm{low}},\bm{f}^{\mathrm{upp}}]\subset\mathbb{R}^{Z\times T}$ a
simultaneous $(1-\alpha)$-credible interval if

$P(\bm{f}\in[\bm{f}^{\mathrm{low}},\bm{f}^{\mathrm{upp}}]|{\bm{y}})\geq
1-\alpha.$ (23)

That is, the multi-dimensional interval
$[\bm{f}^{\mathrm{low}},\bm{f}^{\mathrm{upp}}]$ contains at least a given
percentage of the posterior probability mass. By looking at the plots of both
the lower bound $\bm{f}^{\mathrm{low}}\in\mathbb{R}^{Z\times T}$ and the upper
bound $\bm{f}^{\mathrm{upp}}\in\mathbb{R}^{Z\times T}$, we can obtain a visual
understanding of the posterior uncertainty. This is examplified in the second
row of figure 5. The first row shows a simulated ground truth from which mock
data was generated (see section 2.2). The second row shows the corresponding
MAP estimate $\bm{f}^{\text{MAP}}$, together with the lower bound
$\bm{f}^{\mathrm{low}}$ and upper bound $\bm{f}^{\mathrm{upp}}$ of
simultaneous 95%-credible intervals.

Note that the simultaneous credible intervals defined in equation (23) have to
be interpreted differently than _pixel-wise_ credible intervals which also
appear in the statistical literature. We discuss these differences and the
motivation behind our choice in appendix A.

In section 5.3, we describe how we compute simultaneous credible intervals
using MCMC samples.

###  5.2 Filtered credible intervals

Figure 5: Visualisation of credible intervals from mock data. The upper panel
shows the ground truth from which the mock data was generated. In the second
row, the MAP estimate $\bm{f}^{\text{MAP}}$ (center column) together with the
lower bound $\bm{f}^{\mathrm{low}}$ (left) and upper bound
$\bm{f}^{\mathrm{upp}}$ (right) of the (unfiltered) simultaneous credible
intervals are plotted. Below that are the lower and upper bound of filtered
credible intervals at 3 different scales, together with the corresponding
filtered MAP estimate
$\bm{L}_{t}^{\mathrm{MAP}}=\bm{k}_{t}*\bm{f}^{\text{MAP}}$.

As a naive approach to perform uncertainty-aware blob detection we might try
to identify the blobs directly from the simultaneous credible intervals.
However, if the image of interest is highly uncertain, as is the case in the
present problem of stellar recoveries, such an approach is rather limited.
This can be seen from the first two rows of figure 5. One can observe that the
MAP estimate is quite close to the ground truth. In contrast, the credible
intervals are so wide that they only allow for very rough statements on the
shape of the stellar distribution function. The lower bound is close to zero
for most pixels. Clearly, the credible intervals alone do not accurately
represent the quality of our estimate.

This problem can be resolved by recognising that we are less interested in the
intensity values of specific pixels, but in higher-scale _structures_. Hence,
we are actually interested in coarse-scale information. It is an important
insight that – in contrast to the case where the image of interest is treated
as given – the uncertainty with respect to higher scales cannot be directly
derived from the uncertainties on finer scales. The correlation structure in
the posterior distribution might be more complex and exhibit long-range
dependencies, which means that the overall uncertainty cannot be reduced to
the pixel-by-pixel case.

Consequently, if we want more informative uncertainty quantification, we have
to integrate the notion of _scale_. This is achieved in our work by
integrating the ideas of scale-space representations (see also section 3.1):
Instead of computing credible intervals for the vector of pixel intensities
$\bm{f}$, we can perform uncertainty quantification subject to a scale
parameter $t$ by computing credible intervals for the smoothed image
$\bm{L}_{t}=\bm{k}_{t}*\bm{f}$, where $\bm{k}_{t}$ is the discrete analogue of
the Gaussian kernel, defined in equations (11)-(12). That is, we compute
images
$\bm{L}^{\mathrm{low}}_{t},\bm{L}^{\mathrm{upp}}_{t}\in\mathbb{R}^{Z\times T}$
such that

$P(\bm{L}_{t}\in[\bm{L}^{\mathrm{low}}_{t},\bm{L}^{\mathrm{upp}}_{t}]|{\bm{y}})\geq
1-\alpha.$

We will refer to $[\bm{L}^{\mathrm{low}}_{t},\bm{L}^{\mathrm{upp}}_{t}]$ as a
_filtered credible interval at scale $t$_.

In the last three rows of figure 5, filtered credible intervals at three
different scales are plotted. We see that the intervals tighten with
increasing scale, but loose the ability to represent small structures. This
visualises a very intuitive trade-off between scale and uncertainty: The
uncertainty is lower with respect to larger (high-scale) structures, since
they have less degrees of freedom. Correspondingly, the uncertainty with
respect to small-scale structures is larger, as the given image has much more
degrees of freedom at the small scale. The advantage of filtered credible
intervals is that they allow for an explicit choice. If one is only interested
in high-scale information, they allow to remove the superfluous uncertainty
with respect to small-scale variations. If one moves to a finer scale, one has
to pay the price of potentially larger uncertainty.

###  5.3  Estimating credible intervals from samples

Since in our case there are no analytical expressions for either unfiltered or
filtered credible intervals, we have to estimate both using MCMC samples.
After generating $S$ samples
$\bm{f}^{1},\ldots,\bm{f}^{S}\in\mathbb{R}^{Z\times T}$ from the posterior
distribution $p(\bm{f}|{\bm{y}})$, we estimate simultaneous
$(1-\alpha)$-credible intervals
$[\bm{f}^{\mathrm{low}},\bm{f}^{\mathrm{upp}}]$ by computing a multi-
dimensional interval that contains the desired percentage of samples
$\bm{f}^{i}$. That is, we find $\bm{f}^{\mathrm{low}}\in\mathbb{R}^{Z\times
T}$ and $\bm{f}^{\mathrm{upp}}\in\mathbb{R}^{Z\times T}$ such that at least
$(1-\alpha)\cdot S$ samples satisfy

$\displaystyle\bm{f}^{s}\in[\bm{f}^{\mathrm{low}},\bm{f}^{\mathrm{upp}}].$

Similarly, we estimate filtered credible intervals, at scale $t$, by computing
a multi-dimensional interval that contains the desired percentage of
_filtered_ samples. That is, we find
$\bm{L}^{\mathrm{low}}_{t},\bm{L}^{\mathrm{upp}}_{t}\in\mathbb{R}^{Z\times T}$
such that at least $(1-\alpha)\cdot S$ samples satisfy

$\displaystyle\bm{k}_{t}*\bm{f}^{s}\in[\bm{L}^{\mathrm{low}}_{t},\bm{L}^{\mathrm{upp}}_{t}].$

In both cases, the estimation of credible intervals reduces to the following
computational task: Given a sample of images
$\bm{x}^{1},\ldots,\bm{x}^{S}\in\mathbb{R}^{Z\times T}$, find the smallest
multi-dimensional interval $[\bm{\ell},\bm{u}]\subset\mathbb{R}^{Z\times T}$
that contains at least $(1-\alpha)\cdot S$ of these samples (with
$\bm{x}^{s}=\bm{f}^{s}$, this procedure estimates simultaneous credible
intervals; with $\bm{x}^{s}=\bm{k}_{t}*\bm{f}^{s}$, this procedure estimates
filtered credible intervals at scale $t$). This task is known as the
$k$-enclosing box problem in computational geometry. We use the following
heuristic but fast approach: First, find
$\bm{\ell},\bm{u}\in\mathbb{R}^{Z\times T}$ such that for every pixel $(i,j)$,
there are at least $(1-\alpha)\cdot S$ images for which

$\displaystyle\ell_{ij}\leq x^{s}_{ij}\leq u_{ij}$

(this is done using the order statistics of $(x^{s}_{ij})_{s=1}^{S}$). Since
the resulting multi-dimensional interval $[\bm{\ell},\bm{u}]$ usually contains
too few points, it is then rescaled until it encloses at least
$(1-\alpha)\cdot S$ points.

Note that there are several existing algorithms for solving the $k$-enclosing
box problem exactly or approximately (see Segal & Kedem (1998) and Chan & Har-
Peled (2021)). However, the computational complexity of these algorithms
increases at least quadratically with the number of samples and exponentially
with the underlying dimension, which renders them impractical for our high-
dimensional problem. We also considered a classic method by Besag et al.
(1995). However, in our numerical tests this method often led to credible
intervals that contained 100% of the samples, rendering them useless. For our
purposes, the heuristic approach was sufficient.

###  5.4 Uncertainty-aware Laplacian-of-Gaussians (ULoG)

In contrast to simultaneous credible intervals on the level of individual
pixels, the filtered credible intervals in figure 5 show that there are some
structures that are ”significant”, in the sense that they are visible both in
the image of the lower and the upper bound. This notion of significance is
scale-dependent, since some structures might only be visible in the filtered
credible intervals at sufficiently large scales. Intuitively, we would like to
say that a blob is ”credible” (at a given scale $t$) if all images
$\bm{h}\in[\bm{L}^{\mathrm{low}}_{t},\bm{L}^{\mathrm{upp}}_{t}]$ exhibit a
matching blob. Of course, one has to find a meaningful definition of what
”matching” means, as the exact position and size of detectable blobs will
usually vary among the images in
$[\bm{L}^{\mathrm{low}}_{t},\bm{L}^{\mathrm{upp}}_{t}]$. Since checking all
these images is computationally infeasible, we instead define a suitable
representative
$\bar{\bm{B}}_{t}\in[\bm{L}^{\mathrm{low}}_{t},\bm{L}^{\mathrm{upp}}_{t}]$
that has ”the least amount of blobs” in some quantitative sense. There are
many different loss functions that penalise the amount of blobs in an image.
We experimented with various loss functions of the form
$J(\bm{h})=\left\|\bm{D}\bm{h}\right\|^{2}$ with $\bm{D}$ being a discrete
differential operator. Among all the differential operators we tested, the
discrete Laplacian $\bm{\Delta}$ led to the most favorable properties. We
hence consider the representative
$\bar{\bm{B}}_{t}\in[\bm{L}^{\mathrm{low}}_{t},\bm{L}^{\mathrm{upp}}_{t}]$
given by the solution to the optimisation problem

$\bar{\bm{B}}_{t}=\operatorname*{argmin}_{\bm{h}\in\mathbb{R}^{Z\times
T}}\left\\{\,\left\|\bm{\Delta}\bm{h}\right\|_{2}^{2}:\,\bm{L}^{\mathrm{low}}_{t}\leq\bm{h}\leq\bm{L}^{\mathrm{upp}}_{t}\,\right\\}.$
(24)

This definition can be further motivated by the fact that in the LoG method,
blobs are related to the local minima and maxima of the discrete Laplacian of
an image (see section 3.2). Hence, penalising
$\left\|\bm{\Delta}\bar{\bm{B}}_{t}\right\|$ dampens the Laplacian and yields
a feature-poor representative. On an intuitive level, solving (24) can be
thought of as spanning a two-dimensional ”blanket” between the lower bound
$\bm{L}^{\mathrm{low}}_{t}$ and the upper bound $\bm{L}^{\mathrm{upp}}_{t}$.
In accordance with this intuition, we will refer to the resulting image
$\bar{\bm{B}}_{t}$ as a _$t$ -blanket_. An image of a $t$-blanket generated
from mock data is shown in figure 6.

To illustrate the intuition behind definition (24), a one-dimensional,
continuous example is shown in figure 7. Given a MAP estimate
$f^{\mathrm{MAP}}:[0,1]\to\mathbb{R}$ of an unknown signal, assume that a
posterior filtered credible interval is spanned by two functions
$L^{\mathrm{low}}_{t}:[0,1]\to\mathbb{R}$ and
$L^{\mathrm{upp}}_{t}:[0,1]\to\mathbb{R}$. The $t$-blanket
$\bar{B}_{t}:[0,1]\to\mathbb{R}$ is then the minimiser of
$\left\|h^{\prime\prime}\right\|$ over all functions $h$ that satisfy the
pointwise constraint $L^{\mathrm{low}}_{t}(x)\leq h(x)\leq
L^{\mathrm{upp}}_{t}(x)$ and the boundary conditions
$h^{\prime}(0)=h^{\prime}(1)=0$. The blanket $\bar{B}_{t}$ should not be
confused with the filtered MAP estimate $L^{\mathrm{MAP}}_{t}$, which will in
general be less smooth and might exhibit more blobs than $\bar{B}_{t}$, as
seen in figure 7.

Figure 6: Image of a $t$-blanket $\bar{\bm{B}}_{t}$ ($t=7$). The mock data is
the same as in figures 2 and 3. Two main structures are visible which
correspond to the two significant blobs visualised in figure 3. Figure 7:
Illustration of the idea behind definition (24). Given a filtered credible
interval (shaded region) determined by a lower bound $L^{\mathrm{low}}_{t}$
and an upper bound $L^{\mathrm{upp}}_{t}$, we find a function
$\bar{\bm{B}}_{t}$ between these bounds that minimises the number of
detectable blobs. The dotted, red graph $L^{\mathrm{MAP}}_{t}$ is the filtered
MAP estimate. This is shown for some dummy variable $x$.

If we do not treat the scale $t$ as given, we can study the three-dimensional
scale-space structure of the ”stack” of blankets
$\bar{\bm{B}}=(\bar{\bm{B}}_{t})_{t\geq 0}\subset\mathbb{R}^{Z\times T}$. We
observed in numerical tests that any given blob will typically be visible in
$\bar{\bm{B}}_{t}$ for a range of scales $[t_{1},t_{2}]\subset[0,\infty)$. At
the lowest scale $t_{1}$, it will typically be very weak, then grow stronger
as the scale increases (and the uncertainty decreases, correspondingly) until
it reaches a peak strength at some scale $t^{*}$. Then, as the scale increases
further towards $t_{2}$, the strength of the blob decreases again as it is
”smoothed out”. This phenomenon is somewhat analogous to the deterministic
theory of scale-space blobs (see Lindeberg, 1991). Motivated by this
similarity, we consider the blanket stack $\bar{\bm{B}}$ as an ”uncertainty-
aware” scale-space representation. In a way, it represents our posterior
uncertainty with respect to the unknown image $\bm{f}$ in scale space. In
analogy to the LoG method from section 3.2, we thus propose the heuristic
principle that ”significant” scale-space blobs correspond to local minimisers
of $\Delta^{\text{norm}}\bar{\bm{B}}$ in scale-space.

Instead of considering these significant blobs alone, we can obtain a more
nuanced view of our information on the unknown stellar distribution by
matching the set of candidate blobs detected in the MAP estimate
$f^{\text{MAP}}$ (see e.g. figure 2) to the significant blobs detected from
the blanket stack $\bar{\bm{B}}$. For example, if a blob is detected in the
MAP but does not correspond to a significant blob in $\bar{\bm{B}}$, it does
not mean that it can be ignored, but only that its presence is too uncertain
with respect to the given credibility level $1-\alpha$.

In summary, we arrive at the following uncertainty-aware variant of the LoG
method, which we refer to as ULoG. Given a discrete sequence of scales $0\leq
t_{1}<\ldots<t_{K}$, do the following:

  1. 1.

Perform blob detection on the MAP estimate $\bm{f}^{\text{MAP}}$ using the LoG
method described in section 3.2. Store the detected blobs in a set
$\mathcal{B}\subset\\{1,\ldots,Z\\}\times\\{1,\ldots,T\\}\times\\{1,\ldots,K\\}$.

  2. 2.

For $k=1,\ldots,K$, do:

    1. (a)

Compute a filtered credible interval
$[\bm{L}^{\mathrm{low}}_{t_{q}},\bm{L}^{\mathrm{upp}}_{t_{q}}]$.

    2. (b)

Using $[\bm{L}^{\mathrm{low}}_{t_{q}},\bm{L}^{\mathrm{upp}}_{t_{q}}]$, compute
the $t_{q}$-blanket $\bar{\bm{B}}_{t_{q}}$ by solving (24).

  3. 3.

Stack the blankets in a 3-dimensional array
$\bar{\bm{B}}=(\bar{\bm{B}}_{t_{q}})_{q=1}^{K}$.

  4. 4.

Determine the local minimisers of
$\bm{\Delta^{\text{norm}}}\bar{\bm{B}}\in\mathbb{R}^{Z\times T\times K}$ and
store them in a set $\mathcal{C}$.

  5. 5.

Perform a matching procedure (see section 5.4.1) that matches elements of
$\mathcal{B}$ to elements of $\mathcal{C}$.

  6. 6.

Visualise the results (see section 5.4.2).

####  5.4.1 Matching procedure

We have not yet described the details of the matching procedure mentioned in
the description of ULoG. There are of course many possibilities to define a
criterion for blobs ”matching”. For our purposes, the following procedure
proved successful:

  1. 1.

First, the detected blobs $\mathcal{B}$ in the MAP estimate are post-
processed, as described in section 3.2. Which blobs remain in $\mathcal{B}$
depends on the strength threshold $l^{\mathrm{thresh}}$ and the overlap
threshold $o_{1}^{\mathrm{thresh}}$.

  2. 2.

The scale-space minimisers $\mathcal{C}$ of
$\bm{\Delta^{\text{norm}}}\bar{\bm{B}}$ are cleaned in the same way, the only
difference being that $\bm{\Delta^{\text{norm}}}\bm{L}^{\text{MAP}}$ is
replaced with $\bm{\Delta^{\text{norm}}}\bar{\bm{B}}$.

  3. 3.

Now, for every $\bm{b}\in\mathcal{B}$ we check whether there is an overlap
with any $\bm{c}\in\mathcal{C}$. Here, overlap is defined in the same way as
in step $(ii)$. If there exists a $\bm{c}\in\mathcal{C}$ such that its
relative overlap with $\bm{b}$ is above a treshold $o^{\mathrm{tresh}}_{2}$,
then blobs $\bm{b}$ and $\bm{c}$ are a ”match”.

The ideal configuration for the tunable parameters in this procedure will
typically depend on the application. In our numerical experiments, we obtained
satisfying results using
$l^{\mathrm{thresh}}=r^{\mathrm{thresh}}\cdot\left|\min\Delta^{\text{norm}}\bm{L}^{\text{MAP}}\right|$
with $r^{\mathrm{thresh}}=0.02$ and
$o^{\mathrm{tresh}}_{1}=o^{\mathrm{tresh}}_{2}=0.5$.

####  5.4.2 Visualisation

There are multiple ways in which we could visualise the results of the
outlined algorithm. In figure 3, the following procedure was applied:

  * •

A scale-space blob $\bm{b}=(i,j,q)$ is visualised as a circle with center
$(i,j)$ and radius $\sqrt{2t_{q}}$ (see section 3.2).

  * •

If a blob $\bm{b}\in\mathcal{B}$ does not match a blob in $\mathcal{C}$, it is
visualised as a red circle.

  * •

If a blob $\bm{b}\in\mathcal{B}$ matches a blob in $\mathcal{C}$, it is
visualised as a green circle. Furthermore, the match $\bm{c}$ is visualised as
a dashed, blue circle.

The dashed blue circles have to be interpreted in the following way: If there
is a second structure inside the circle, then we cannot distinguish it with
confidence from the first blob. Hence, the radius of the circle indicates how
well the blob can be localised for the given credibility level $1-\alpha$. For
example, in figure 3 we see that two blobs on the right correspond to the same
blue circle. This means that we cannot separate them with credibility above
95%. We are confident that there is at least one blob in the true image
$\bm{f}$ that corresponds to the blue circle, but we cannot say with
credibility that there are two. Hence, these circles can be used as a tool to
visualise localisation uncertainty in the blob detection.

##  6 Application to M54

M54 is the nuclear star cluster of the Sagittarius dwarf galaxy and hosts a
mixture of stellar populations. Alfaro-Cuello et al. (2019) studied these
populations using resolved data of individual stars in the cluster. They
combined photometry from the Hubble Space Telescope with spectroscopic
measurements from the Multi-Unit Spectroscopic Explorer (MUSE) to measure ages
and metallicities of 6600 M54 stars, whose distribution is shown in the top-
left panel of Figure 8. Two distinct clumps are clearly identifiable: an old
metal-poor clump and young metal-rich clump, both of which display a trend of
increasing metallicity for younger stars. The authors hypothesised that the
old metal-poor population is formed from the accretion and merger of old and
metal poor globular clusters, while the young population is some combination
of stars formed in-situ in M54 and interlopers from the Sagittarius dwarf. We
test our blob detection techniques on M54, with these resolved measurements
providing a _ground truth_ distribution to compare with.

Boecker et al. (2020) (hereafter B20) used the M54 MUSE data to compare age-
metallicity recoveries from resolved measurements against those from
integrated spectra. They created an integrated M54 spectrum (SNR = 90) by co-
adding the spectra of individual M54 stars and analysed this using pPXF. The
resulting recovery is shown in the upper right panel of Figure 8. It matches
the ground-truth distribution in coarse statistics, with mean age and
metallicities from the two distributions separated by just 0.3 Gyr and 0.2 dex
$[\mathrm{Fe}/\mathrm{H}]$ respectively. The integrated recovery also displays
distinct populations, well separated in age and metallicity. We see a young
metal-rich population (single bright pixel at the top), and old metal-poor
contributions also.

One further difference between the B20 result and the ground truth
distribution is that the integrated recovery splits the old population into
two clumps. The authors quantified the uncertainty in their pPXF recovery by
Monte-Carlo resampling the data 100 times, randomly perturbed by the residual
from the baseline pPXF fit. They calculated 100 MAP estimates from the
resampled data, interpreting these as posterior samples. The pixel-wise 84%
credible intervals of these samples enveloped the two clumps in the old
population, suggesting that the split may not be significant. Our uncertainty-
aware blob detection techniques are well placed to address this question in a
principled way: our MCMC samples are directly drawn from the posterior, and
ULoG is a well-motivated way to do blob detection using these samples. We
therefore apply these techniques to the integrated M54 spectrum from B20 to
assess the significance of the recovered blobs.

Figure 8: Age-metallicity distributions for the M54 star cluster. Top left:
the ground truth number density from individual measurements of resolved stars
in M54, from Alfaro-Cuello et al. (2019). Top right: a MAP revovery using PPXF
full spectrum fitting of an integrated-light M54 spectrum, from B20. Bottom
row: Results of ULoG on M54 data using low-regularisation. Both panels show
the same MAP estimate of the recovered posterior probability density. The
coloured circles show the location of significant blobs with 95%-credibility,
with colours as in Figure 3. The left shows results using SVD-MCMC, the right
shows Full MCMC. In both cases, ULoG identifies 3 distinct blobs with
95%-credibility.

###  6.1 Notes on the previous works

In the Alfaro-Cuello et al. (2019) resolved study, the authors dissect the
metal-rich clump into two distinct populations. Whilst this extra split is
supported by spatial and kinematic information (Alfaro-Cuello et al., 2020),
based purely on the age-metallicity distribution seen in the upper left panel
of Figure 8 this split is not visible.

We noticed a small systematic effect in the B20 integrated-light study. In
their recovery, a non-zero weight was assigned to the youngest SSP (0.03 Gyr)
at an intermediate metallicity, with no counterpart in the resolved
distribution. This very young contibution was not visible in the B20 paper
since they use a linear axis for age, squashing very young SSPs to an
inperceptably narrow range. Our visualisation - with a image pixel per SSP -
highlighted this spurious contibution much more significantly. We believe this
is a systematic effect. The spectra of the very youngest SSPs, hosting the
hottest stars, are relatively featureless i.e. showing relatively few and weak
absorption features. A small contribution of very young stars may therefore be
used to offset residuals from incorrect flux calibration and/or treatment for
dust absorption. In order to focus our experiments on blobs with counterparts
in the ground truth, we therefore exclude SSPs below 0.1 Gyr in our
recoveries. After this cut, no weight is allocated to the youngest available
SSP (i.e. 0.1 Gyr), supporting the hypothesis that this very young
contribution was a systematic effect.

###  6.2 Setup

The inputs for our recoveries mirror those used in B20. We use the same
library of SSP templates i.e. E-MILES library Vazdekis et al. (2016), with an
bimodal IMF of slope 1.3, BaSTI isochrones (Pietrinferni et al., 2004) using
all available metallicities, but only retaining ages above 0.1 Gyr (see
Section 6.1). We mask the wavelengths ranges $[5850,5950],[6858.7-6964.9]$ and
$[7562.3-7695.6]$ (Angstrom) to avoid large sky-residuals. Since our analysis
routines are specialised to the case of the purely linear part of spectral
modelling, we first perform a fit with pPXF and fix the non-linear parts to
the pPXF best-fit values. We find best fitting kinematic parameters of
$[V=146.\mathrm{km/s},\sigma=3.\mathrm{km/s},h_{3}=0.088,h_{4}=0.187]$ and the
multiplicative polynomial of order 21 used by pPXF for continuum correction.
We pre-process the SSP templates to account for these effects, convolving them
with the LOSVD and multiplying them by the polynomial. Our Bayesian model was
set up using a truncated Gaussian prior (6) with an Ornstein-Uhlenbeck
covariance (8). The M54 integrated spectrum and noise estimates were provided
by the B20 authors. We try two choices for the regularisation parameter
$\beta$: a _low regularisation_ choice which gives a MAP reconstruction
similar to the B20 result ($\beta=1\cdot 10^{3}$), and a _high regularisation_
option for comparison ($\beta=5\cdot 10^{5}$).

We make the following specific choices for ULoG. To compute filtered credible
intervals, we Full MCMC and SVD-MCMC with $q=15$ (chosen as described in
section B). In both cases we follow the procedure laid out in section 4.5 i.e.
taking 10000 samples for $\bm{f}$ and checking that convergence diagnostics
are satisfied. To compute filtered credible intervals we consider 10 scales
$t_{k}=1+(k-1)\cdot 1.5$. The computation times (using Intel MKL on 8 3.6-GHz
Intel Core i7-7700 CPUs and 64 GB RAM) were as follows: for low/high
regularisation, SVD-MCMC required 1785/998 seconds, while Full MCMC required
12908/5295 seconds. This implies a speedup factor of roughly 5 to 10 for SVD-
MCMC over the full model, with greater speedup for less regularisation.

###  6.3 Results

Fits to the data for the low-regularisation case are shown in Figure 9. The
residuals are smaller than 10%, showing that the data is well-fit. The fact
that residuals are correlated between neighbouring wavelengths suggests
inaccuracies at the 10% level in the assumed, uncorrelated noise model. The
similarity between the SVD and Full-MCMC data-fits gives us confidence that
the SVD truncation has not biased the results. The data fits for the high-
regularisation case look very similar.

The recoveries for the low-regularisation case are shown in the bottom panels
of Figure 8. Results for both SVD-MCMC and Full MCMC again agree well. In both
cases, ULoG detects three distinct blobs: one young and metal rich, and two
old and metal poor. Importantly, we infer that these two old and metal poor
blobs are not likely to be part of a single, combined blob. Justification for
this can be seen clearly in Figure 10 which shows marginalised constraints on
the star formation histories: for old ages ¿5 Gyr, the 95% credible band are
inconsistent with a single blob. For this level of regularisation, our
computations clearly separate the old contribution into two metal-poor blobs
whereas B20 does not. This conclusion strongly depends on the regularisation
used, however.

Figure 9: Plots of the integrated M54 spectrum ${\bm{y}}$ and the data fit
${\bm{y}}^{\text{pred}}$ (in red), both for SVD-MCMC (left) and Full MCMC
(right). The data fit is related to the respective MCMC posterior mean
$\bm{f}^{\text{mean}}$ through
$\bm{y}^{\text{pred}}=\bm{G}\bm{f}^{\text{mean}}$. The plots in the second row
show the relative residuals
$r^{\text{rel}}_{i}=(y_{i}-y^{\text{pred}}_{i})/y_{i}$, where $\bm{y}$ is the
measurement data. The corresponding band (in blue) shows the simultaneous
posterior predictive 95%-credible intervals. Grey vertical bands show
wavelength regions which have been masked to avoid sky-residuals. Figure 10:
Marginal posterior for the age distribution (see also appendix A.1), for the
low-regularisation case. We show 95%-credible intervals (shaded area) and
marginal posterior mean (solid red line). Left: Based on SVD-MCMC. Right:
Based on Full MCMC.

Figures 11 and 12 show the inferred distributions using the high-
regularisation option. Though the young and metal rich blob is still detected
as significantly distinct, the two old and metal poor blobs have now merged
together. Additionally, an insignificant blob (red circle) is now seen at an
old and metal-rich population. The origin of this third blob is unclear,
however it is not present when we extend the SSP grid to include younger ages.
This is likely to be some systematic effect, however it is small and
ultimately deemed statistically insignificant by ULoG. This detection of only
two significant blobs more accurately matches the ground truth (top left panel
of Figure 8) which displays two rather than three distinct blobs. As with the
low-regularisation case, in both Figures 11 and 12 we see that the results of
both SVD-MCMC and Full MCMC agree well. Thus, while our experiments have
hopefully demonstrated the utility of ULoG (sped up with an SVD-MCMC backend)
the dependence of blob detection on regularisation parameter remains an open
question, which is discussed further in the next section.

Figure 11: Results of ULoG on M54 data with higher regularisation. Figure 12:
As for Figure 10 but with higher regularisation.

##  7 Discussion

The ideas explained in section 5 relate concepts from different fields such as
mathematics, statistics and image processing. Since ULoG uses scale-space
representations to perform statistical inference, it can be viewed as an
example of a statistical scale-space method (see e.g. Holmström & Pasanen,
2017). To our knowledge, the uncertainty-aware scale-space representation
based on ”blankets” (see section 5.4) is new to the literature. Nevertheless,
similar ideas have been used in variational imaging for a long time, although
not in the context of uncertainty quantification, and are known as ”taut
string” or ”tube methods” (see Mammen & van de Geer, 1997; Davies & Kovac,
2001; Hinterberger et al., 2003; Scherzer et al., 2009). Clarifying all of
these connections is a promising starting point for theoretical work, since
the resulting insights could be used to address various inference problems in
astronomical imaging.

Our implementation of ULoG is based on SVD-MCMC (see section 4). In this case,
the overall computational cost of the method is dominated by the generation of
posterior samples. Given these samples, the additional cost of computing
filtered credible intervals for a range of scales is negligible. Having access
to posterior samples, one might also consider alternative ways to perform
uncertainty-aware blob detection. One could try to circumvent the computation
of filtered credible intervals altogether and instead apply the LoG method (or
any other known blob detection algorithm) directly to each sample. However,
such a direct approach has to tackle a number of challenges: First, there is
no obvious way how to combine the results of the blob detection from each
individual sample into meaningful statistics. Different samples might have
blobs of widely differing size at non-matching locations. Furthermore,
individual posterior samples tend to be quite noisy. This is problematic,
since a blob detection method applied to a noisy sample will detect a large
number of ”spurious” blobs, making it very challenging to extract meaningful
information from its result.

However, ULoG has the key advantage that it only requires computation of
filtered credible intervals. This means that it does not explicitly need
access to posterior samples and hence can be combined with any method of
estimating filtered credible intervals. For example, in the context of linear
imaging problems Pereyra (2017) proposed to estimate posterior credible
regions using concentration-inequalities, circumventing the need for sampling.
Approaches based on this idea have been shown to be very scalable (see e.g.
Cai et al., 2018a). Using concentration-inequalities to directly estimate
filtered credible intervals could make ULoG applicable to very large imaging
problems. We will investigate this idea in future work.

###  7.1 SVD-MCMC

We have seen in section 6 that the choice of the regularisation parameter
$\beta$ in equation (5) affects the inferences. There are various strategies
to address the choice of $\beta$ in the Bayesian framework, which would be
straightforward to combine with SVD-MCMC. In a hierarchical approach (see e.g.
Molina et al., 1999; Pereyra et al., 2015), we could assign a prior on
$\beta$. Alternatively, we could estimate $\beta$ from the data using
empirical Bayesian methods (e.g. as in Pereira et al., 2015; Vidal et al.,
2020). We could also try to eliminate the presence of the regularisation
parameter altogether, and use a simulation-based prior (for an example in the
context of weak gravitational lensing, see Remy et al., 2022).

While in this paper we implemented SVD-MCMC with HMC, one could consider
different MCMC methods. For example, SVD-MCMC could be combined with
specialised samplers for nonsmooth probability densities, such as proximal
MCMC (Durmus et al., 2018), to allow e.g. the use of $\ell^{1}$-Laplace
priors, which is a common choice in imaging (see e.g. Caselles et al., 2015).

In our numerical tests, we observed that the HMC sampler ran into numerical
problems when the regularising prior was assigned to $\bm{f}$ in mass-weights
instead of light-weights. Specifically, a significant number of divergent
transitions occurred (Betancourt, 2016), and $\hat{R}$ statistics were much
greater than 1, signalling unconverged chains. This is likely due to the fact
that the observation operator $\mathcal{G}$ is much more ill-conditioned with
respect to mass-weighted distribution functions. This problem could be
resolved by replacing HMC with a more robust sampler, such as e.g. Randomize-
Then-Optimize (Bardsley et al., 2014), though it is unclear how efficient this
specific option would be for such a high-dimensional problem.

###  7.2 Applications beyond the current problem

In this work we have presented ULoG and SVD-MCMC on an idealised problem.
Several other effects may be important for spectral modelling, including
convolution with the LOSVD, treatment for reddening, and accounting for gas
emission. For the M54 experiment presented here, we have pre-processed the
data to remove these effects before proceeding with the linear stellar
population model. We note that in the future, these other effects could be
incorporated into a more complete model which is still fully compatible with
SVD-MCMC and ULoG.

Our method for uncertainty-aware blob detection is very general and may find
applications beyond spectral modelling. One specific example which we intend
to explore in the future is the task of blob detection in orbit distributions
inferred from orbit-superposition modelling of galaxies (e.g. Zhu et al.,
2022).

In our application to M54, we have seen that SVD-MCMC can lead to a reduction
of computation time up to a factor 7 compared to Full MCMC. For larger
problems, this factor might even increase as the full potential of the
truncated SVD comes truly into play in higher dimensions. Potentially, this
could make uncertainty quantification feasible not only for single-spectrum
fitting, but even for full IFU datacube modelling (see e.g. Hinterer et al.,
2023).

##  8 Conclusions

In this work we have adressed the question of blob detection in images with
uncertainties. We have proposed a novel, uncertainty-aware variant of the
classic Laplacian-of-Gaussians method for blob detection (ULoG). ULoG detects
significant blobs considering various scales, and requires as input
uncertainty estimates for the image in the form of filtered credible
intervals. Furthermore, we describe an efficient computational method for
calculating filtered credible intervals which combines truncated SVD and MCMC
sampling (SVD-MCMC).

These methods have been developed in the context of stellar population
modelling of integrated-light spectra. We illustrated the methods first on
mock data, and then applied it to real data of the nuclear star cluster M54.
Results using SVD-MCMC and Full MCMC match well. The observed speed-up using
SVD-MCMC is on average a factor 7, with sampling lasting 15-30 minutes as
opposed to 2-4 hours. The number of significant blobs detected ranged between
2-3, depending on the amount of regularisation used. This work focussed on
uncertainty-aware blob detection at fixed regularisation, but in the future
our methods can be combined with Bayesian techniques for selecting or
marginalising over the regularisation parameter.

###### Acknowledgements.

The authors thank the referee for a constructive and helpful report. PJ and
GvdV acknowledge funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme under grant
agreement No 724857 (Consolidator Grant ArcheoDyn). PJ and GvdV also
acknowledge support by the Austrian Science Fund (FWF): F6811-N36. FP and OS
were funded by the Austrian Science Fund (FWF): F6807-N36. The financial
support by the Austrian Federal Ministry for Digital and Economic Affairs, the
National Foundation for Research, Technology and Development and the Christian
Doppler Research Association is gratefully acknowledged. AB was funded by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID
138713538 – SFB 881 (“The Milky Way System”, subproject B08) and acknowledges
support from grant ProID2021010080 in the framework of Proyectos de I+D por
organismos de investigación y empresas en las áreas prioritarias de la
estrategia de especialización inteligente de Canarias (RIS-3), FEDER Canarias
2014-2020.

## References

  * Alfaro-Cuello et al. (2020) Alfaro-Cuello, M., Kacharov, N., Neumayer, N., et al. 2020, ApJ, 892, 20
  * Alfaro-Cuello et al. (2019) Alfaro-Cuello, M., Kacharov, N., Neumayer, N., et al. 2019, ApJ, 886, 57
  * Anderssen (1986) Anderssen, R. S. 1986, in Inverse Problems: Proceedings of the Conference Held at the Mathematical Research Institute at Oberwolfach, ed. J. R. Cannon & U. Hornung (Birkhäuser Basel), 11–30
  * Bardsley (2012) Bardsley, J. M. 2012, SIAM J. Sci. Comput., 34, A1316
  * Bardsley et al. (2014) Bardsley, J. M., Solonen, A., Haario, H., & Laine, M. 2014, SIAM J. Sci. Comput., 36, A1895
  * Besag et al. (1995) Besag, J., Green, P., Higdon, D., & Mengersen, K. 1995, Stat. Sci., 10
  * Betancourt (2016) Betancourt, M. 2016, arXiv e-prints, arXiv:1604.00695
  * Boecker et al. (2020) Boecker, A., Alfaro-Cuello, M., Neumayer, N., Martín-Navarro, I., & Leaman, R. 2020, ApJ, 896, 13
  * Breth (1978) Breth, M. 1978, Ann. Stat., 6
  * Cai et al. (2018a) Cai, X., Pereyra, M., & McEwen, J. D. 2018a, MNRAS, 480, 4170
  * Cai et al. (2018b) Cai, X., Pereyra, M., & McEwen, J. D. 2018b, MNRAS, 480, 4154
  * Calvetti & Somersalo (2018) Calvetti, D. & Somersalo, E. 2018, Wiley Interdiscip. Rev. Comput. Stat., 10, e1427
  * Cappellari (2017) Cappellari, M. 2017, MNRAS, 466, 798
  * Cappellari & Emsellem (2004) Cappellari, M. & Emsellem, E. 2004, PASP, 116, 138
  * Carnall et al. (2018) Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, MNRAS, 480, 4379
  * Caselles et al. (2015) Caselles, V., Chambolle, A., & Novaga, M. 2015, in Handbook of mathematical methods in imaging, 3
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Chan & Har-Peled (2021) Chan, T. M. & Har-Peled, S. 2021, Discrete & Computational Geometry, 66, 769
  * Cid Fernandes et al. (2005) Cid Fernandes, R., Mateus, A., Sodré, L., Stasińska, G., & Gomes, J. M. 2005, MNRAS, 358, 363
  * Dashti & Stuart (2017) Dashti, M. & Stuart, A. M. 2017, in Handbook of Uncertainty Quantification, ed. R. Ghanem, D. Higdon, & H. Owhadi (Springer International Publishing), 311–428
  * Davies & Kovac (2001) Davies, P. L. & Kovac, A. 2001, Ann. Stat., 29, 1
  * Despali et al. (2022) Despali, G., Vegetti, S., White, S. D. M., et al. 2022, MNRAS, 510, 2480
  * Durmus et al. (2018) Durmus, A., Moulines, É., & Pereyra, M. 2018, SIAM J. Imaging Sci., 11, 473
  * Falcón-Barroso & Martig (2021) Falcón-Barroso, J. & Martig, M. 2021, A&A, 646, A31
  * Florack et al. (1992) Florack, L. M. J., ter Haar Romeny, B. M., Koenderink, J. J., & Viergever, M. A. 1992, Image Vis. Comp., 10, 376
  * Ganesh (2009) Ganesh, N. 2009, J. Multivar. Anal., 100, 1610
  * Gelman & Rubin (1992) Gelman, A. & Rubin, D. B. 1992, Stat. Sci., 7, 457
  * Geyer (2011) Geyer, C. J. 2011, in Handbook of Markov Chain Monte Carlo
  * Gladders et al. (2013) Gladders, M. D., Oemler, A., Dressler, A., et al. 2013, ApJ, 770, 64
  * Groetsch (1983) Groetsch, C. W. 1983, in Improperly Posed Problems and Their Numerical Treatment, ed. G. Hämmerlin & K.-H. Hoffmann, International Series of Numerical Mathematics (Birkhäuser), 97–104
  * Hansen & O’Leary (1993) Hansen, P. C. & O’Leary, D. P. 1993, SIAM J. Sci. Comput., 14, 1487
  * Hanson (1993) Hanson, K. M. 1993, in Medical Imaging 1993: Image Processing, ed. M. H. Loew, Vol. 1898 (SPIE), 716–731
  * Hinterberger et al. (2003) Hinterberger, W., Hintermüller, M., Kunisch, K., von Oehsen, M., & Scherzer, O. 2003, J. Math. Imaging Vis., 19, 219
  * Hinterer et al. (2023) Hinterer, F., Hubmer, S., Jethwa, P., et al. 2023, SIAM J. Imaging Sciences, 16, 192
  * Holmström & Pasanen (2017) Holmström, L. & Pasanen, L. 2017, Int. Stat. Rev., 85, 1
  * Kaipio & Somersalo (2007) Kaipio, J. & Somersalo, E. 2007, J. Comput. Appl. Math., 198, 493
  * Koleva et al. (2008) Koleva, M., Prugniel, P., Ocvirk, P., Le Borgne, D., & Soubiran, C. 2008, MNRAS, 385, 1998
  * Leja et al. (2019) Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, J. S. 2019, ApJ, 876, 3
  * Li et al. (2015) Li, Y., Wang, S., Tian, Q., & Ding, X. 2015, Neurocomputing, 149, 736
  * Lindeberg (1990) Lindeberg, T. 1990, IEEE Trans. Pattern Anal. Mach. Intell., 12, 234
  * Lindeberg (1991) Lindeberg, T. 1991, PhD thesis, Royal Institute of Technology, Stockholm and Stockholm University
  * Lindeberg (1994) Lindeberg, T. 1994, Scale-Space Theory in Computer Vision (Boston: Springer)
  * Lindeberg (1998) Lindeberg, T. 1998, Int. J. Comput. Vis., 30, 79
  * Lindeberg (2013) Lindeberg, T. 2013, J. Math. Imaging Vis., 46, 177
  * Lower et al. (2020) Lower, S., Narayanan, D., Leja, J., et al. 2020, ApJ, 904, 33
  * Mammen & van de Geer (1997) Mammen, E. & van de Geer, S. 1997, Ann. Stat., 25, 387
  * Martin et al. (2012) Martin, J., Wilcox, L. C., Burstedde, C., & Ghattas, O. 2012, SIAM J. Sci. Comput., 34, A1460
  * Masias et al. (2012) Masias, M., Freixenet, J., Lladó, X., & Peracaula, M. 2012, MNRAS, 422, 1674
  * Mikolajczyk et al. (2005) Mikolajczyk, K., Tuytelaars, T., Schmid, C., et al. 2005, Int. J. Comput. Vis., 65, 43
  * Molina et al. (1999) Molina, R., Katsaggelos, A., & Mateos, J. 1999, IEEE Transactions on Image Processing, 8, 231
  * Neal (2011) Neal, R. 2011, in Handbook of Markov Chain Monte Carlo, 113–162
  * Ocvirk et al. (2006) Ocvirk, P., Pichon, C., Lançon, A., & Thiébaut, E. 2006, MNRAS, 365, 46
  * Parno et al. (2021) Parno, M., Davis, A., & Seelinger, L. 2021, JOSS, 6, 3076
  * Pereira et al. (2015) Pereira, A., Antoni, J., & Leclère, Q. 2015, Applied Acoustics, 97, 11
  * Pereyra (2017) Pereyra, M. 2017, SIAM J. Imaging Sci., 10, 285
  * Pereyra et al. (2015) Pereyra, M., Bioucas-Dias, J. M., & Figueiredo, M. A. T. 2015, in 2015 23rd European Signal Processing Conference (EUSIPCO), 230–234, ISSN: 2076-1465
  * Phan et al. (2019) Phan, D., Pradhan, N., & Jankowiak, M. 2019, arXiv preprint arXiv:1912.11554
  * Pietrinferni et al. (2004) Pietrinferni, A., Cassisi, S., Salaris, M., & Castelli, F. 2004, ApJ, 612, 168
  * Price et al. (2020) Price, M. A., Cai, X., McEwen, J. D., et al. 2020, MNRAS, 492, 394
  * Price et al. (2021) Price, M. A., McEwen, J. D., Cai, X., et al. 2021, MNRAS, 506, 3678
  * Remy et al. (2022) Remy, B., Lanusse, F., Jeffrey, N., et al. 2022, Probabilistic Mass Mapping with Neural Score Estimation, Tech. Rep. arXiv:2201.05561, arXiv
  * Repetti et al. (2018) Repetti, A., Pereyra, M., & Wiaux, Y. 2018, in 2018 26th European Signal Processing Conference (EUSIPCO), 2668–2672
  * Scherzer et al. (2009) Scherzer, O., Grasmair, M., Grossauer, H., Haltmeier, M., & Lenzen, F. 2009, Variational Methods in Imaging, subsequent edition edn. (Springer)
  * Segal & Kedem (1998) Segal, M. & Kedem, K. 1998, Inf. Process. Lett., 65, 95
  * Simha et al. (2014) Simha, V., Weinberg, D. H., Conroy, C., et al. 2014, arXiv e-prints, arXiv:1404.0402
  * Tipping & Bishop (1999) Tipping, M. E. & Bishop, C. M. 1999, J. R. Stat. Soc., B: Stat. Methodol., 61, 611
  * Tojeiro et al. (2007) Tojeiro, R., Heavens, A. F., Jimenez, R., & Panter, B. 2007, MNRAS, 381, 1252
  * van der Marel & Franx (1993) van der Marel, R. P. & Franx, M. 1993, ApJ, 407, 525
  * Vazdekis et al. (2015) Vazdekis, A., Coelho, P., Cassisi, S., et al. 2015, MNRAS, 449, 1177
  * Vazdekis et al. (2016) Vazdekis, A., Koleva, M., Ricciardelli, E., Röck, B., & Falcón-Barroso, J. 2016, MNRAS, 463, 3409
  * Vazdekis et al. (2010) Vazdekis, A., Sánchez-Blázquez, P., Falcón-Barroso, J., et al. 2010, MNRAS, 404, 1639
  * Vidal et al. (2020) Vidal, A. F., De Bortoli, V., Pereyra, M., & Durmus, A. 2020, SIAM J. Imaging Sci., 13, 1945
  * Weickert (1998) Weickert, J. 1998, Anisotropic Diffusion in Image Processing (Stuttgart: B. G. Teubner)
  * White & Rees (1978) White, S. D. M. & Rees, M. J. 1978, MNRAS, 183, 341
  * Yip et al. (2004) Yip, C. W., Connolly, A. J., Vanden Berk, D. E., et al. 2004, AJ, 128, 2603
  * Zhu et al. (2022) Zhu, L., van de Ven, G., Leaman, R., et al. 2022, arXiv e-prints, arXiv:2203.15822

##  Appendix A Simultaneous vs. pixel-wise credible intervals

In section 5.1 we defined credible intervals as _simultaneous_ intervals.
These need to be distinguished from _pixel-wise_ credible intervals that are
also often used for uncertainty visualisation (for scalar parameters, these
notions coincide). Given a pixel $(i,j)$, we call an interval
$[\ell_{ij},u_{ij}]\subset\mathbb{R}$ a _pixel-wise posterior
$(1-\alpha)$-credible interval_ if it is a $(1-\alpha)$-credible interval with
respect to the _marginalised_ posterior distribution $p(f_{ij}|{\bm{y}})$.
That is, if there holds

$P(f_{ij}\in[\ell_{ij},u_{ij}]|{\bm{y}})\geq 1-\alpha.$

If we assemble all lower bounds $(\ell_{ij})$ and upper bounds $(u_{ij})$ into
images $\bm{\ell},\bm{u}\in\mathbb{R}^{T\times Z}$, these images define then a
multi-dimensional interval $[\bm{\ell},\bm{u}]\subset\mathbb{R}^{T\times Z}$
via equation (22). However, this interval will in general not satisfy equation
(23), as pointed out e. g. by Ganesh (2009). For our purposes, simultaneous
credible intervals are more suitable, since by definition they satisfy
equation (23), thus allowing the intuitive interpretation that the signal of
interest lies ”between” the lower bound $\bm{f}^{\mathrm{low}}$ and the upper
bound $\bm{f}^{\mathrm{upp}}$ with posterior probability greater or equal
$1-\alpha$.

The difference between simultaneous and pixel-wise credible intervals can also
be stated by referring to the concept of _posterior credible regions_. We call
a set $C^{\alpha}\subset\mathbb{R}^{T\times Z}$ a posterior
$(1-\alpha)$-credible region for $\bm{f}$ if

$p(\bm{f}\in C^{\alpha}|{\bm{y}})\geq 1-\alpha.$

By equation (23), a simultaneous $(1-\alpha)$-credible interval is always a
$(1-\alpha)$-credible region, while this might fail for a pixel-wise
$(1-\alpha)$-credible interval.

We also remark that the terminology in this area is not completely
standardised. Credible regions and credible intervals are sometimes called
Bayesian confidence regions and Bayesian confidence intervals, respectively,
while the multi-dimensional simultaneous credible intervals mentioned above
have been referred to as Bayesian confidence bands or credible bands in the
area of time series analysis (see Breth 1978).

###  A.1 Credible intervals for the age distribution

The credible intervals for the age distribution plotted in figures 10 and 12
were computed as follows: For each MCMC sample $\bm{f}^{s}$, we sum over the
metallicity,

$\displaystyle a_{j}^{s}=\sum_{i=1}^{Z}f_{ij}^{s},\qquad s=1,\ldots,S.$

Then, we estimated simultaneous $95\%$-credible intervals for
$\bm{a}\in\mathbb{R}^{T}$ using the method outlined in section 5.3 (that is,
setting $\bm{x}^{s}=\bm{a}^{s}$).

##  Appendix B Choice of $q$

How can we choose an appropriate size of the latent dimension $q$? One
possibility is to look at the rate of decay of singular values. These are
shown for our problem in the top panel of Figure 13. Singular values decay
more slowly beyond $q=15$, hence the accuracy of the truncated SVD
approximation grows more slowly beyond this value. Choosing a value of $q$
around this inflection point provides a good balance between accuracy and
possible speed-up from dimensionality reduction. Note however that this choice
is not _guaranteed_ to be sufficient in all circumstances e.g. for very high
signal-to-noise data it is conceivable that a larger $q$ may be required. We
therefore recommend to check convergence of results against the choice $q$.

Figure 13: The influence of the latent dimension $q$ on SVD-MCMC. The top
panel shows singular values of the matrix of MILES SSPs. The lower left panel
shows the error between filtered credible intervals computed with SVD-MCMC and
Full MCMC. The error is measured in the mean Jaccard distance (see equation
(25)). Beyond $q=15$ the error stabilises, signalling converged results. The
lower right panel shows the corresponding computation time for SVD-MCMC, where
the dashed horizontal line marks the computation time for Full MCMC. In all
panels, the vertical dotted line is at $q=15$, the value adopted for SVD-MCMC
inferences in this work.

To assess convergence we will run SVD-MCMC over a range of $q$ and quantify
changes in the results. The metric we will use to compare two sets of MCMC
samples is the _mean Jaccard distance_ between filtered credible intervals
computed from the samples. For two intervals
$[u^{1},v^{1}],[u^{2},v^{2}]\subset\mathbb{R}$, we define their _Jaccard
distance_ as

$\quad
D_{J}([u^{1},v^{1}],[u^{2},v^{2}])=1-\frac{\left|[u^{1},v^{1}]\cap[u^{2},v^{2}]\right|}{\left|[u^{1},v^{1}]\cup[u^{2},v^{2}]\right|}.$

The value of the Jaccard distance is always between $0$ and $1$, where $0$
corresponds to identity and $1$ corresponds to disjointness (up to sets of
measure zero). To measure the distance between filtered credible intervals
$[\bm{u}^{1},\bm{v}^{1}],[\bm{u}^{2},\bm{v}^{2}]\subset\mathbb{R}^{T\times Z}$
in image space, we take the mean Jaccard distance over all pixels:

$\overline{D}_{J}([\bm{u}^{1},\bm{v}^{1}],[\bm{u}^{2},\bm{v}^{2}]):=\frac{1}{T\cdot
Z}\sum_{i=1}^{T}\sum_{j=1}^{Z}D_{J}([u^{1}_{ij},v^{1}_{ij}],[u^{2}_{ij},v^{2}_{ij}]).$
(25)

The lower left panel of figure 13 shows the results of our convergence test on
mock data. It shows the mean Jaccard distance between filtered credible
intervals computed using Full MCMC and SVD-MCMC with 5000 warmup and 5000
sampling steps as a function of $q$. The filtered credible intervals were
computed for scale $t=5$ and credibility parameter $\alpha=0.05$, but we
observed that the choice of these parameters did not significantly affect the
resulting plot. Beyond $q=15$ the distance does not decrease, signalling that
filtered credible intervals do not systematically change and our results are
converged. The minimum distance does not reach zero, but instead flattens
around $0.05$ due to the inherent stochasticity of samples (i.e. two
independent sets of samples from the Full MCMC model would also show this
level of difference). The bottom right panel of figure 13 shows the
computation time required for sampling the SVD-MCMC models as a function of
$q$. At $q=15$, we still benefit from a factor 15 speedup compared to sampling
the full model.

For general usage of SVD-MCMC (where it may be infeasible/undesirable to
compare against Full MCMC), we suggest the following strategy for picking $q$.
Pick some initial $q_{0}$ smaller than the inflection point of the singular
values and a step-size $\Delta q\approx 5$. Perform SVD-MCMC inferences at
$q=q_{0}+i\Delta q$, first for $i=0$ and $1$, and compute the Jaccard distance
between these two sets of results. If $q_{0}$ is small enough, there will
likely be a significant decrease in Jaccard distance between $i=0$ and $1$.
Finally, increment $i$, computing the Jaccard distance between the cases $i$
and $i-1$, till the distance stabilises to a constant value.