Shared Stochastic Gaussian Process Latent Variable Models: A Multi-modal Generative model for Quasar spectra
Authors: Vidhi Lalchand
Anna-Christina Eilers
Venue: Accepted by TMLR ()

Abstract: This work proposes a scalable probabilistic latent variable model based on Gaussian processes in the context of multiple observation spaces. We focus on an application in astrophysics where it is typical for data sets to contain both observed spectral features as well as scientific properties of astrophysical objects such as galaxies or exoplanets. In our application, we study the spectra of very luminous galaxies known as quasars, and their properties, such as the mass of their central supermassive black hole, their accretion rate and their luminosity, and hence, there can be multiple observation spaces. A single data point is then characterised by different classes of observations which may have different likelihoods. Our proposed model extends the baseline stochastic variational Gaussian process latent variable model (GPLVM) to this setting, proposing a seamless generative model where the quasar spectra and the scientific labels can be generatedsimultaneously when modelled with a shared latent space acting as input to different sets of Gaussian process decoders, one for each observation space. Further, this framework allows training in the missing data setting where a large number of dimensions per data point may be unknown or unobserved. We demonstrate high-fidelity reconstructions of the spectra and the scientific labels during test-time inference and briefly discuss the scientific interpretations of the results along with the significance of such a generative model.

Published in Transactions on Machine Learning Research (Jan/2025)
Shared Stochastic Gaussian Process Latent Variable Models: A Multi-modal Generative Model for Quasar Spectra
Vidhi Lalchand
vidrl@mit edu
Eric and Wendy Schmidt Center Broad Institule of MIT and Harvard; MIT
Anna-Christina Eilers
cilers Cemit edu
MIT Kooli Institute for Astrophysics and Space Research
Reviewed on OpenReview: https: // shorturl. at/ at/SUgLP
Abstract
This work proposes scalable probabilistic latent variable model based on Gaussian processes Lawrence 2004) in the context of multiple observation spaces. We focus on an application in astrophysics where it is typical for data sets to contain both observed spectral features as well scientific properties of astrophysical objects such galaxies or exoplanets. In our application, study the spectra of very luminous galaxies known as quasars, and their properties, such as the mass of their central supermassive black hole, their accretion rate and their luminosity, and hence, there can be multiple observation spaces. single data point is then characterised by different classes of observations, which may have different likelihoods. Our proposed model extends the baseline stochastic variational Gaussian process latent variable model (GPLVM) (Lalchand et al., 2022) to this setting, proposing scamless generative model where the quasar spectra and scientific labels can be generated simallaneeusly when modelled with a shared latent space acting as input to different sels of Gaussian process decoders, one for each observation space. In addition, this framework allows training in the missing data setting whore a large number of dimensions per data point may be unknown or unobserved. We demonstrate high-fidelity reconstructions of the spectra and the scientific labels during test-time inference and briefly discuss the scientific interpretations of the results along with the significance of such a generative model. The code for this work is available at htttps://githii cow wwwp:///////////////// Baa.000//uaaar-GGPPGPG
1 Introduction
Many challenges in the contemporary physical sciences arisc from the analysis of largo-scale, noisy, het- eroscedastic, and high-dimensional datasets (Clarke et al. 2016). Hence, there is increasing consensus that addressing the challenges posed by large-scale data in the experimental pipelines for discovery, forecasting, and prediction warrant scalable machine learning. Modern experiments in physics, chemistry, and astronomy are capable of producing extremely complex high-dimensional data, but it not just the sheer volume the data, but also the velocily, referring to the rate of production (Carleo et al.] 2019) of data, which poses an additional challenge. Further, an additional axis is the variety or heterogencity inherent in scientific datasets in the form of multiple outputs or observation spaces. For instance, images (pixels) can be accompanied with attributes like illumination, pose, and resolution. In addition, some of these attributes observed features are often unknown or unobserved, resulting in incomplete data sets with missing components. The multiple output setting is inadequately addressed in the probabilistic machine learning literature; in this work we tackle precisely this setting proposing a multi-output. scalable probabilistic model for a scientific application in astronomy.
A preponderance of literature in unsupervised learning focuses on the setting of a single large-scale homogeneous dataset, for instance, images, text, speech, continuous numerical values. Further, most parametric models

---

Published in Transactions on Machine Learning Research (Jan/2025)
also assume a complete dataset where each dimension of every data point is observed. Real-world data are often only partially observed and in some cases yield very sparse data matrices with majority of the dimensions missing (Corduneanu & Jaakkola, 2012). Robust unsupervised learning in these settings is a challenge and one needs to account for sensible epistemic uncortaintics. In the running application we consider in this work, handling missing data during training is crucial, as spectral pixels real observations can be missing due to absorption features in the Earth's atmosphere or different wavelength coverage of the spectra when observed with different telescopes or spectrographs.
Dimensionality reduction is standard precursor to further analysis in scientific datasets, as the intrinsic dimensionality of the data might actually be quite low. It is usually possible to summarise the key axis of variation in the data with very few dimensions (Facco et al.) 2017) side-stepping the curse of dimensionality and facilitating further downstream analysis. Projection techniques like PCA, multidimensional scaling (MDS) and independent component analysis (ICA) utilise eigenvalue decomposition (Murphy 2012) while other non-lincar techniques like t-SNE (van der Maaten & Hinton 2008 SNE (van der Maaten 2009) and UMAP (McInnes et. al. 2020) construct a probability distribution on high dimensional points and replicate a similar distribution on low-dimensional points iteratively using the KL divergence. These methods, however, are not generative in their traditional incarnation. They cannot be used to generate instances of the high-dimensional data points; hence, they are not very useful in many astrophysical settings where the ability to generate points in high-dimensional data space is crucial.
Generative latent variable models supplement traditional dimensionality reduction techniques as they offer the simultancous benefits of probabilistic interpretation and data generation while learning faithful embedding of the high-dimensional training data in low-dimensional latent space. A generative probabilistic framework like the GPLVM (Lawrence) 2004) works by optimising the parameters of a Gaussian process decoder from a low dimensional latent space (Z E RNXQ) to high-dimensional data space (X E R RND)) such that Q < D and points close in latent space are nearby in data space, Since the decoder is a non-parametric Gaussian process, the kernel function controls the inductive biascs of the function mapping like smoothncss and periodicity. There typically is no encoder mapping hence, these models are also called Gaussian process decoders. is possible to additionally incorporate a back-constraint or an encoder which maps from the data to latent space putting GPLVMs on the same footing as variational auto-encoders (VAEs; Lawrence & Quiñonero Candela, 2006
Bui & Turner 2015). This amortises the cost of variational inference in very large-scale datasets but we don't employ this setting as it is not straightforwardly applicable to missing data contexts.
This work proposes a novel formulation of the GPLVM based on the idea of a shared latent space. The earlier work by Ek (2009) was the first to propose the idea of a shared data generation process but prechuded truly scalable inference due to the standard ((N3) scaling. We extend this framework two important ways. First, we show that the shared GPLVM is compatible with stochastic variational inference (SVI) (Hoffman et al. 2013) where we derive a joint evidence lower bound which factorises actoss multiple observation spaces due to conditional independence but share predictive strength though inducing locations and latent variables. Secondly, we train the entire model in presence of missing dimensions in one or both of the obscrvation spaces.
We demonstrate this scalable model in an astrophysical application using data of quasars. Quasars are the most luminous galaxies in the universe, powered by accretion onto a central supermassive black hole (SMBH) with millions to billions of solar massos in size. Understanding the formation, growth, and evolution across cosmiic time of quasars and their SMBHs is one of the major goals of observational cosmology today. To this end, precise measurements of the physical properties of quasars are crucial, but they typically demand very expensive and time-intensive observations, multiple epochs are needed to accurately determine, for instance, the quasar's SMBH mass. The high-dimensional data used in this work contains ~ 22, 000 quasars from the Sloan Digital Sky Survey (SDSS: Lyke et al.) 2020: Wu & Shen, 2022) with high-quality spectral information (binned to 590 spectral
pettnssionspiielee along with four scientific labels per quasar, i.e. their black hole mass, luminosity, redshift and so-called Eddington ratio a measure of the quasar's accretion rate. By modelling the spectra and scientific labels through a generative model acting on shared latent space we aim to reason about the physical properties of the quasars just through its "single-epoch" spectral information, thus circumventing the time-intensive multi-epoch observations. Earlier work applying probabilistic generative modelling using a GPLVM to high-dimensional quasar spectra
EEilers et al. 2022
2

---

Published in Transactions on Machine Learning Research (Jan/2025)
have been constrained scalability and examine less than 50 astronomical objects. We demonstrate our framework on datasets over 400x larger. We summarise our key contributions below:
Contributions We propose a probabilistic generalive framework called the Shared slochastic GPLVM which is designed for use cases with multiple outputs/obsent /obervaaieen spaces. We seamlessly account for missing dimensions both at training and test time due to the probabilistic nature of the model. We demonstrate through astrophysical experiments that we can reconstruct previously unseen spectral pixels to a high degree of fidelity, interpolate missing or unobserved spectral regions and predict scientific labels. Crucially, we demonstrate that it possible share predictive strength by learning a common latent variable space Z across multiplo-outputs (X,Y) where Y RNx is additional view of the data, refer to this additional observation space with L dimensions. In this way we indirectly model the relationships and correlation structure between the different observation spaces. We demonstrate this concrctely with an pererice/predictt experiment where we generate/ predict all the scientific attributes at test-time by using latent variables (Z) only informed by the quasar spectra (X), we can denote this cross moosl prediction as X + Z - Y. To the best of our knowledge, predicting multiple outputs with stochastic variational GPs for scalability is methodologically novel. Further, this is the first demonstration of a scalable probabilistic latent variable model in astrophysical settings. In Section 10 we discuss how the development of these unsupervised learning frameworks can instigate novel insights into some of the major open questions in astronomy today,
2 Stochastic Variational GPLVM with Shared latent space
In this section, we first describe background on the stochastic variational GPLVM Lalchand et al. 2022). We then develop the idea of shared latent space and inducing points within the stochastic variational GPLVM framework. The fundamental contribution of this work is develop inference scheme to show that shared latent space does not preclude scalable inference through SVI.
2.1 SV-GPLVM: Stochastic Variational GPLVM
In the traditional formulation underlying GPLVMs we have a training set comprising of N D-dimensional real-valued observations X = ( E RND. These data are associated with N Q-dimensional latent variables, Z = { 1 NxQ where Q < D provides dimensionality reduction (Lawrence) 2004).. The
forward mapping (Z - X) is governed by GPs independently defined across dimensions D. The sparse GP
formulation describing the data is as follows:
p(Z) II N( 0, Il),,
D N(fa Knm Kmm uni Qun);
U( 2,0) II
(1)
d=1
N
D
P(XF,2 = II
II
N (Z fa(zn),
n=1d d=11
where us E R and va is the dth column of X. Knn is the N x N covariance matrix corresponding to a user-chosen positivo-defnite kernel function ko(z, 2²) evaluated at latent points { and parametrised by shared hyperparameters 0. Inducing variables per dimension {ur)2 are distributed with GP prior usl ~ Nu(,, 0, Kmm) (where K Kmm is M x M) computed on inducing input locations .2,...EM Z E R M which live in latent space with Z and have dimensionality Q (matching 2,). In addition, Kum is the N x M cross-covariance computed on the latents z, and the inducing locations žm-
(uglo where Qmn Knn Knm K-1 mn NxN, F= (fal where fa UU
The core idea of Stochastic Variational Inference (SVI) applied to sparse variational Gaussian processes (GPs), as introduced by Hensman al. (2013), that the inducing variables us can be variationally marginalized. This is achieved by learning their variational distribution q(ua) NMmm, using stochastic gradient.
3

---

Published in Transactions on Machine Learning Research (Jan/2025)
methods. Crucially, this approach retains uncollapsed representation of ua throughout the optimization process. While Hensman et al. (013)) proposed SVT for GP regression, Lalchand et al. (2022) extended this
work to GPLVMs where the inputs Z to the GPs are unobserved and each dimension of the high-dimensional output space za is modelled by an independent GP fa with shared kernel hyperparameters.
Succinctly, posterior inference entails minimising the KI-divergence wetween the posterior over unknowns p(F, U; X and the variational approximation q Z). Following Lalchand et al. (2022) when the variational approximation admits the factorisation below in Titsias 2009
(F.UU.) (
(F(FFF)))(WW))( ()a(U)
II d=1
p(fa ua. ()a( II (zx) a=1
(2)
the evidence lower bound (ELBO) can be derived as,
logp(X) 2 (log (log((X]E, -K((@ (LL((V))I)(()
(3)
where q(.) denotes the variational approximation as in Eq. 12) For brevity, we suppress the conditioning over the inducing imputs in the prior p(U)) and the kernel hyperparameters in p(FlU, Z).
If we choose to optimize the latent variables Z as point estimates rather than variationally integrate them out (basically we do not introduce (()) we end up with the following simplification,
p(X) x(x) P(FFU,Z))(U) ()) log P( ()p(U)p(() dFdU L.
(4)
q(U)
Re-writing the lower bound as sum of terms across data points N and outputs/ciimensions D we get,
Dx Σ (logp(zn Zn, oz)) pisala n.d
Σ (Tug(aa)ll((uu1))) + Σ logp(zz)
(5)
Latent point estimates { J==1 can be learnt along with and variational parameters (Z, ma Sa) by taking gradients of the ELBO in Eq. (). However, important constraint is that this formulation assumes single-kernel matrix (single set of kernel hyperparameters) underlying all the independent GPs D. In the next section, we introduce the concept of an additional observation space with dimensions L, which can be modelled using an independent stack of Gaussian processes (GPs), ft. These GPs are equipped to learn their own set of hyperparameters for added ficxibility while sharing the latent embedding Z and inducing inputs Ž capture correlations across the different output spaces.
2.2 Shared joint variational lower bound
In the astrophysical application we focus on in this work we have two observation spaces corresponding N quasars. We denote the quasar spectra (pixels) with the matrix X Nxb and the scientific labels corresponding to the N objects with Y NxL The GPLVM construction models each column (pixel dimension and label dimension) with independent GP, with the GPs corresponding to the pixel dimensions (a) and label dinensions ((t)+ modelled with their own independent kernels and kernel hyperparameters, O₂ and fy. Within each observation space the kernel hyperparameters are shared, so we learn two sets of hyperparameters corresponding to two observation spaces.
fa ~ GP(0P
f GimP((O ko,)
(6)
The priors over the function values are given by,
p(falOz) No k(4)) p(f( NN k!!))
(7)
(8)
4

---

Published in Transactions on Machine Learning Research (Jan/2025)
where K/m and 2 denote the N x N kernel mitrices which rely on their own set of hyperparameters. The two observation spaces also yield two data likelihoods given by,
N
N D
p(X|fnnn 2) = II
p a Sa: 2) = II II
((nndi Ja(zn), ,
(9)
n=1
2 d=1
D+L
p(YIfD+LD+L.2)
I ply
II
II
Ni o((z)),
(10)
1- 1711
In the absence of sparsity the log-marginal likelihood of the joint model compartmentalises nicely due to the assumed factorisation in the likelihoods. We marginalise out the latent function values S1.D and JD+1.D+L per dimension,
p(X, 110 0y 2) = S /
(11)
L.D D+1.D+E
F(X ()p dfIID
(12)
D+1:D+L
1.D
D+L II (=D+1
D+L
plaa Z) II
plun Z)
II
K + +o
II D=++1
N( nn + og)
(13)
d=1
where æa and y! denote single likelihood objective is then given by the following,
of the observation spaces X and Y. The log marginal
D
D+L Σ 1= D7+1
log (( By 2) = Σ log p( Z)+ +
log p(yn 1y.
(14)
In order to induce sparsity we introduce inducing variables ut and u, for each individual dimension in the observation spaces; however, they are underpinned by shared inducing inputs Z which live in the shared latent space Z and share the same dimensionality, Q. With sparse GPs each of the terms in the decomposition above can be bounded by Lz, while the inducing points Z can be shared between the terms yielding the joint evidence lower bound.
log p(X, Y00,, p(X,Y10,,
0y 2) 2 Σ (logp(en Zu; a o3))p(tauus Σ
(L(q(ua)lp(uua)
(15)
n,d
2o Zn, ofP
Σ (L(g(u)) + logp(2)
n.l
Overall we optimisc the shared variational lower bound w.r.t kernel hyperparameters for the two groups of GPs, Og and 0y, variational parameters {ma, Sa) and (m,, and single set of shared latent point estimates Z and inducing inputs 2. We include the full training algorithm in Algorithm
3 Predictions & Reconstructions
High-dimensional points can arrive in different formats for the test data, cither we observe both modalitics (a",y") where x* = " "t] and similarly for y" or only one of the modalities with the other missing, i.e. {x*} or {8*} only. The prediction exercise then involves inferring the latent x" corresponding to the unseen test point.
Since the GPLVM is decoder only model, we cannot obtain the latent embedding x" deterministically; instead we re-optimize the ELBO with the additional test data point (a")"") while keeping all the global and model hyperparameters frozen at their trained values. Note that since the ELBO factorises across data N+1 D+L
points, L((#n, yn 1 x²,y") = Σ Σ
Cn,,, the gradients to derive the new latent point x" only depend
5

---

Published in Transactions on Machine Learning Research (Jan/2025)
0
fa
D
fi
Yn
N
By
o
Figure 1: The graphical model of the shared GPLVM with two sets of independent GPs and their respective hyperparameter sets.
on the respective component terms connected to the data point. In the event of a missing dimension d" for a new data point x", the augmented ELBO for the spectra dimensions can written as follows:
E Cu + logp(z")+ dzfe-
+
Σ (p(=z) Ifa, o2))(a))
(16)
Once we infer z" we can compute the full reconstruction distributions (we may be interested in this the test data point y² had any missing dimensions, for example y [vi v2. ?,.............................. ........77] ,,)""))"))))) which are just the GP posterior predictive for each column or dimension, without loss of generality, for dimension yf:
/ plyit [i, un, x" )(f( uilyrdd
p(uil="
/ plui Ii,, un, =" )o
(17)
f plyi 20,,"11"000 )))(z))))
where p(yilu,2 2") WKKmKKon Un K.. Kanknn Km. + og) and g(ui) N(m St) refers to the optimised variational distribution. The final integral is tractable and gives the following form:
p(vil=*) N(Km mi, Kmnknt - Kmm) Km= + og)
(18)
and similarly for any 2x- For cross-modal reconstruction (where we only observe one modality of the test. data) the latent z" acts the information bottleneck, hence, the same posterior prodictive distributions be derived, xt + -² +
VI = I,... D+ L.
4 Schematic of the Model
In Fig. we present schematic of the model architecture with two observation spaces (X,Y), the corre- sponding stacks of individual GPs {fa} and {fi} which model the individual columns of the spectra X and scientific attributes Y, respectively, and the low-dimensional latent space Z. The dimensionality of the latent and observation spaces is denoted by Q, D, L, respectively, and N denotes the nunber of objects / data points (quasars). Note that the correlation between the two observation spaces is not explicitly but implicitly modelled through a shared latent space. Generating a single data point (æn, yn (a row across and Y) entails forward pass through the GPs, where In [..., Ind. generated as [fi (z,), f.......ff(.,,,,, and yn [--., (.-., ... is generated [sD+1 (zn), [D+2(Zn), ... f (2.))

---

Published in Transactions on Machine Learning Research (Jan/2025)
Sing'e column of spectro ccross dbjects
(faler
Single column aeross objects
D
L
Latent space
Observations Data space
Figure 2: Shared GPLVM with multiple observation spaces. The blocks on the right-hand side denote the double observation spaces (X, Y) of quasar spectra and scientific labels respectively. the center are two stacks of GPs, one for each observation space which control the data generation process through the shared latent space. In the figure above assume Q (for case of visualisation) since denote the GPs are two dimensional surfaces, however, lypically Q can be higher than 2 corresponding to higher dimensional GPs.
5
Algorithm
We enclose the pseudo-code in Algorithm 2 for stochastic variational inference in the context of the shared model for clarity. Let Cx and Cy denote the ELBO's for each of the observation spaces and let _ and C denote the ELBOs formed with randomly drawn mini-batch of the data (across all dimensions). For mini-batch (subset) of the data XB C X, the mini-batch ELBO is given by,
C2 (B)
N
Σ h.d
[logp(zb.alfa 2s, 2 o3) + Σ log p(z))()))))))))))))))))
Σ KLlg(u))ll(uul2)) (19)
B
where the scaling term is important for the mini-batch ELBO to be estimator of the full-dataset ELBO.
5.1 Computational Cost
The training cost of the canonical stochastic variational GPLVM is dominated by the number of inducing points O(MDD) (free of N) where M & N and D is the data-dimennionalltt (we have D GP mappings fa, one per output dimension). The practical algorithm is made further scalable with the use of mini-batched learning. In our shared model with two sets of GPs, the dynamics of the training cost are the same except that they increase linearly in the number of additional dimensions (L), making the cost (M + L)). The number of global variational parameters to be updated each step (parameters of ((U)) MQ M(D L) + M²(D+ where MQ are the M Q-dimensional inducing inputs Z (shared), MDD+ + L) is the size of the mean parameters of the inducing variables uad u, and M L) the full-rank covariances of the inducing variables. The local variational parameters Z (the latent embedding shared across GPs) are of size NQ and model hyperparameters (kernel hyperparametent are of size 2Q + 4, which account for Q input lengthscales, scalar signal variance and noise variance per GP group {fa} and {ffi. We use the squared exponential kernel with automatic relevance determination across both sets GPs.
6 Related methodological work
In this section we present related work in multi-output Gaussian processes in more detail. Canonical multi-output Gaussian processes almost always refer to the supervised framework or multi-output regression

---

Published in Transactions on Machine Learning Research (Jan/2025)
Algorithm 1: Training Framework
Input: ELBO objective -= Lz + Ly: gradlent based optimiser optim(), observation spaces X (spectra) and Y (scientific
labels)
Initial model params:
0 (0x, By) (covariance hyperparameters for GP mappings fDD.ID++
2² (og.0] (variance of the noise model for each likelihood),
(=n)==l (point estimates latent embeddingg)
Initial variational params:
E R XM (inducing locations),
(m,, Ss. L
(global variational params for inducing variahles per dimension 1t))
while not converged do
Choose random mini-batch of the data from both the observation spaces XB C C,YDCC
Form mini-batch estimate of the ELBO: _ + c
Gradient step for global parameters
2,,02,2,1 (8) 1/ Vz (C2 (B) + c (where Zz are the latent embeddings corresponding to
Gradient step for local parameters 1
points in the niti-bbtch)
Update all parameters 2,0,00 02 2B [= < optim() using gradients g.
end
return 0,02, 2,
Algorithm 2: Prediction Framework
(Predict z" corresponding to test X*,Y*)
Input: Trained global and local parameters 0, 0, 2, 2, A, Z, test observation spaces X (spectra) and Y* (sclentific labels). Initialise latent embedding = (Ep hx corresponding the N additional data points.
points.
Extend the joint, ELBO include
correspanding
C t- Lx End n. flog plone aa o pa (uy. 21((aa) + Σ.. logp(=n*)
C C > n (logplyn' .lff, En,, , (u + >.. logp(zn")
C +
Freeze all global and local paramcters except for Z* while not converyed do Gradient step for Z": 1" <- Vz- L"
Update 2 + optim() using gradients 1"
end return 2
(Note that the gradients of La and Ly with respect to z" are 0 and the only terms that are optimised are the additional terms corresponding to the new data points.)
Álvarez et al. 2010) where the targets {2a(x))2 multidimensional. continnous and Gaussian distributed corresponding to inputs z E RR. The main focus is on defining suitable cross-covariance function between the outputs, this allows treatment of the outputs as single GP with a suitable covariance function (Alvarez
itt al. 2012). The intrinsic and linear co-regionalization models (Goovaerts) 1993 Journel & Huijbregts/ 1976
are & popular approach in which each output is modelled as weighted sum of shared latent. functions (GP) where typically the number of latent processes is smaller than the number of outputs enabling efficiencics. To
some extent, multitask learning with GPs can be viewed as an instance of multi-oulput learning where we want to avoid tabula rasa learning for each task and evolve framework for sharing information between multiple tasks [Bonilla et al. 2007p. Further, there is the simplistic paradigm where each output is modelled with its own independent single-output GP and no correlation between the outputs is assumed. This approach, while easy to implement, is severely limited in its ability to jointly model the outputs.
In the unsupervisod paradigm, the starting point a high-dimensional data matrix N D. The conventional Gaussian process latent variable model (GPLVM) operates like a multi-output model by default where each column of the data is modelled by independent GP the same shared set of inputs, kernel function
8

---

Published in Transactions on Machine Learning Research (Jan/2025)
Metrics
Mean Absolute Error
Models
Baseline
Shared (ours)
Specialised (ours)
Attributes (4)
Spectra
0.1103 0.0035
0.0972 0.0024
0.0969 = 0.0013
Bolometric Luminosity
0.2444 + 0.0043
0.2144 + 0.0033
0.2069 + 0.0014
Black hole mass
0.2712 + 0.0031
0.2319 + 0.0028
0.2362 + 0.0025
Eddington Ratio
0.2525 + 0.0070
0.2118 + 0.006
0.2226 = 0.0015
Table 1: Summary of test-time reconstruction abilities. Mean absolute error on denormalised data (# standard error of mcan) evaluated on average of 5 splits with 75% of the data used for training in the 22k dataset. The shared and specialised versions of the model outperform the baseline model in all the reconstruction tasks. In the shared vs. specialised comparison the performance difference is not statistically significant for the spectra roconstruction but overall favours the shared model with hyperparameter sharing across the scientific label reconstruction.
and hyperparameters. The GPLVM is decoder-only model, and some of their prominent variants are the constronnt GPLVM (Lawrence & Quiñonero Candela 2006), the supervised GPLVM (Jiang et al.)
2012), the discriminative GPINM (Urtasun & Darrell, 2007] and the shared GPIVM (Ek 2009 The latter considers the task of dealing with multiple views or observation spaces (each of which high-dimensional), This work is built on the idea of a shared latent space underlying the multiple observation spaces similar Ek et al. (007)) but adapts it for scalable inference using stochastic variational inference (SVI). The stochastic
variational GPLVM proposed in Lalchand et al. (022), but only considered single obscrvation space in canonical GPLVM.
7 Experiments
In this section, demonstrate range of experiments aimed at assessing the reconstruction quality of test (unseen) quasar spectra and scientific attributes, well as the robustness of uncertainty quantification by computing the negative log predictive density (NLPD) of test labels - -o=')("" under the predictive distribution where z" has been informed by both modalities (spectra and labels) and only spectra.
The data used in this work are quasar spectra observed as part of the Sloan Digital Sky Survey (SDSS) DR16 (Lyke et al. 2020). We chose all quasars with spectra that have signal-to-noise ratio (SNR) per pixel > 10, which results in total of 22844 quasar spectra. The observed spectra are shifted into the rest frame wavelength space, re-binned onto common wavelength grid, and flux normalized to unity at around 2500 Å. We mask strong absorption lines that might arisc in the spectra due to foreground galaxies along our line-of-sight to the quasar, and are thus not intrinsic spectral features of the quasar. The four scientific labels for these quasars are (1) their SMBH mass, (2) their bolometric luminosity, i.e. the total power output across electromagnetic wavelengths, (3) their redshift which denotes the factor by which the emitted wavelengths have been "stretched" due the expansion of the universe, and (4) their Eddington ratio, which is a measure of the accretion and growth rate of the SMBH. All measurements were previously uniformly determined by Wu & Shen 2022). We conduct experiments across two data sets with 1k (sce Table and 22k points (as ablation to assess performance with a smaller dataset). The "Baseline" model in the experiments refers to the canonical stochastic variational GPLVM (Lalchand et al.) 2022) which treats multiple observation spaces using the same set of independent GPs learning a single set of kernel hyperparametors. The 'Specialised' model refers to the shared model but with individual hyperparameters for each scientific label.
7.1 Reconstructing Quasar spectra
We assess the quality of our probabilistic generative model in reconstructing test quasar spectra. At test. time we deal with spectra and scientific labels from test quasars denoted by Ya) (we use the index 'gt' to denote ground truth). The 2-stop prediction Iearns the low-dimensional shared latent variables Zt (point estimate per data point) which acts as input to the GP decoder predicting the corresponding spectra 7' X eat- Note that the ground truth spectra contains several missing pixels (dimensions) and
9

---

Published in Transactions on Machine Learning Research (Jan/2025)
the probabilistic decoder provides reasonable reconstruction at those locations. In Fig. g we visualise the reconstruction (posterior predictive mean) of four test quasars along with ground truth measurements and 95% prediction intervals. We achieve remarkably good reconstruction estimates from the latents; further, the prediction intervals capture the ground truth spectra providing robust coverage at peaks and extrapolated regions.
7.1.1 Reconstructing missing spectra
In this experiment, we test the generative model's ability to learn from massively missing chunks of the spectra at test-time. We observe partial window of the spectra in each plot (given by the shaded region in Fig. 4 hence the latent variables corresponding to these points are only informed by the observed region. We then reconstruct the whole spectrum from the latent variables informed by the partial spectra. We enclose our results in Fig. A Reconstruction entails the following inference steps: X partial + 2" + Xhull We note that the quality of the mcan prediction deteriorates compared to the fully observed test point predictions. However, the coverage of the prediction intervals robust even as we move away from the shaded observed regions. Furthermore, if the model is given a spectral region with very little information (e.g., in the bottom panel of Fig. 4 the shaded region contains information only about the quasar's continuum), the uncertainties increase significantly, as one would expect.
7.2 Predicting scientific labels only from spectra X*
The dimensions I. corresponding to the scientific labels in the dataset governed by their own GP decoders [fi))441 are a critical prediction quantity. The ability to reconstruct these quantities from learned latent variables an important test of the generalizability of our model. Very often astronomers want reason about the scientific altributes of quasars just by analysing their spectra. In this experiment, we demonstrate precisely this use case in which latent variables z are informed only by spectra X*, computing the cross- modal prediction involves learning Z" from the ground truth spectra and then using just Z" to predict the scientific labels Y", succinctly, we can write these steps as: X + Z +Y*
In Fig. 5 we demonstrate the accuracy of our reconstructions by plotting each of the dimensions against ground truth held-out data. We show reconstructions for 200 test points randomly sampled from the full test set. Each point in the scatter denotes quasar, and the x-axis denotes the ground truth measurement. The orange vertical error bars denote the intervals lo computed by extracting the diagonals from the predictive posterior GP for cach dimension. We can observe high degrec of prediction accuracy across the three scientific labels and, furthermore, the reconstruction quality is weakly correlated to the spectral signal-to-noise ratio (SNR) (at least beyond the data quality cut of SNR > 10 which we apply as a preprocessing step); this is mainly due to acute data imbalances where more than 90% of the objects appear at the lower SNRs affecting prediction quality and generalization for the quasars at higher SNRs (see Appendix D). Note that we did not attempt to reconstruct the redshift label from the spectra X. This because the spectra have already been normalized account for their redshift. Specifically, the spectra were shifted to rest-frame wavelength space by dividing the observed wavelengths by 1 + 2 (where is the redshift). As a result, normalized spectral shapes no longer contain significant information about their redshift (as shown in Yang et al. 2021] and
Onorato et al. 20244) ).
However, we still include the redshift information in our multimodal GPLVM. This is because excluding it marginally weakens the results, both for spectra reconstruction and label prediction. We hypothesize that, while the spectral shapes themselves do not carry redshift information, the patterns of missing pixels may still encode some redshift-related information.
From Table 2 we can see that the variability of NLPD estimates is much higher under the cross-moda protocol of decoding scientific labels on the basis of the spectrum. In other words, reconstructing scientific labels from spectra alone has a much higher uncertainty than reconstructing scientific labels from both modalities. As a result, the shared model offers only a subtle improvement over the baseline model in terms of uncertainty quantification. In the fully observed protocol, the shared model more confident (Table 2 and accurate (Table D in its predictions relative to the baseline model. It is also important to note that under the much harder modal protocol, the performance of the shared model surpasses that of the baseline model.
10

---

Published in Transactions on Machine Learning Research (Jan/2025)
Spectra reconstruction
10
predicted spectrum data; log10 (M./Ma) 9.1045
normalised flux
+1.96a
4-
2-
1250
1500
1750
2000
2250
2500
2750
3000
predicted spectrum
data; log10 (M./Mg) 9.4316
normalised flux
1-
-1
1250
1500
1750
2000
2250
2500
2750
3000
predicted spectrum
data; log1o (M./Ma) = 9.4158
normalised flux
11.960
1250
1500
1750
2000
2250
2500
2750
3000
10
predicted spectrum
data; logı0 (M./Ma) = 8.9762
#1.96a
normalised flux
1250
1500
1750
2000
2250
2500
2750
3000
rest-frame wavelength [A]
Figure 3: Reconstruction plots of Lest quasar spectra with +1.96a intervals. The blue curve denotes the posterior predictive mean at each dimension.
11

---

Published in Transactions on Machine Learning Research (Jan/2025)
Spectra Reconstruction [from partial observations]
predicted spectrum
data: logio (M./Ma) = 8.9762
Observed Region
normalised flux
11.960
1250
1500
1750
2000
2250
2500
2750
3000
normalised
1250
1500
1750
2000
2250
2500
2750
3000
normalised flux
1250
1500
1750
2000
2250
2500
2750
3000
7-
normalised
4-
3-
2-
1250
2000
2750
3000
rest-frame wavelength [A]
Figure 4: Reconstruction of a single spectra from the latent informed by partially observed spectrum. The shaded orange regions demote the "observed" wavelength regions for this experiment. Note in the 4th panel the 95% prediction intervals are wider as they were observed over a shorter and less informative wavelength window.
12

---

Published in Transactions on Machine Learning Research (Jan/2025)
08.00
Buluemeds luminadiy
10.35
Black hale muse
Lddington Patio
47:15
00.00
-25
23
Prodicedd 1eirr
07.50
Predictor logo sggehau
07.49
800-
07.00
46.75
Pre.dicted
9:25
28
+8.50
0.75
-1.50
-5
4.50
-1.75
46.0
40.5
41.0
47.5
40.0
0.30
8.20
5.75 10.00
20.25
-2.0
-1.
-1.0
-0.5
0.0
Messured lode (erpor
Messured (M.
Nessured 16022
Figure 5: Scientific label prediction for the quasars' bolometric luminosity (left), black hole mass (middle) and Eddington ratio (right) colored by the SNR of their spectra based on test X only. The dashed black line (---) denotes the 1-to-1 line to aid visualisation of reconstruction accuracy. The vertical and horizontal errorbars (-) denotes posterior predictive standard deviation.
Experiment
Luminosity
Black hole mass
Eddington Ratio
Baseline (Fully observed)
0.076 + 0.0012
0.148 = 0.0043
0.1178 + 0.0031
Baseline (Spectra observed)
0.1143 + 0.042
0.1742 + 0.0621
0.1376 + 0.0646
Shared (Fully observed)
-0.029 + 0.0018
0.1009 0.0013
0.0702 + 0.018
Shared (Spectra observed)
0.0862 + 0.051
0.1604 + 0.0245
0.1294 + 0.0461
Table Summary of test-time uncertainty quantification under the full and partial reconstruction framework. Mean negative log predictive density (lower is better) on de-normalised data across train/teer splits. A higher NLPD indicates lower confidence in the predictions and the model indicates this when reconstructing the labels of quasar just the basis of its spectrum. The fully observed model corresponds to predictions based on test data when the all modalities are observed and reconstructed, (X gt Yga) + Z + Y or when only the spectrum of the quasar is informing the latent, (X* gt + Z" - Yoa) The latents are in turn used to decode the scientific labels.
7.3 Generating spectra corresponding to synthesized labels: an ablation study
In this experiment, we demonstrate the ability of the model to generate spectra corresponding to synthesized scientific labels. Concretely, we simulate artificial labels by systematically varying only one of the labels within a reasonable range each plot in Fig. G The range of variation for each label is summarized by the colorbar in each plot, for instance, in the black hole mass simulation we generate 100 spectra (X*) corresponding simulated labels (Y*) of black hole masscs log (M. in the range [7.9, 10.0]; in order to ablate the influence of other labels (redshift, bolometric, luminosity, and the Eddington ratio) on the spectra we keep their values fixed to mean values computed from the training dataset. Lastly, we can summarize the inference steps as the inverse of the previous section: Y* + 2 + X*. The ability for prediction and generation of cross- modalliolsss is an important strength of our design.
Reassuringly, the model exhibits the expected behaviour. For instance, the emission lines are broader for quasars with higher black hole masses, and the spectral dependence with bolometric luminosity shows e.g. in the Mg, II line at 2800Å the well-known Baldwin effect (Baldwin 1977), which indicates that quasar spectra show decreasing equivalent width of their UV and optical emission lines with increasing bolometric luminosity.
addition, the gray bars denote the 95% prediction intervals averaged across the 100 spectra at each dimension. The uncortainty intervals are wider at higher wavclengths, as there a high concentration of missing pixels in the training data at those wavelengths; hence, there is greater uncertainty about the generated spectra.
13

---

Published in Transactions on Machine Learning Research (Jan/2025)
Spectra generated for simulated Black hole mass
10
Nomalize Flux
9.2
9.0
1250
1500
1750
2250 Restframe wavelength [A]
2000
2500
2750
3000
Spectra generated
luminosity
20
47.6
15
47.4
Juminosyy
Normalised Flux
17.2
46.8
46.6
0-
46.4
1250
1500
1750
2000
2250 Restframe wavelength [A]
2500
2750
3000
Spectra generated for simulated Eddington ratio
0.25
25
0.00
Normalsed
Eddimgtom
11.25
-1.50
1250
1500
1750
2250 Restframe wavelength [A]
2000
2500
2750
3000
Figure An experiment demonstrating cross-modal prediction Y* - ---+"* - In the plots we show generated quasar spectra for simulated labels for black hole mass, bolometric luminosity and Eddington ratio. In each plot we vary the respective scientific label in reasonable range (shown by the range the colorbar) while keeping the other labels to fixed values.
+
14

---

Published in Transactions on Machine Learning Research (Jan/2025)
7.4 On benchmarking with traditional regression approach
The dependency between the spectra and labels is modelled through the shared latent space, which generates both views of the data. Without a generative model it would not be possible to generate spectra corresponding to scientific labels (generate X corresponding to Y), we would only be able to predict labels based on fixed observed spectra. With a joint generative model, it is possible to generate both spectra and labels or generate one conditioned only on observing the other view through a jointly learned compression. In other words, we can extract the shared latent (Z) corresponding to X (spectra) to predict reconstructed estimate of the spectra and the labels (X, Y), or the other way round Y optimise, Z decade, (X, Y) as shown in experiments. Further, the spectra are only partially observed with several missing pixels making canonical regression less missing/uncertain straightforward as one would need to account for missing/ uncertain inputs. In the case of training GP to
predict labels directly from fixed high-dimensional spectra, one would have to impute the missing pixels at the outset, as the prior covariance matrix cannot be computed on missing inputs.
8 Limitations
The primary limitation of the model arises from the fact that GPLVMs are decoder-only models, they constitute a (potentially) smooth mapping from the latent space to the data space. This means that points close in latent space will be close in data space but not. the other way round. Data space similarities can be preserved by including a back-constraint or an encoder that additionally maps the data to latent space [Lawrence & Quiñionero Candcla, 2006 Bui & Turner 2015). The encoder usually takes the form of a neural network, but other choices are possible. As the title of the manuscripl emphasizes, the model we propose is decoder-only model. The absence of encoder complicates the test-time inference as there is no deterministic way to access the latent points z" corresponding to the new unseen observation y")). Inferring the latent point corresponding to the unseen test point entails freezing the model and variational parameters post-training and re-optimizing the ELBO objective subject to (=*) with the new data point(s) included (see Algorithm 2
This optimization procedure for test inference ends up being too inconvenient in contrast to encoder- decoder models such as VAFs (Kingma & Welling 2013), where inferring a latent point corresponding to unseen data point entails a forward pass through the trained encoder network (constant-tins predictions O(1)). The main reason an extension to autoencoded shared stochastic GPLVM is not straightforward is due to the presence of missing data. The encoder network must be capable of handling arbitrarily missing dimensions. Mothods like the partial VAE (Ma et al., 2018a b) address this challenge in the context of VAEs where each observation is augmented with a column index indicating the observed dimension; the encoder network then processes these tuples as a 'set' with a permutation invariant set encoding function. A similar approach can be straightforwardly adapted from our shared latent space setup with Gaussian process decoders. The combination of a set encoder (to process missing dimensions) and non-parametric decoder has not appeared in the literature to the best of our knowledge. We are currently working on incorporating this feature into our framework. Note that the presence of an encoder acts as a constraint in the model, and there is an inherent trade-off in terms of faster test inference and marginally weaker reconstructions predictions.
9 Methodological extensions for future work
A natural extension of the shared GPLVM proposed here is to incorporate encoder, bringing the whole architecture more in line with modern boop-generative autoencoder models while simultancously preserving the advantage of uncertainty quantification in the outputs. Bui & Turner (2015) propose the canonical GPLVM with an encoder and train it with SVI, therefore, a natural extension is to demonstrate it for the shared setting. The shared setting in the data context presented here opens up questions about the encoding of partially observed vectors along with fully observed ones, and masked encoders (He et al. 2022) might be relevant. architectural paradigm to explore.
15

---

Published in Transactions on Machine Learning Research (Jan/2025)
Another extension is related to the structure of the kernels, since the kernels factorise over dimensions and are closed over the multiplicative operation, one approach is to select composite product kernel over dimensions of the latent space instead of two sets of GPs with their own respective kernels. The former approach (not studied here) induces a partitioning of the latent space into dimensions which would cater to specific views of the data and can be modelled with individual kernel functions. In our context, we did not see the gain the mixture kernel formulation; note that this approach raises an additional question as to how many latent dimensions to allocate each kernel function. Nevertheless, partitioning the latent space and modelling with mixed composite kernels could be beneficial in specific settings. More experiments are required to understand the exact settings in which this approach may outshine the alternative framework presented here.
10 Scientific Interpretation and Significance
Our new generative model allows to simultaneously model the spectral properties of quasars as well as their scientific labels, thus opening up novel possibilities to study the evolution of quasars and their dependency on physical parameters governing the SMBH growth. In the following, we will highlight just two possible and exciting future applications of this work.
Astronomers observe SMBHs with billions of solar masses size in the center of very distant, high-redshift. quasars, at a time when the universe still in its infancy and only a few hundred Myr (million years) old. This rapid growth of SMBHs in the very short amounts of available cosmic time has been open puzzle for decades, and it has been argued that very high accretion rates in excess of the theoretical upper limit, the so-called Eddington limit, are required to explain the rapid black hole growth. However, obtaining precise black hole mass measurements of quasars is challenging and time-intensive as it requires multi-epoch observations of certain emission lines in the quasar spectra (Peterson 1993 Barth et al. 2013). This procedure becomes increasingly challenging for very distant, high-redshift. quasars, as relativistic time-dilation effects require longer timespans of these observations. Furthermore, the traditionally used rest-frame optical emission lines to calibrate the black hole masses unobscrvable with ground- based observatorics, these optical wavelengths have been shifted to the infrared at larger distances, which require space-based telescopes, such as NASA's 's recently launched James Webb Space Telescope. Using the new generative model provides us with the possibility to determine the masses of SMBHs for quasars at all redshifts from single-epoch data observed with ground-based observatories, as it does not require wavelength coverage of specific emission lines; further, it can handle missing data.
Furthermore, we have shown that our model is able to predict other physical properties of quasars such their bolometric luminosities (see Fig. 5). This suggests that we obtain a measurement of the quasars' absolute luminosities using their spectral information alone, which provides new opportunity to quasars as the so-called "standard candles." Standard candles are incredibly valuable for astronomy, as knowing the luminosity of an object allows one to determine its distance. Famously, supernovae have been used as slandard candles, which led to the Nobel Prize winning discovery of the expansion of our universe and the existence of dark energy (Riess et al. [998). The use of quasars as standard candles has previously been suggested leveraging the relationship between the X-ray and UV luminosities of quasars to determine their distances (e.g- Tusso & Risaliti 2017 Sacchi et al., 2022). Our generative model enables predictions of the quasars' bolometric luminosities leveraging these spectral dependencies on luminosity, enabling the use of quasars standard candles (or "standardizable candles") to map the expansion of the universe to larger distances than possible with supernovae. Further testing of the required precision with which the quasars' luminosities and thus distances can be predicted by the generative model will be necessary in order to determine whether the constraints on the universe's 's expansion rate obtained by the more numerous quasars can improve on the current constraints by the fewer but very accurately measured distances from supernovae.
16

---

Published in Transactions on Machine Learning Research (Jan/2025)
References
Mauricio Álvarez, David Luengo, Michalis Titsias, and Neil D Lawrence. Efficient multioutput gaussian processes through variational inducing kernels. In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, pp. 25-32. JMLR Workshop and Conference Proceedings, 2010.
Mauricio Alvarez, Lorenzo Rosasco, Neil D Lawrence, et al. Kernels for vector-valued functions: review. Foundations and Trendso in Machine Learning, 4(3):195 266, 2012.
Jack A. Baldwin. Laminosity Indicators in the Spectra of Quasi-Stellar Objects. The Astrophysical Journal, 214:679 684, June 1977. doi: 10.1086/155294
Aaron J. Barth, Vardha Bennert, Gabriela Canalizo, Alexei V. Filippenko, Elinor L. Gates, Jenny E. Greene, Weidong Li, Matthew A. Malkan, Anna Pancoast, David J. Sand, Daniel Stern, Tommaso Treu, Jong-Hak Woo, Roberto J. Assef, Hyun-Jin Bae, Brendon J. Brewer, S. Bradley Cenko, Kelsey T. Clubb, Michael C. Cooper, Aleksandar M. Diamond-Stanic, Kyle D. Hiner, Sebastian F. Hönig, Eric Hsiao, Michael T. Kandrashoff, Mariana S. Lazarova, A. M. Nierenberg, Jacob Rex, Jeffrey M. Silverman, Erik J. Tollerud, and Jonelle L. Walsh. The Lick AGN Monitoring Project 2011: Spectroscopic Campaign and Emiission-line Light Curves. The Astrophysical Journal Supplement Series, 217(2))26, April 2015. doi: 10.1088/0005 /0/0000607000/7///////
Edwin V Bonilla, Kian Chai, and Christopher Williams. Multi-task gaussian process prediction. Advances in neural information processing systems, 20, 2007.
Thang D. Bui and Richard E. Turner. Stochastic variational inference for Gaussian process latent variable models using back constraints. In Black Bor Learning and Inference NIPS workshop, 2015.
Giuseppe Carleo, Ignacio Cirac, Kyle Cranmer, Laurent Daudet, Maria Schuld, Naftali Tishby, Leslie Vogt- Maranto, and Lenka Zdeborová. Machine learning and the physical sciences. Reviews of Modern Physics, 9144):045002, 2019.
P. Clarke, P. V. Coveney, A. F. Heavens, J. Jäykkä, B. Joachimi, A. Karastergiou, N. Konstantinidis, A. Korn, R. G. Mann, J. D. McEwen, S. de Ridder, Roberts, Scanlon, F. S. Shellard, and A. Yates. Big data in the physical sciences: challenges and opportunities. ATI Scoping Report, 2016.
Adrian Corduncanu and Tommi S Jaakkola. Contimuation methods for mixing heterogenous sources. arXiv preprinl ar orxi::1301. 0562, 2012.
Anna- Caaa-Chiistina Eilers, David W Hogg, Bernhard Schölkopf, Daniel Mackey, Froderick B Davies, and Jan-Torge Schindler. A generative model for quasar spectra. The Astrophysical Journal, 938(1):17, 2022.
Carl Henrik Fk. Shared Gaussian process lalent variable models. PhD thesis, Citeseer, 2009.
Carl Henrik Ek, Philip H. S. Torr, and Neil D. Lawrence. Gaussian process latent variable models for human pose estimation. In International workshop on machine learning for multimodal interaction, pp. 132-143. Springer, 2007.
Elena Facco, Maria d'E d'Errico, Alex Rodriguez, and Alessandro Laio. Estimating the intrinsic dimension of datasets by a minimal neighborhood information. Scientific reports, 7(1):12140, 2017.
P Goovaerts. Spatial orthogonality the principal components computed from coregionalized variables. Mathematical Geology, 25:281 302, 1993.
Kaiming He, Xinlei Chen, Saining Xie, Yanghao Li, Piotr Dollár, and Ross Girshick. Masked autoencoders are scalable vision learners. In Proceedings of the IEEE/CVE conference on computer vision and pottern recognition, pp. 16000 16009, 2022.
James Hensman, Nicoló Fusi, and Neil D. Lawrence. Gaussian processes for big data. Proceedings of the T Nindh Conference Uncertoinly in Artificial Intelligence (UAT2013), 2013.
17

---

Published in Transactions on Machine Learning Research (Jan/2025)
Matthew D. Hoffman, David M. Blei, Chong Wang, and John Paisley. Stochastic variational inference. Journal of Machine Learning Research, 14(4): 100111033-3333 1347, 2013. URL http:://hhaal..ooaaapppt/t/t//// html
Xinwei Jiang, Junbin Gao, Tianjiang Wang, and Lihong Zheng, Supervised latent linear gaussian process latent variable model for dimensionality reduction. IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cyberneties), 42(6):1620 2012.
Andre G Journel and Charles Huijbregts. Mining geostalisties. 1976.
Diederik P Kingma and Max Welling. Auto-encoding variational Bayes. In International Conference on Learning Represenlations, 2013.
Vidhi Lalchand, Aditya Ravuri, and Neil D. Lawrence. Generalised GPLVM with Stochastic Variational Inference. In Gustau Camps- Valls, Francisco J. R. Ruiz, and Isabel Valera (eds.), Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, volume 151 of Proceedings of Machine Learning Research, pp. 7841 7864. PMLR, 28 30 Mar 2022. URL https: ptttps://proceedingg. maapeedddinnmmm.d........mm.meee mlr. p v151/1alchand22 html
Neil D. Lawrence. Gaussian process latent variable models for visualisation of high dimensional data. In Advances in neural information processing systems, pp. 329-336, 2004.
Neil D. Lawrence and Joaquin Quiñonero Candela. Local distance preservation in the GPLVM through back constraints. In Proceedings of the 23rd international conference on Machine learning, pp. 513 520, 2006.
E. Lusso and G. Risaliti. Quasars standard candles. I. The physical relation between disc and coronal emission. Astronomy and Astrophysics, 602:A79, June 2017. doi: 10.1051/000146331/22066363
Brad W. Lyke, Alexandra N. Higley, J. N. McLane, Danielle P. Schurhammer, Adam D. Myers, Ashley J. Ross, Kyle Dawson, Solène Chabanier, Paul Martini, Nicolás G. Busca, Hélion du Mas des Bourboux, Mara Salvato, Alina Streblyanska, Pauline Zarrouk, Etienne Burtin, Scott F. Anderson, Julian Bantista, Dmitry Bizyaev, W. N. Brandt, Jonathan Brinkınann, Joel R. Brownstein, Johan Comparat, Paul Green, Axel la Macorra, Andrea Muñoz Gutiérrez, Jiamin Hou, Jeffrey A. Newman, Nathalie Palanque-Delabrouille Isabelle Pâris, Will J. Percival, Patrick Petitjean, James Rich, Graziano Rossi, Donald P. Schneider, Alexander Smith, M. Vivek, and Benjamin Alan Weaver. The Sloan Digital Sky Survey Quasar Catalog: Sixteenth Data Release. The Astrophysical Journal Supplement Series, 250(1):8, September 2020. doi: 10.3847 008887/1138-4 11/1111888488 13366//bb622)
Chao Ma, Wenbo Gong, José Miguel Hernández-Lobato, Noam Koenigstein, Sebastian Nowozin, and Cheng Zhang. Partial vac for hybrid recommender system. In NIPS Workshop on Bayesian Deep Learning, volume 2018, 2018a.
Chao Ma, Sebastian Tschiatschek, Konstantina Palla, José Miguel Hernández-Lobbato, Sebastian Nowozin, and Cheng Zhang. Eddi: Efficient dynamic discovery of high-value information with partial vae. arXiv preprint Xix: 11142, 2018b.
Leland McInnes, John Healy, and James Melville. UMAP: Uniform manifold approximation and projection for dimension reduction, 2020.
Kevin P Murphy. Machine learning: probabilistic perspective. MIT press, 2012.
Silvia Onorato, Joseph F. Hennawi, Jan- Torge Schindler, Jinyi Yang, Feige Wang, Aaron J. Barth, Eduardo Bañados, Aana-Christinn Eilers, Sarah E. I. Bosman, Frederick B. Davies, Bram P. Venemans, Chiara Mazzucchelli, Silvia Belladitta, Fabio Vito, Emanuele Paolo Farina, Irham T. Andika, Xiaohui Fan, Fabian Walter, Roberto Decarli, Masafusa Onoue, and Riccardo Nanni. Optical and near-infrared spectroscopy of quasars at z > 6.5: public data release and composite spectrum. Xiv e-prints, art. arXiv:2406.07611 June 2024. doi: 10.48550 Xiv 2406.07612.
18

---

Published in Transactions on Machine Learning Research (Jan/2025)
Bradley M. Peterson. Reverberation Mapping of Active Galactic Nuclei. Publications of the Astronomical Society of the Pacific, 105:247, March 1993. doi: 10.1086/133140,
Adam G. Riess, Alexei V. Filippenko, Peter Challis, Alejandro Clocchiatti, Alan Diereks, Peter M. Garnavich, Ron L. Gilliland, Craig J. Hogan, Saurabh Jha, Robert P. Kirshner, B. Leibundgut, M. M. Phillips, David Reiss, Brian P. Schmidt,, Robert A. Schommer, R. Chris Smith, J. Spyromilio, Christopher Stubbs, Nicholas B. Suntzeff, and John Tonry. Observational Evidence from Supernovae for an Accelerating Universe and a Cosmological Constant. The astromomical journal, 116(3):1009-1003 September 1998. doi: (0.109//300499
A. Sacchi, G. Risaliti, M. Signorini, E. Lusso, E. Nardini, G. Bargiacchi, S. Bisogni, F. Civano, M. Elvis, G. Fabbiano, R. Gilli, B. Trefoloni, and C. Vignali. Quasars high-redshift standard candles. Astronomy and Astrophysics, 663:L7, July 2022. doi: 10.1051/0004-636) (360000446661/100 / /202243111
Michalis Titsias. Variational learning of inducing variables in sparse Gaussian processes. In Artificial Intelligence and Statistics, pp. 567-574 667,574, 2009.
Raquel Urtasun and Trevor Darrell. Discriminative Gaussian process latent variable model for classification. In Proceedings of the 24th international conference on Machine learning, pp. 927-934, 2007.
Laurens van der Maaten. Learning parametric embedding by preserving local structure. In David van Dyk and Max Welling (eds.), Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics, volume of Proceedings of Machine Learning Research, pp. 384-391, Hilton Clearwater Beach Resort, Clearwater Beach, Florida USA, 16-18 Apr 2009. PMLR. URL https://proceedngg..t press//rsmmaatee00ahhht
Laurens van der Maaten and Geoffrey Hinton. Visualizing data using t-SNE. Journal of Machine Learning Research, 9(86):2579-2600 2008. URL /ttp/////////mm.roor////////aae http:///w.ao.rreaooooo///////aaa
html
Qiaoya Wu and Yue Shen. A Catalog of Quasar Properties from Sloan Digital Sky Survey Data Release 16. The Astrophysical Journal Supplement Series, 263(2):42, December 2022. doi: 10.3847 1111--383555 /acccad.
Jinyi Yang, Feige Wang, Xiaohui Fan, Aaron J. Barth, Joseph F. Hennawi, Riccardo Nanni, Fuyan Bian, Frederick B. Davies, Emanuele P. Farina, Jan-Torge Schindler, Eduardo Bañados, Roberto Decarli, Anna- Christina Eilers, Richard Green, Hengxiao Guo, Linhua Jiang, Jiang- Tao Li, Bram Venemans, Fabian Walter, Xuc-Bing Wu, and Minghao Yue. Probing Early Supermassive Black Hole Growth and Quasar Evolution with iefrrinfrared Spectroscopy of 37 Reionization-era Quasars at 6.3 < S 7.64. The Astrophysics Journal, (2:22). December 2021. doi: 01884771538-4371 ab22..
19

---

Published in Transactions on Machine Learning Research (Jan/2025)
A The evidence lower bound 2
In this section, explicitly show the derivation of expected log-likelihood term in Eq. 5
p(X) 2 f / ((F101/2)((V))0gp(XX1 (KL(())))((()()))))) log
(20)
f
p P(F|U, (U) log ( Z)dFdU = ((uu)) / plffluus Z) log
II n=
N D II (Zn,d (z(),, on dfuuu
/ g(un)) / PlSaluas 2)
Σ
log, ( fa(zn), o o
n.d
Σ flog (logp(Fn,dlfas 2n, p2
(21)
n.d
A.1 Tractability of the expectation term
The expectation term above not just factorisable but also tractable, we show this explicitly below:
f g(uu)) / p(falus, Z)
log N o2 Σ
(22)
n.d
= = g(ua)
Wogna kn Km Un, 02) In.n
2or2
n.d
En,d is scalar hd dimension of point xn), k. is matrix 1 x M the th of Knm, know that
p(fal ua, Z) = N(Knm K mm ua Knn Kmm K K=1 Kmn). Further, Sa(zn) is scalar, denoting the value at index zn of the vector fa. Inn is the th entry in the diagonal of the matrix Qun Knn Knm K=1 Kmn- Next, performing the integration with respect to quua) = =N(u,; ma Sa) yiclds,
f ((ua)) Σ log (I k Kmm us oz) 4n,n n.d
dud = Σ [[1og (og,N(nnntt ke Kmi ma, 0.)
20 In.n
20
n.d
Tr(S2An)
(23)
20
Bringing it all together, the final lower bound can be written out explicitly as,
C2 Σ
(log Alogp(an.de plan.al Rn. oz)) fallua KL(((U))lP(U)) logg)())))))))((( log
Σ logNN(zn.t R. K- mdd m., 03)
Tr(SaAn)
Σ (K(g(u)lll(aa)
20 In
20
Σ logr(zz)
(24)
where all the data dependent terms factorise enabling mini-batching of gradients. The shared model uses a sum of ELBOs L. + E, with shared latent embedding Z, which we call the joint evidence lower bound (Eq. (15p in the main paper). The additive structure ensures that the joint ELBO is also factorisable across data points N.
20

---

Published in Transactions on Machine Learning Research (Jan/2025)
Modules
Paraneters
inducing Inputs z.Z
2500
208440
nodel_spectra _swpppreeer..ereee www.raatooooa.ttaatee.eeerrrrrr hrt.iaaiinnnddssriuuutttttii wwwttribbuioo.vvriiiitttoeooeeee
model
164250
vartationll htt.aattooootttttttgt...ooomm varategy._vrrtationn wtw../rrraatonn..osrrotoooooooo. www.rttbutoo..ooooomo chol,
41062500
nodel_spectra mean hwwpe/ooddddssseeeeeee
model_
peeel_spectrreccvarrreoodde base kernel. lengthscale
10
model_ labels. wariiaaionallsssstttee.......... variiationaaaa v wis.rrbbttionnnoooooooooooooooit iisttrrbutoon..vrttttonnnl mean
1080
model_ labels. wtr.ratttooonnttttteeeeettttttto htt.:/varrrttiolloossstrrtttttoe ddsttributtoo.ccoo model
250000
Looel__aablssmman module.constame
4
model__lablls.covrr
mwdumo.d.eaaaooooo raw www_pussttttaeeeeeeeeeeee
1
model_labels.ce _ll_aalleeloa_eeeeoe bartuueebbeeeeeeee.... raw Lengthscale
10
Total Trainable Params: 41689373
Figure 7: Number of trainable parameters for the 22k model where model_ model_spectra refers to the GPs corresponding to X and model_labeil refers to the GPs corresponding to Y observation space.
Experimental set-up
In this section we detail the configuration of the experiments in Section 7 of the main paper. For each of the datasets, we repeat every experiment with five random seeds yielding different splits of the training data. The attributes of the data and sparse GP sct-up are given in Table 3 We used learning rate of 0.001 across
Dataset
D
Num inducing M Latent dim. Q
1K
996
590
120
10
22K
22844
657
250
10
Table 3: Experimental configuration to reproduce experiments in section of the main paper.
all parameters and ran the mini-batch loop with batch size of 128 for 15,000 itcrations on an Intel Core i7 processor with a GeForce RTX 3070 GPU with 8GB RAM memory. In order to give an estimate of the scale of the model for the 22k dataset we enclose a summary snapshot of the number trainable parameters our shared model.
C Testing on smaller sampled dataset
Metrics
Mean Absolute Error
Models
Baseline
Shared (ours)
Attributes
Spectra
0.0731 + 0.0020
0.0707 + 0.0033
Black hole mass
0.1882 + 0.0046
0.1765 + 0.0054
B. Luminosity
0.1731 + 0.0043
0.1658 + 0.0084
Eddington Ratio
0.1707 + 0.0024
0.1653 + 0.0026
Table 4: Summary of test-time reconstruction abilities for a smaller 1k dataset sampled with the same SNR data quality cut of SNR 10. Mean absolute error denormalised data (# standard error of mcan) evaluated on average of 5 splits with 90% of the data used for training in the 1k dataset. The performance improvement is not statistically significant for the spectra reconstruction and the bolometric luminosity prediction but favours the shared model for black hole mass and Eddington ratio.
21

---

Published in Transactions on Machine Learning Research (Jan/2025)
D Sensitivity to SNR
In Fig. we essentially plot the SNR per data point (Lest quasars) vs. the absolute error. The triangular scatter denotes a density imbalance: more data points are clustered at lower SNR values 10-15) compared to higher SNR values (>20). This just an artefact the data quality cut we deploy. The points seem evenly scattered across the range of absolute errois and SNR values. There doesn't appear to be consistent trend where absolute error systematically increases decreases with SNR. There is a very weak negative correlation (more discernable in the zoomed in plots in the bottom row) for the black hole mass and weak positive correlation for the luninosity. The unusual weak positive correlation can be betler understood in the context of the count and missing pixels analysis (Fig, 9 in each SNR bin. There are very few high SNR (> 40) quasars the dataset and further, these quasars also have more missing pixels than quasars at lower SNRs. Since the training dataset biased towards lower SNR quasars, the model generalises poorly to higher SNR examples, for luminosity and Eddington ratio. Tt highly likely that the pattern of missing pixels is similar among brighter objects and in turn among less luminous objects, making difficult for the model to adapt to this shift at test time due to the acute data imbalance.
Lummnnsity
mass
Eddington
2.4
Mean dos
0.30
0.7
02-
0.-
0.0
0.00
20
30
10
20
30
50
10
33
40
50
Luminasity
hole mazs
Eddington
Hewn 2i6 emer
Mean 122
5.0
10
20
22
30
16
10
72
16
20
22
340
Figure 8: Top: SNR error for test quasars, Bottom: SNR error zoomed in higher density regions
E Error vs. Uncertainty Calibration
In Fig. 10 quantify uncertainty calibration by computing the mean predictive uncertainty for increasing error levels. We bin test points into five evenly spaced intervals over the range of the mean absolute error for cach label. The x-axis denotes the centers of the bins for each label and the predictive uncertainty is computed as the square rool of the variance - the diagonal of the GP posterior predictive distribution. In each label, we observe the trend of predictive uncertainty increasing with absolute error, although this effect very subdued at lower error levels. Further, important to note that there is inherent noise this calculation emanating from the differing number of points in each bin; this might explain some of the variability in the bolometric luminosity uncertainty estimates.
22

---

Published in Transactions on Machine Learning Research (Jan/2025)
Bolometric Luminosity
mass
Eddington
0.5
1.0
Mean Abs. error
Moan abs. eatar
Moan abs. eitor
0.8
0.6
0.3
0.6
0.4
0.2
0.1
0.0
0.0
20
30
40
50
20
30
40
50
30
40
30
ShIt
SNI
SNIf
Average Missing Pixels per SNR Bin
100
Count per
10000
7500
5000
2500
(9.939. 15.528)
21.050]
dc (ttttt (89)))
IRc2 (95 122z (8ss (3or77 Leerr
(00s 70 1az2 (9)
(21 B1.56 26
48 uurt
54
(37
SNR Bins
Figure 9: The SNR vs. error visualised with one slandard deviation in each bin. The boltom plot shows the average number of missing pixels for the same SNR bins well as the total number of quasars in each bin
23

---

Published in Transactions on Machine Learning Research (Jan/2025)
Bolommeticc Luminooiity
Black hilte
2.24
0.13
0.35
0.3
0.100
Proitted
0072
0.44
0.2
6009
3,9755
0.157
65534
Drmed 4bokdt Dror
Figure 10: Uncertainty calibration across different error levels for each label computed over test data. The x-axis denotes centers of evenly spaced bins across the range of the MAE.
F Cross-validating M and Q
The two main parameters of our shared framework which need to be fixed at the outset are: the number of inducing points M and the latent space dimensionality Q. We set at 250 all 22k experiments after cross-validation to M 1000 and found that . 250 gave the best possible trade-off in terms of speed and accuracy. The reconstruction results with M 1000 were only marginally better than with inducing points M 250, but the compute significantly increased due the cubic scaling in inducing points.
Joint ELBOs for varying latent dim. Q
Q=2
Q=5
1000
Q=10
Q=15
600
2000
4000
6000
8000
10000
Figure 11: Sensitivity to Q: The negative ELBO objective for varying latent space dimensionality (lower is better)
In Fig. we visualize the evolution of the ELBOs across varying latent dimensionality. We notice a meaningful improvement in increasing the dimensionality from Q 2 but very marginal gains beyond this sctting in experiments. It may be important highlight that due automatic relevance
Q 10; delermination of the squared exponential kernel, setting high latent dimensionality should not degrade results as the model automatically prunes redundant dimensions by driving the corresponding inverse lengthscales to O. However, they do increase the compute cost, hence, it is important to set Q at reasonable value which is ficxible enough for structure discovery and not too constrained, while simultancously minimising the computational burden.
24


Reviews:
[{"review_id": "t1He8FeSwv", "reviewer": "Reviewer_8g1J", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "I believe the main strengths of the paper are as follows:\n1. The paper is a novel extension of Lalchand, and the changes made are well motivated by the goal of fitting a SVI-GPLVM to multimodal data. \n2. The authors experimentally evaluate their model in a wide variety of settings that appear strongly motivated by the astrophysics applications of interest. They show the ability of the model to reconstruct from a variety of corrupted samples, as well as the ability of the model to learn to predict scientific labels from quasar spectra alone. \n\nThe main weaknesses I find relate to the notation and are as follows:\n\n1. Are $\\mathbf{f}_d$ functions or scalars? $\\mathbf{f}_d(z_n)$ suggests they are functions being evaluated at $z_n$ but then I find the expressions $\\mathcal{N}(\\mathbf{f}_d, \\cdots)$ confusing as we cannot have a normal distribution over functions, it would be something along the lines of $GP(\\mathbf{f}_d, \\cdots)$.\n2. How does $ K _ {nn} $ depend on  $n$ and what is $n$? In the paper, $ K _ {nn} $ is defined as  $k _ \\theta$  evaluated on the set of points $ \\mathbf{x}_n$ for $n=1$ to $N$, however $n$ is a dummy variable in the definition so its hard to see its involvement in $ K _ {nn} $. Also, the kernel is defined on $Z$ so should these not be $ \\mathbf{z}_n$?\n3. Relatedly, I don't see definitions for $K _ {mm} $ or  $K _ {mn}$? I think from the part the play in the posterior I understand what they are but they should be defined.\n4. What is the dimension of  $K _ {mm} $? Again from reading between the lines I believe $m=M$ and so this would be an $MxM$ matrix, in which case I do not see how $\\mathbf{u}_d$ can have distribution $\\mathcal{N}(0,K _ {mm})$ unless $D=M$.\n5. In equation (3), why is the final line not $P(X|F,Z,U)$ instead of  $P(X|F,Z)$? The dependence on $Z$ here seems to go.\n6. In equation (7), is this the same kernel with different hyper-parameters instead of different kernels with different hyper-parameters? Again the notation make this aspect confusing. \n7. For the experiment showing we can predict scientific properties from spectra, how does this approach compare with a more traditional approach that just tries to predict the properties from spectra without learning a latent representation? Strong statements are made on the value of predicting such properties from spectra in the conclusion but it would seem these could be equally well supported by a standard predictive approach as opposed to a multimodal latent variable model.\n8. Finally, whilst astrophysics lies outside of my area of expertise I would say that some conclusions in the \"Scientific Interpretation and Significance\" seem very strong and if more evidence cannot be provided for these claims it would seem they should be toned down slightly. The section lays claim to what seem like two large open astrophysics problems, claiming they are resolved by this method. It would seem to me much more experimental evidence would be required for this claim. However, as stated this is not my area of expertise.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In this work the authors develop a novel version of the stochastic variational gaussian process latent variable model of Lalchand et al which is specifically for multimodal data. Whilst the original model could be applied without adjustment to multimodal data, it would require using the same kernel and hyper-parameters across modalities. If the form of the data is very different across modalities, for example visual and standard continuous covariates, this restriction could lead to poor predictive performance as the same kernel will struggle to capture both. Therefore, the authors develop an approach which is jointly trainable and allows different kernels to be used across modalities. They evaluate their approach on a selection of astrophysics applications for which the methodology was developed.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "My requested changes are:\n- The notation in section 2 needs to be cleaned up greatly as in its current state it prohibits the full understanding of the technical portion of the paper. \n- Add a comparison to a standard predictive approach in the section on predicting properties from spectra. \n- Consider toning down or further justifying claims in the conclusion.", "broader_impact_concerns": "No concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "M9nvQNX9c0", "reviewer": "Reviewer_T1nP", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\n\nThe paper considers an important application area (multi-view astronomical data). The model seems well suited for this data and can lead to interesting scientific results in astronomy, as demonstrated by the extensive experiments in the paper. \n\n## Weaknesses\n\n- To my understanding, the model is a very simple extension of Lalchand et al. 2022 (see \"summary of contributions\"). In particular, the model is the same as Lalchand et al. 2022, just with different hyperparameters for the different data dimensions. As such, the paper should be much clearer about this. As a result, the paper may be better suited to an astronomy focused journal.\n- Some of the exposition of the model is confusing. For example, in eq. 14, the GPs $f_d$ are integrated out. But then in eq. 16, the actual objective, the expectation is taken with respect to the variational distributions. Further, eq. 16 is not evaluated wrt to this expectation. To understand the ultimate objective function, I had to go back to Lalchand et al. 2022 (e.g. their eq. 6). \n- Lalchand et al. 2022 consider point, MAP and Bayesian options for estimating $z$. In many cases, they find the Bayesian version (i.e. specifying a variational distribution $q(z)$ and learning the variational parameters) has better reconstruction. Why did the authors not consider this version?\n- p. 2: \"We don't employ this setting as it is not straightforwardly applicable to missing data contexts.\" Masked autoencoders are often used to learn representations - can they not be applicable here? Moreover, in multimodal settings, autoencoders are now regularly used to predict one modality from another modality. Can the authors speak more to why they only consider GPLVMs, and not another VAE-based multi-view model?\n- Is the per-column conditional independence assumption on $x$ appropriate for pixel data? \n- Some of the details in the experiments are difficult to understand as a non-astronomer (see requested changes).", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes an extension of the stochastic variational Gaussian process latent variable model (GPLVM) (Lalchand et al. 2022) to the multi-view setting.\n\nThe application is in astrophysics - specifically, the authors consider very luminous galaxies, for which they have both spectral pixel data and tabular data measurements such as black hole mass, accretion rate, and luminosity. \n\n### Model\n\nIn GPLVM, the coordinates of the data $x = (x_1,\\dots, x_D)$ are modeled independently. That is, \n$$p(x|f, z) = \\prod_{d=1}^D N(x_d; f_d(z), \\sigma^2)$$\nwhere $f_d$ are GPs and $z$ is the latent representation of $x$. Because the coordinates are modeled independently, the GPVLM model naturally can incorporate another data view $y=(y_1,\\dots, y_L)$, simply by concatenating the data $(x_1,\\dots, x_D, y_1,\\dots, y_L)$ and adding additional GPs $f_l$, $l=1,\\dots, L$. This is what this paper proposes (where the $f_l$ are allowed different kernel hyperparameters than $f_d$).\n\nFor optimization, the paper follows Lalchand et al. 2022. Specifically, the paper uses the \"point\" version of the Lalchand et al. 2022 algorithm, which obtains point estimates for the latents $z$. \n\nThe authors conduct extensive experiments demonstrating:\n- on test data, the model reconstructs quasar properties similarly or better than the Lalchand et al. 2022 model\n- on test data, the model reconstructs quasar spectra well even when the data is partially unobserved\n- on test data, the model predicts quasar properties well from only the quasar spectra\n- the model can generate example spectra for given scientific properties\n\n### References\n- Vidhi Lalchand, Aditya Ravuri, and Neil D. Lawrence. (2022) Generalised GPLVM with Stochastic Variational Inference. Artificial Intelligence and Statistics.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "No", "audience_explanation": "", "requested_changes": "### Critical changes\n- To my understanding, the model is a very simple extension of Lalchand et al. 2022 (see \"summary of contributions\"). In particular, the model is the same as Lalchand et al. 2022, just with different hyperparameters for the different data dimensions. As such, the paper should be much clearer about this.  This would also improve exposition of the model.\n- eq. 15 and eq. 16: this is confusing - in one, you integrate over $f$, in the other, you don't. (see also comments in weaknesses)\n- Can the authors consider the MAP and Bayes version of stochastic variational GPLVM? \n- Fig 1: I'm confused by the self referential arrow on the $f_d$ - aren't the coordinates conditionally independent given the data?\n- Experiments: some paragraphs are a little confusing for readers without astronomy knowledge. For example: \"Note that we omit the reconstruction of the redshift label z, as we do not expect any redshift evolution in the spectral shape of quasars, but rather the GPLVM is likely learning where spectral pixels are measured vs. missing. Since we are using data from SDSS, all spectra have the same observed wavelength coverage, which is then shifted to the rest-frame wavelengths by dividing by (1 + z), and thus the observed pixel dimensions in our data set contains the redshift information, rather than any spectral differences.\" As a non-astronomer, I do not understand this. Why do you not expect redshift evolution in spectra? What are the rest-frame wavelengths and why are they shifted? You then say that pixels contain redshift information, but earlier you said that you don't expect redshift evolution in the spectra pixels?\n- Table 2: can you compare to Lalchand et al. 2022 too here?\n- p. 14: \"Our new model allows us to determine the masses of SMBHs for quasars at all redshifts from single-epoch data and does not require coverage of specific emission lines, since it can handle missing data.\" How does the model handle all redshifts? What if there are redshifts missing in the training dataset?\n- p. 2: \"We don't employ this setting as it is not straightforwardly applicable to missing data contexts.\" Masked autoencoders are often used to learn representations - can they not be applicable here? Moreover, in multimodal settings, autoencoders are now regularly used to predict one modality from another modality. Can the authors speak more to why they only consider GPLVMs, and not another VAE-based multi-view model?\n\n\n### Other changes\n- I'm confused by the inducing point distribution $u_d|\\widetilde{Z} \\sim N(u_d; 0, K_{mm})$. Is the kernel computed with respect to $\\widetilde{Z}$? This should be made clear.\n- eq. 13: what happened to conditioning on $Z$?\n- Experiments: can you say \"test\" $X^*, Y^*$ instead of \"unseen\"? \n\n### Typos\n- p. 2: \">400 x bigger\" re-write: datasets over 400 times larger\n- p. 3: \"Essentially by keeping the representation of $u_d$ uncollapsed.\" This sentence isn't complete and is also confusing - in what sense is $u_d$ uncollapsed and why does that help?\n- p. 4: equation (3) seems superfluous as it was derived in Lalchand et al. 2022\n- p. 4: \"yiels\" -> \"yields\"\n- p. 4: \"we introduce the idea of an additional observation space...\" this could be made more clear e.g. \"we extend this model to accommodate an additional data view\"\n- p. 6 \"variaitonal\" -> \"variational\"\n- change your `\\citep` to `\\citet` when for in-text references\n- p. 14: \"few hundred Myr old\" - million years?", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "UKkdOKZweU", "reviewer": "Reviewer_qTSY", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\n- The main strength of this work is in its application, which seems to be a novel and impactful place to apply recent GPLVM work. In particular, the improvement of using SVI over existing work in the astrophysics literature such as (Eilers et al., 2022) is apparent w.r.t. scalability.\n- The technical contributions are clearly stated and developed.\n- The provided code is well-written and self-documented, which I view as a big plus. To this end, the pseudo-code in Algorithm 1 is also very clear.\n## Weaknesses\n - (General) Compared to (Lalchand et al., 2022), there isn't much technical novelty, mainly showing that sharing inducing points still leads to a factorized ELBO. This result is unsurprising, since the GPs between dimensions are still conditionally independent, just with different hyperparameters, and indeed, this is how the derivation proceeds.\n   Nevertheless, this is a well-motivated approach to a novel application, and technical novelty isn't a criterion for acceptance, so this is not a big issue.\n - (General) Many of the claims made about the scientific significance of this work are not easily verifiable by someone in the machine learning community. In particular, some of the comments in Section 8 read as claims rather than speculation; they may be straightforward to an astronomer, but would not be known to the TMLR audience. \n   \n    For example, the final claim that \"[the] new model allows us to now also use quasars as standard candles, [...] enabling new constraints on the dark energy content and Hubble constant [in the future]\" seems like a strong claim, but has no further elaboration or reference. If this claim relies on the accurate reconstruction of the absolute luminosity, it seems very important to experimentally verify this. \n   \n   Similarly, it is claimed that \"[the new method] allows us to determine the masses of SMBHs for quasars at all redshifts from single-epoch data and does not require coverage of specific emission lines, since it can handle missing data. This opens up [new possibilities...].\" My non-expert understanding (based mostly on Section 7.1.1 of (Eilers et al., 2022) and the references therein) is that single-epoch estimators of the SMBH mass do exist, but are unreliable or inapplicable with high redshifts. But in the current manuscript, it doesn't seem like any results are evaluating the prediction quality of the SMBH mass w.r.t. redshifts, which seems central to the claim. Moreover, if the claim is that the proposed method reconstructs the black hole mass *more* reliably than existing methods, this should be verified.\n- (Section 7.2) It is hard to judge the quality of uncertainty quantification here, especially with respect to black hole mass. It would be nice to see a calibration plot.\n- (Section 7.2) Likewise, the claim that the error is effectively independent of SNR is not easy to verify from Figure 5. This claim is counterintuitive to me, and should be explicitly verified, e.g., by plotting the average error vs. SNR. When making this claim, it also seems important to note that a cut of SNR > 10 was applied to the data.\n- (Appendix B) It's great that this work acknowledges some of its limitations in Appendix B, as well as some future directions to address them. I think it would strengthen the paper to include this brief discussion somewhere in the main text.\n## Other Questions\n- (General) The proposed method is not constrained to just (X, Y) -- is there a reason to consider all scientific labels to have the same kernel hyperparameters?\n- (Code) It seems like `spectra_dims` in `baseline_gplvm` (which is only used in the `partition` option, and doesn't seem to appear in the paper) refers to parts of the latent space are given a \"rougher\" kernel, hoping this corresponds to the label outputs. Is this a correct interpretation? If so, the variable name `spectra_dims` may be a bit confusing/misleading, and it might help clarity to just get rid of this unused code.\n- (Equation 12) Why is there a $\\cdots$ in the integrals?\n- (Section 2.1) Should the paragraph after Eq. (1) be referring to latent points $z, z'$ instead of $x, x'$?\n- (Section 4) Why are there circular arrows for $f_d$ and $f_l$? And should both be in their own plates (over $d=1,\\dots, D$ and $l = D+1, \\dots, D+L$, respectively)?\n- (Figure 5) The caption figure says the error bars denote one standard deviation, while the text says $2\\sigma$. Looking at `visualization.py`, it seems like it may actually be $1.96\\sigma$. Please clarify what is plotted, as this is crucial to interpreting the figure.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces a scalable generative model for quasars based on GPLVMs with SVI. The core technical observation is that \"scientific labels\" such as the black hole mass or Eddington ratio of a quasar and the quasar's spectra may be modeled with separate GP decoders using the same latent embedding as spectra, i.e., have a \"shared\" latent representation. The proposed method is shown to hold in general for observation spaces (X, Y), and then applied specifically to quasars, where it is shown that (1) the spectra are effectively captured in latent space, (2) the spectra may be reasonably interpolated, even with large amounts of missing data, and (3) the scientific labels may be estimated using only the spectra. These claims are evaluated on a large dataset consisting of ~22,000 astronomical measurements, and their scientific implications are discussed.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "## Critical\n- (General) Regarding claims of scientific interpretation and significance, please temper claims or provide evidence, as discussed in the weaknesses section above.\n- (Section 7.2) Please provide more direct evidence for the claim that the reconstruction accuracy of scientific labels is robust and independent of the SNR. In particular, please also comment on the relevance of the SNR > 10 cut to the data.\n- (Figure 5) Please clarify what the uncertainty bands represent (i.e., $1\\sigma/1.96\\sigma/2\\sigma$).\n\n## Other Suggestions\n- (General) All figures are a bit blurry. I think it would enhance the presentation to save figures in vector format or with a higher DPI.\n- (General) Please consider incorporating the limitations discussed in Appendix B into the main text.\n- (Section 1) The last sentence of the dimensionality reduction paragraph (i.e., \"These methods however, are not generative [...]\") was a bit difficult to parse -- I'd suggest rephrasing it for clarity.\n- (Section 1) I'm not convinced of the claim that \"generative latent variable models of late have supplanted traditional dimensionality reduction techniques\"; e.g. the UMAP paper is recent and very highly cited, but is not a generative latent variable model. I would suggest using a word other than \"supplanted\" here, as this claim is not otherwise important to the paper.\n- (Section 3) I think I understand that missing dimensions are algorithmically handled by simply not computing the corresponding gradients -- i.e. if the dimension $d^*$ is missing data point $x_*$, then the ELBO is updated as $\\mathcal{L}_X \\gets \\mathcal{L}_X + \\log p(z^{*}) + \\sum_{d \\neq d^{*}} \\mathbb{E}_q \\big[ p(x^{*}_{d} \\,|\\, f_d, z^{*}, \\sigma_x^2) \\big]$. This seems pretty clear by the use of `GaussianLikelihoodWithMissingObs` in the code, but explicitly highlighting this in Section 3 might improve clarity.\n- (Algorithm 1) Semantically, the Training Framework and Prediction Framework parts of Algorithm 1 should probably be split into two separate algorithms, since they each have separate inputs and outputs.\n- (Figure 6) Please update the caption of Figure 6 so that the content of the plot is self-contained, particularly since its description in Section 7.3 succeeds the plot itself.\n\n## Minor Typos\nI also noticed a few possible typos in the text, which I hope are useful to the authors.\n- (Section 3) The second arrow in $x^* \\to z^* \\to p(y_l^* \\, | \\, z^*)$ is longer than the first.\n- (Algorithm 1) In Step 2 of Prediction, \"$\\gets$\" and \"$=$\" are used inconsistently.\n- (Section 7.3) I suggest that \"black hole masses in the range $\\log_{10}(M_\\cdot /  M_\\odot) = 7.9-10.0$\" should read \"black hole masses $\\log_{10}(M_\\cdot /  M_\\odot)$ in the range $[7.9, 10.0]$\".\n- (Section 2.1) Before Equation (4), there is a typo in \"this yields\".\n# References\nEilers, A. C., Hogg, D. W., Schölkopf, B., Foreman-Mackey, D., Davies, F. B., & Schindler, J. T. (2022). A Generative Model for Quasar Spectra. _The Astrophysical Journal_, _938_(1), 17.\n\nLalchand, V., Ravuri, A., & Lawrence, N. D. (2022, May). Generalised GPLVM With Stochastic Variational Inference. In _International Conference on Artificial Intelligence and Statistics_ (pp. 7841-7864). PMLR.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: All three reviewers agree that the final version deserves to be accepted. I share their opinion, especially with the additional experimental evaluations and improvements in the text. As all these have already been included in the current version a minor revision is not necessary and I recommend acceptance.
Decision: Accept (poster)