# Learning Audio Concepts from Counterfactual Natural Language

###### Abstract

Conventional audio classification relied on predefined classes, lacking the
ability to learn from free-form text. Recent methods unlock learning joint
audio-text embeddings from raw audio-text pairs describing audio in natural
language. Despite recent advancements, there is little exploration of
systematic methods to train models for recognizing sound events and sources in
alternative scenarios, such as distinguishing fireworks from gunshots at
outdoor events in similar situations. This study introduces causal reasoning
and counterfactual analysis in the audio domain. We use counterfactual
instances and include them in our model across different aspects. Our model
considers acoustic characteristics and sound source information from human-
annotated reference texts. To validate the effectiveness of our model, we
conducted pre-training utilizing multiple audio captioning datasets. We then
evaluate with several common downstream tasks, demonstrating the merits of the
proposed method as one of the first works leveraging counterfactual
information in audio domain. Specifically, the top-1 accuracy in open-ended
language-based audio retrieval task increased by more than 43%.

Index Terms—  sound event detection, audio understanding, multimodal
representations, free-form text, counterfactual representation learning, audio
captioning

##  1 Introduction

Conventional audio processing in machine learning relies on predefined
categories, limiting the ability of the models to understand audio nuances
using descriptive text. These constraints limit open-ended and contrastive
training for better audio-text alignment. New trends in the field improve
classic models using audio-text learning from audio data and matching natural
language descriptions that have become successful models in image-text tasks
[1]. Learning audio representations from pairs of audio and their textual
descriptions facilitates the development of foundational models for audio
tasks, enabling audio-text models to generalize beyond the confines of
predefined classes, leveraging natural language descriptions [2].

Fig. 1: Counterfactual reasoning helps to distinguish various sound sources in
an audio signal captured by a microphone. We identify isolated sound sources
using GPT-3.5-Turbo. Subsequently, we intervene to alter one or more of the
sources of sounds to construct an imaginative linguistic representation in an
alternative world that could have happened if there were other objects
instead, thereby eliminating dependence on empirical audio data to reason
objects driving the acoustic waves.

Advancements in audio-text representation include AudioCLIP [3], a tri-modal
model, and Wav2CLIP [4], which extends CLIP [1] to audio. Subsequently,
Elizalde et al. proposed CLAP (Contrastive Language-Audio Pretraining) [5], a
method that takes its inspiration from successful image-text models [1]. CLAP
is able to train on audio and text directly without relying on image data for
the learning process. However, human annotated audio captioning data is
expensive and time-consuming to acquire. Recently, with large language models
(LLMs) such as ChatGPT, an upgrade version of GPT-3 [6] fine-tuned to follow
human instructions [7], have become popular and been utilized to augment
learning for various domains, both LAION-Audio-630K [8] and WavCaps [9]
leverage the powerful textual re-writing and editing capabilities in order to
acquire more audio-text pairs. Furthermore, methods such as Pengi [2] and
listen, think, and understand [10] have built upon these foundational
methodologies.

(a) CLAP [5]

(b) Our method

Fig. 2:  (a) The CLAP method [5] utilizes audio captions to train audio
embeddings. (b) Our method leverages a prompt, represented by
$p=\\{p_{1},p_{2}\\}$, using the GPT-3.5-Turbo model as LLM, that elicits
counterfactual captions. This model guides interventions on the existing
captions. The overarching goal is to pinpoint the sources of acoustic waves
through human-narrated captions, shown as identification. After identifying
these sources, controlling modification and intervention are performed by
$p_{2}$ to incorporate the resulting counterfactual scenarios into a causal
learning framework. This technique significantly improves the capability of
audio-text models to distinguish subtle variations in text to align with
sounds emanating from different objects.

Distinguishing between sounds in similar conditions, such as the sound of a
firework and a gunshot in the same concert event, requires controlled trials
of both sounds and a learning algorithm to differentiate between them in the
same context. However, existing audio-text datasets lack alternative scenarios
for ethical and practical reasons. Counterfactual reasoning has been utilized
to improve multimodal models involving vision modality [11] and to aid
grounding concepts within visual objects [12]. The semantic differences
between pairs of similar but slightly different audio clips have been
discussed in [13] and [14], addressing audio captioning. To the best of our
knowledge, , our work is pioneering in utilizing the knowledge base and
capabilities of LLMs to integrate counterfactual reasoning. This provides data
and methods to enhance the learning of audio-text correspondence with
counterfactual information. The proposed method utilizes counterfactual
language and multimodal embeddings to improve textual discriminability in
audio processing tasks, as shown in Figure 1. Subsequently, we develop a
composite loss function that integrates the concept of triplet angular spaces
and enforces factual audio-text consistency. We are inspired by recent
developments in causality using LLMs [15] that offer new possibilities, such
as augmentation of meaningful counterfactuals in situations where audio
learning data may not be available. The potential of the proposed method is to
influence a wide array of applications, including automatic speech
recognition, sound event detection, and audio-visual scene understanding. The
core innovation of our work lies in being the first to integrate causal and
counterfactual models into the analysis of audio data. Generated
counterfactual sentences and prompts attempted will be available at
https://github.com/ali-vosoughi/counterfactual-audio.

##  2 Learning Audio from Counterfactual

Causality describes the relationship between a cause and its subsequent
effect. This study focuses on the causal relationships within audio samples
that can be inferred from their corresponding text captions. Specifically, we
address scenarios where data acquisition of counterfactual audio is impossible
or costly by leveraging natural language as a substitute for imaginative data.
Our method aligns with recent advancements in causality in natural language
that aim to extend domains of causality to LLMs [16, 15]. We first explore
prior work of CLAP [17] that motivated us and then propose our method as the
first work to extend counterfactuals to audio learning.

Contrastive Language-Audio Pretraining (CLAP): Given a batch of $N$ pairs
$(x_{i},y_{i})$, $x_{i}$ represents an audio sample and $y_{i}$ is its
corresponding text caption for $i=1,\ldots,N$. We define the audio and text
encoders as $E_{a}(i)=\phi_{\text{audio}}(x_{i})$ and
$E_{t}(i)=\phi_{\text{text}}(y_{i})$, respectively. These encoders transform
each $x_{i}$ and $y_{i}$ into vectors in $\mathbb{R}^{d}$, which are
aggregated to form matrices $E_{a}\in\mathbb{R}^{N\times d}$ and
$E_{t}\in\mathbb{R}^{N\times d}$. The similarity matrix $C$ in the CLAP
framework is given by (1) [5]:

$C=\tau\left(E_{t}\cdot E_{a}^{\top}\right)$ (1)

The CLAP loss minimizes the discrepancy between audio and text
representations, is denoted by $\mathcal{L}$ and formulated as:

$\mathcal{L}=0.5(\ell_{\text{text}}(C)+\ell_{\text{audio}}(C))$ (2)

Here, $\tau$ serves as a scaling factor in (1), modulating the effect of the
similarity scores. While the original CLAP framework, shown in Figure 2(a), is
valuable for audio-text pretraining, it does not allow causality expression.
Addressing this gap, we introduce a counterfactual natural language strategy
to infuse the CLAP framework with causality. This addition is beneficial as it
circumvents the difficulties associated with data collection for
counterfactual audio analysis.

Causal Identification and Intervention: The counterfactual sentences in our
model are generated through a prompt-based intervention on an observed caption
$y$, represented as $y^{\star}=\text{do}(y|p)$. This prompt $p$ is designed to
fulfill three key aspects: it is factually grounded, identifies acoustic
sources in captions to serve as identifying causes, and manipulates these
sources to alter the caption, to serve as causal interventions.
Identifiability means that all causal sources of acoustic waves can be
obtained only from the language [15]. Examples of these transformations are
depicted in Table 1.

Table 1: Samples of original captions from the Clotho [18], MACS [19], and AudioCaps [20] datasets and their counterfactual pairs. Dataset | Original Caption
---|---
Generated Counterfactual
Clotho | A gun is loaded, then loaded by hand some more
A piano is played, then played by hand some more.
A few gunshots are fired at the target shooting range
A few fireworks light up the night sky at shooting range.
AudioCaps | An adult male speaks and a crash occurs
An adult male speaks and a thunderstorm rumbles.
Large group of people clapping
Flock of birds chirping in unison.
Idling car, train blows horn and passes
Dogs barking, train blows horn and passes.
MACS | A crowd of people indoors talking
A group of cars honking on a busy street.
Adults and children are walking and talking
Cars and trucks are honking and zooming.
Adults talking and some footsteps coming across
Dogs barking and some footsteps coming across.

Control Mechanisms for Counterfactual Language: We utilize prompts based on
the Chain-of-Thought (CoT) method [21] to align with objectives of causal
identification and intervention and generating counterfactuals. Figure 2(b)
shows our prompt design’s rationale, objectives, and context. We introduce a
two-step prompting mechanism, denoted as $p=\\{p_{1},p_{2}\\}$, as shown in
Fig. 2(b). In this mechanism, $p_{1}$ anchors the discussion in factual
elements, dissects acoustic objects, and plays the role of causal identifier.
Concurrently, $p_{2}$ governs the generation of counterfactual statements by
intervening in the identified causal acoustic sources. The decomposition of
the audio captions to their acoustic sources and objects based on the physics
of the acoustic waves using $p_{1}$ will ground the generated output of the
LLM to physically possible scenarios. The counterfactuals can range from full
negative examples, as found in hard negative sampling [22], to minor,
physically plausible counterfactual scenarios. The steering of counterfactuals
will influence the learning of the audio-text embeddings. For instance, in the
example of "children and adult voices with footsteps and birds singing in the
background," $p_{2}$ can be used to change the primary sound, inferred from
the caption, "children and adult voices," background sound, "with footsteps,"
or ambient sound, "birds singing in the background." In this paper, we use a
combination of these assumptions without restricting the assumption for the
sake of generalizability.

Loss Functions to Incorporate Counterfactuals: The angle loss aims to minimize
the angular difference between factual and counterfactual captions. It is
defined as:

$\displaystyle
L_{\text{angle}}=\frac{1}{N}\sum_{i=1}^{N}\max\Big{(}0,\cos(\phi_{\text{audio}}(x_{i}),\phi_{\text{text}}(y^{\star}_{i}))$
(3)
$\displaystyle-\cos(\phi_{\text{audio}}(x_{i}),\phi_{\text{text}}(y_{i}))+\mu\Big{)}$
(4)

Here, $\mu$ represents the angular margin. To encourage factual consistency
between audio samples and their corresponding captions, we define the factual
consistency loss as follows:

$\displaystyle
L_{\text{factual\\_consistency}}=\frac{1}{N}\sum_{i=1}^{N}\lVert\phi_{\text{audio}}(x_{i})-\phi_{\text{text}}(y_{i})\rVert^{2}_{2}$
(5)

The total loss, combining both the factual consistency loss and angle loss, is
expressed in (6):

$L_{\text{total}}=w_{1}\cdot L_{\text{angle}}+w_{2}\cdot
L_{\text{factual\\_consistency}}$ (6)

Choices of hyperparameters in (6) are pragmatic, serving as a trade-off to
best exploit the capabilities of counterfactuals while ensuring factual
consistency.

##  3 Experimental design

###  3.1 Encoders

Audio encoder: We use PANNs [23] encoder, specifically ResNet-38 has been used
with pretrained weights loaded with adapter layers to fine-tune model and
align the embeddings.

Text encoders: We use the same CLIP text encoder modules [1] provided from
HuggingFace111https://huggingface.co/docs/transformers/model_doc/clip [24] for
encoding captions and counterfactuals. The weights of the encoders were frozen
in all stages. Therefore, we excluded the effect that may arise form the
encoding performance of the language encoders.

We employed logarithmic Mel spectrograms of audio sampled at 32kHz. The hop
size is set to 320 frames, the window size to 1024 frames, and we utilized 64
Mel bins spanning the frequency range of 50-14000 Hz. Audio clips were
randomly truncated to contiguous 10-second segments for training purposes,
with zero padding applied for shorter clips. The captions remained unaltered.
During training, batches containing pairs of audio and text are randomly
selected.

###  3.2 Data

Training datasets: Total of 44,292 from AudioCaps [20], 29,646 pairs from
Clotho [18] (each audio has five captions, so we created five pairs per clip),
and 17,276 pairs from MACS [19] have been used during pretraining. The reason
that we picked these three datasets are that these are purely annotated by
humans.

Test datasets: We use the Clotho dataset [18] for evaluating the model’s
performance in language-based audio retrieval task. We generated
counterfactual captions using the GPT-3.5-Turbo, similar to the training sets.
For evaluating the model’s performance as zero-shot classification in
conventional problems with limited classes, we use the Environmental Sound
Classification 50 (ESC-50) [25] with 50 predefined categories related to audio
classes, and UrbanSound8K (US8K) [26] with ten classes as air conditioner, car
horn, children playing, dog bark, drilling, engine idling, gunshot,
jackhammer, siren, street music.

###  3.3 Baseline

We adopt the approach from CLAP [5], and train our version with the same
datasets we used in generating counterfactuals, including only AudioCaps,
Clotho, and MACS.

##  4 Results and Discussions

###  4.1 Evaluation on Downstream Tasks

Results on Clotho: We use Clotho [18] to test the performance of our method on
language-based audio retrieval task. As listed in Table 2, our method’s
performance on text-based retrieval tasks yields a 43% improvement in top-1
accuracy, reinforcing its superior precision. Notably, the performance has
slight improvement for top-10 retrieval tasks. We used the cosine similarity
of text and audio embeddings to measure the performance. Therefore, the text
encoder plays a crucial role in capturing small spelling nuances that have
significant effects, and the limitations of the text encoder in capturing
nuances will impact the model’s performance. For instance, the cosine
similarity of embeddings of "This is the sound of a dog" vs. "This is the
sound of a cat" using BERTscore text encoder [27] is 0.99, while these two
sentences are carrying two significantly distinct animals. While this study
excludes the effects of different text encoders to focus on the merits of
counterfactual audio analysis, a comprehensive token dictionary and text
encoder can improve sentence distinction.

Table 2: Performance (out of 1) on the Clotho evaluation set for text-to-audio retrieval is listed. Method | Top-1 | Top-10
---|---|---
CLAP | 0.088 | 0.395
Our method | 0.126 | 0.423

Results on ESC-50 and US8K: We evaluated the zero-shot classification
performance of our proposed model on two benchmark datasets, ESC-50 and US8K,
and summarized in Table 3. As the table shows, our model performs commendably
on the ESC-50 dataset, which features many classes. This performance is
slightly better than that of the CLAP method. Conversely, the performance lags
on the US8K dataset when compared to CLAP. One reason might be that the number
of classes in US8K is much lower, in contrast to what our model learned during
the training. As the class labels of the US8K lack sufficient textual detail
about the data, the classification becomes harder as all labels get relatively
similar and high similarity scores. For instance, the BERTscore distance [27]
between class label "siren" and all other class labels is 0.844 $\pm$ 0.020,
and the trend consistently repeats for all class labels of the dataset.
Contrastively, for the ESC-50 dataset, the value of similarity of the class
label "siren" to all other classes is 0.821 $\pm$ 0.018. Further analysis is
necessary to ascertain the reasons behind the reduced performance on the US8K
dataset and to refine our model for enhanced overall accuracy.

Table 3: Comparison of zero-shot performance between our proposed method and existing methods on ESC-50 [25] and US8K [26]. Method | ESC-50 | US8K
---|---|---
Wav2CLIP [28] | 0.414 | 0.404
AudioClip [29] | 0.694 | 0.653
CLAP | 0.729 | 0.798
Our method | 0.744 | 0.475

###  4.2 Ablation Studies

Figure 3 shows the t-SNE visualizations of the embeddings for various settings
of the coefficients. Evidently, starting from a random guess, the embeddings
evolve with each subsequent addition of the loss terms. The plots in Fig. 3
show that when all loss functions are set to zero, the audio encoders are
based on the PANN embeddings in ours, and text encoders which are frozen and
are showing the fixed space of captions and counterfactuals. One particular
observation is that audio embeddings of the CLAP are closer to the fact
(orange dots) than counterfactuals (blue dots). This observation, by itself,
shows that CLAP is favorably learning to stay closer to facts. One reason is
that our generated counterfactuals combine various types of counterfactuals,
and some of the counterfactuals are similar to negative, so CLAP can
successfully work on them. By incorporating counterfactuals via angle loss, we
observe that the audio embeddings of our method get distant from the
counterfactuals, staying closer to the facts. Alternatively, by only including
factual loss, audio embeddings tend to align fully with factual embeddings,
staying closer to the facts. Noteworthy is the trade-off that adding the
factual consistency loss reduces the distance to the counterfactuals; however,
the distance to the facts becomes much less, and we observe this as a trade-
off. Having both factual consistency loss and angle loss ensures that the
counterfactuals are distant enough while facts are kept closer to the audio
embeddings.

$w_{1}$ | $w_{2}$ | Top-1 | Top-10 | $sim(x,y)>sim(x,y^{\star})$
---|---|---|---|---
0 | 0 | 0.0002 | 0.01 | $517.80\pm 237.36$
1 | 0 | 0.0819 | 0.3365 | $1000.80\pm 5.76$
0 | 100 | 0.1328 | 0.4379 | $861.80\pm 293.76$
1 | 100 | 0.1102 | 0.3782 | $967.80\pm 43.36$

Table 4: Ablations on different combinations of $w_{1}$ and $w_{2}$. The last
column represents the effect of each loss on the variation in audio-text
similarity versus audio-counterfactual by changing the loss terms.

In order to quantify the closeness of audio embeddings to original text and
counterfactual text, we use the cosine similarity of audio-text and audio-
counterfactual embeddings. Based on these similarities, we count the number of
times over 1044 samples of the Clotho evaluation set that the audio embeddings
are closer to facts than counterfactuals. As expected in the first row, the
random guess is 517, roughly about half of the 1044 samples in the test set.
In Table 4, adding angle loss $w_{1}$ and using counterfactuals improve our
audio embeddings. Another pattern is that having factual consistency loss
improves accuracy; however, the closeness of audio-text, as compared to
counterfactuals, is not as good as when we add angle loss. Therefore, by
observing patterns in Fig. 4 and Table 4, there is a trade-off between facts
and counterfactuals. By getting closer to facts through factual consistency
loss, the accuracy of Clotho retrieval increases, while the model does not
learn to distinguish well between the facts and counterfactuals. In contrast,
adding angle loss incorporates counterfactuals in training, helping to
distinguish between descriptions that might be directly (as facts/captions) or
indirectly (as counterfactual captions) related to the same audio. This trade-
off is conformal with intuition. Imagine a judge evaluating two cases. One
case presents strong evidence (facts), and the other presents a scenario that
might have happened but did not (counterfactual). If a judge focuses solely on
ensuring the judgments align perfectly with the facts, then might struggle to
distinguish between the two cases effectively. However, considering both the
facts and the counterfactual scenario, the judge becomes better at discerning
their nuances. This balance between staying close to facts and acknowledging
counterfactuals is what our model learns during training.

(a) $w_{1}$: 0, $w_{2}$: 0

(b) $w_{1}$: 1, $w_{2}$: 0

(c) $w_{1}$: 0, $w_{2}$: 100

(d) $w_{1}$: 1, $w_{2}$: 100

Fig. 3: t-SNE visualization of audio embeddings with ours and CLAP compared to
original and counterfactual caption embeddings under different parameter
configurations. Visualization keys: factual text (red dots), counterfactuals
(blue), our audio embeddings (green), and CLAP audio embeddings (orange). As
loss terms are incrementally introduced, our audio embeddings consistently
align more closely with factual data and distance from counterfactuals. Our
audio embeddings get closer to facts and distance from the counterfactuals for
various combinations of angle loss and factual consistency loss terms.

It is important to note that while some existing methods may outperform ours,
e.g., on the US8K dataset due to the use of larger datasets. We argue that
these methods often overlook the importance of causal reasoning, tending to
improve accuracy rather than introducing counterfactuals as an insightful and
novel method in bridging counterfactuals to the audio community. In contrast,
our method considers both the data and its underlying causal factors,
providing a more robust and insightful audio-text representation by borrowing
identification and intervention from the science of causality.

##  5 Conclusion and Future Direction

For the first time, we incorporate a counterfactual framework into the audio
domain. We leverage LLMs for counterfactual reasoning by prompt-based
intervention on the identified acoustic objects. This integration aims to
identify variations in audio-text representations by focusing on natural
language as a surrogate for the lack of alternative audio-text data when
acoustic waves’ origin and root cause vary. Our method exploits human-
generated reference captions for surrogate counterfactuals and adopts them to
the audio-text pretraining in a triplet model with factual consistency.
Counterfactual natural language effectively compensates for the scarcity of
comprehensive counterfactual audio data for ethical or feasibility reasons and
enhances the distinguishability of audio-text models. Empirical evaluations
using the various datasets substantiated the effectiveness of our method. In
particular, our method yields a 43% improvement in top-1 accuracy for open-
text tasks. Future research may explore the efficacy of these counterfactuals
in challenging existing factual representations and their subsequent impact on
audio-text correlation. Another avenue for future work could involve examining
various levels of counterfactual reasoning.

## References

  * [1] Alec Radford, Jong Wook Kim, Chris Hallacy, Aditya Ramesh, Gabriel Goh, Sandhini Agarwal, Girish Sastry, Amanda Askell, Pamela Mishkin, Jack Clark, et al.,  “Learning transferable visual models from natural language supervision,”  in International Conference on Machine Learning. PMLR, 2021, pp. 8748–8763.
  * [2] Soham Deshmukh, Benjamin Elizalde, Rita Singh, and Huaming Wang,  “Pengi: An audio language model for audio tasks,”  arXiv preprint arXiv:2305.11834, 2023.
  * [3] Andrey Guzhov, Federico Raue, Jörn Hees, and Andreas Dengel,  “Audioclip: Extending clip to image, text and audio,”  in ICASSP 2022 - 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2022, pp. 976–980.
  * [4] Ho-Hsiang Wu, Prem Seetharaman, Kundan Kumar, and Juan Pablo Bello,  “Wav2clip: Learning robust audio representations from clip,”  in ICASSP 2022 - 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2022, pp. 4563–4567.
  * [5] Benjamin Elizalde, Soham Deshmukh, Mahmoud Al Ismail, and Huaming Wang,  “Clap learning audio concepts from natural language supervision,”  in ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2023, pp. 1–5.
  * [6] Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al.,  “Language models are few-shot learners,”  Advances in neural information processing systems, vol. 33, pp. 1877–1901, 2020.
  * [7] L. Ouyang _et al._ ,  “Training language models to follow instructions with human feedback,”  Advances in Neural Information Processing Systems, vol. 35, pp. 27730–27744, 2022.
  * [8] Y. Wu, K. Chen, T. Zhang, Y. Hui, T. Berg-Kirkpatrick, and S. Dubnov,  “Large-scale contrastive language-audio pretraining with feature fusion and keyword-to-caption augmentation,”  in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2023, pp. 1–5.
  * [9] X. Mei _et al._ ,  “Wavcaps: A chatgpt-assisted weakly-labelled audio captioning dataset for audio-language multimodal research,”  arXiv preprint arXiv:2303.17395, 2023.
  * [10] Yuan Gong, Hongyin Luo, Alexander H Liu, Leonid Karlinsky, and James Glass,  “Listen, think, and understand,”  arXiv preprint arXiv:2305.10790, 2023.
  * [11] Mehdi Zemni, Mickaël Chen, Éloi Zablocki, Hédi Ben-Younes, Patrick Pérez, and Matthieu Cord,  “Octet: Object-aware counterfactual explanations,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2023, pp. 15062–15071.
  * [12] Lisa Anne Hendricks, Ronghang Hu, Trevor Darrell, and Zeynep Akata,  “Grounding visual explanations,”  in Proceedings of the European conference on computer vision (ECCV), 2018, pp. 264–279.
  * [13] Shunsuke Tsubaki1 Yohei Kawaguchi, Tomoya Nishida, Keisuke Imoto, Yuki Okamoto, Kota Dohi, and Takashi Endo,  “Audio-change captioning to explain machine-sound anomalies,”  in Proceedings of the DCASE 2023 Workshop, 2023.
  * [14] Daiki Takeuchi, Yasunori Ohishi, Daisuke Niizumi, Noboru Harada, and Kunio Kashino,  “Audio difference captioning utilizing similarity-discrepancy disentanglement,”  in Proceedings of the DCASE 2023 Workshop, 2023.
  * [15] Emre Kıcıman, Robert Ness, Amit Sharma, and Chenhao Tan,  “Causal reasoning and large language models: Opening a new frontier for causality,”  arXiv preprint arXiv:2305.00050, 2023.
  * [16] Judea Pearl,  Causality,  Cambridge university press, 2009.
  * [17] Benjamin Elizalde, Soham Deshmukh, Mahmoud Al Ismail, and Huaming Wang,  “Clap learning audio concepts from natural language supervision,”  in ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2023, pp. 1–5.
  * [18] Konstantinos Drossos, Samuel Lipping, and Tuomas Virtanen,  “Clotho: An audio captioning dataset,”  in ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2020, pp. 736–740.
  * [19] Irene Martín-Morató and Annamaria Mesaros,  “What is the ground truth? reliability of multi-annotator data for audio tagging,”  in 2021 29th European Signal Processing Conference (EUSIPCO). IEEE, 2021, pp. 76–80.
  * [20] Chris Dongjoo Kim, Byeongchang Kim, Hyunmin Lee, and Gunhee Kim,  “Audiocaps: Generating captions for audios in the wild,”  in Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), 2019, pp. 119–132.
  * [21] Jason Wei, Xuezhi Wang, Dale Schuurmans, Maarten Bosma, Fei Xia, Ed Chi, Quoc V Le, Denny Zhou, et al.,  “Chain-of-thought prompting elicits reasoning in large language models,”  Advances in Neural Information Processing Systems, vol. 35, pp. 24824–24837, 2022.
  * [22] Joshua Robinson, Ching-Yao Chuang, Suvrit Sra, and Stefanie Jegelka,  “Contrastive learning with hard negative samples,”  arXiv preprint arXiv:2010.04592, 2020.
  * [23] Qiuqiang Kong, Yin Cao, Turab Iqbal, Yuxuan Wang, Wenwu Wang, and Mark D Plumbley,  “Panns: Large-scale pretrained audio neural networks for audio pattern recognition,”  IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 28, pp. 2880–2894, 2020.
  * [24] Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac, Tim Rault, Rémi Louf, Morgan Funtowicz, et al.,  “Transformers: State-of-the-art natural language processing,”  in Proceedings of the 2020 conference on empirical methods in natural language processing: system demonstrations, 2020, pp. 38–45.
  * [25] Karol J Piczak,  “Esc: Dataset for environmental sound classification,”  in Proceedings of the 23rd ACM international conference on Multimedia, 2015, pp. 1015–1018.
  * [26] Justin Salamon, Christopher Jacoby, and Juan Pablo Bello,  “A dataset and taxonomy for urban sound research,”  in Proceedings of the 22nd ACM international conference on Multimedia, 2014, pp. 1041–1044.
  * [27] Tianyi Zhang, Varsha Kishore, Felix Wu, Kilian Q Weinberger, and Yoav Artzi,  “Bertscore: Evaluating text generation with bert,”  arXiv preprint arXiv:1904.09675, 2019.
  * [28] Ho-Hsiang Wu, Prem Seetharaman, Kundan Kumar, and Juan Pablo Bello,  “Wav2clip: Learning robust audio representations from clip,”  in ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2022, pp. 4563–4567.
  * [29] Andrey Guzhov, Federico Raue, Jörn Hees, and Andreas Dengel,  “Audioclip: Extending clip to image, text and audio,”  in ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2022, pp. 976–980.


# Disentangling Multiple Emitting Components in Molecular Observations with
Non-negative Matrix Factorization

Damien de Mijolla Department of Physics and Astronomy, University College
London, Gower Street, London WC1E 6BT  Jonathan Holdship Leiden Observatory,
Leiden University, PO Box 9513, NL-2300 RA Leiden, The Netherlands  Department
of Physics and Astronomy, University College London, Gower Street, London WC1E
6BT  Serena Viti Leiden Observatory, Leiden University, PO Box 9513, NL-2300
RA Leiden, The Netherlands  Department of Physics and Astronomy, University
College London, Gower Street, London WC1E 6BT  Johannes Heyl Department of
Physics and Astronomy, University College London, Gower Street, London WC1E
6BT

###### Abstract

Molecular emission from the galactic and extragalactic interstellar medium
(ISM) is often used to determine the physical conditions of the dense gas.
However, even from spatially resolved regions, the observed molecules are not
necessarily arising from a single component. Disentangling multiple gas
components is often a degenerate problem in radiative transfer studies. In
this paper we investigate the use of the non-negative matrix factorization
(NMF) approach as a means to recover gas components from a set of blended line
intensity maps of molecular transitions which may trace different physical
conditions. We run a series of experiments on synthetic datasets designed to
replicate conditions in two very different environments: galactic pre-stellar
cores and the ISM in high redshift galaxies. We find that the NMF algorithm
often recovers the multiple components resembling those used in the data-
generating process, provided that the different components have similar column
densities. When NMF fails to recover all the individual components it does
however group together the most similarly emitting ones. We further found that
initialisation and regularisation are key factors in the efficiency of the NMF
algorithm.

Molecular gas (1073), Interstellar phases (850), Submillimeter astronomy
(1647), Astrostatistics techniques (1886), Astrochemistry (75)

##  1 Introduction

Understanding the chemical and physical processes regulating the interstellar
medium (ISM) is necessary for a thorough understanding of many astronomical
processes, such as star formation (van Dishoeck & Blake, 1998) and galaxy
quenching (Colombo et al., 2020). As molecular lines are the primary means
through which the conditions of the interstellar medium can be probed, their
interpretation is of vital importance for accomplishing these goals.

Radiative-transfer modelling is the framework through which radiation emitted
by the ISM can be connected back to the physical conditions of the gas it is
emitted from. Through computational radiative-transfer modelling codes, such
as RADEX (van der Tak et al., 2007), the strength of molecular transitions can
be used to constrain the chemistry, temperature, density and other
astrophysical conditions of the emitting interstellar medium. This, in turn,
then enables a better understanding of the local environment in which the ISM
is embedded.

Unfortunately, the radiative-transfer modelling problem needing solving when
analysing molecular lines is typically degenerate. Because of the large number
of free parameters in radiative-transfer models and their similar influence on
intensities, usually multiple set of input parameters are capable of
reproducing any given set of molecular line intensities. In Tunnard & Greve
(2016), it was found that radiative transfer modeling was typically incapable
of recovering parameters to better than half a dex. Although such degeneracies
can be partially addressed through using rigorous statistical methodologies
(Tunnard et al., 2015; James et al., 2021) or through incorporating
astrochemical knowledge (Viti, 2017; de Mijolla et al., 2019), the
interpretation of molecular lines still remains difficult.

Notably, the analysis of emission arising from blended components is
particularly challenging. In sub-mm astronomy, the spatial resolution limits
of modern telescopes mean that we typically do not observe flux from a single
gas component but rather from multiple non-resolved components, each with
unique physical conditions. This is particularly the case for extra-galactic
observations where objects as large as giant molecular clouds, or even whole
sections of a galaxy, will be contained in a single beam. When this is the
case, interpreting the physical conditions of the ISM becomes especially
difficult as it is unclear what fraction of the measured emission originates
from each component or even how many components are responsible for the
observed emission. Moreover, as observers increasingly turn their attention
towards the analysis of external galaxies, where multiple blended components
are expected, mitigating this blending is likely to become an ever more
pressing concern.

To address the difficulties associated with the radiative-transfer problem,
there is value in approaching the task of disentangling gas components within
observations from a more data-driven perspective. In this study, we
investigate the use of data-driven approaches for separating the flux emitted
by molecular lines into the individual gas components within observations.
More precisely, we investigate the use of NMF (Lee & Seung, 1999) which is a
matrix-factorization algorithm used to decompose a non-negative matrix into a
product of non-negative matrices. In NMF, components are estimated globally
from all pixels available. As there will typically be far fewer components
than pixels, NMF can then exploit the redundancy across pixels to constrain
the possible shapes of the components.

Our work is not the first attempt at using data-driven approaches to
understand the ISM. For example, there has been a long line of research into
using Principal Component Analysis to interpret observations of the ISM (Lo et
al., 2009; Gratier, Pierre et al., 2017; Neufeld et al., 2007). Perhaps more
similar to our work, there have also been efforts into using NMF to interpret
interstellar spectra. For example, NMF has been used to analyze interstellar
bands (Berné et al., 2007; Foschino et al., 2019). Efforts at using NMF in
astronomy date as far back as 1997, where positive matrix factorization, a
precursor to NMF, was used to interpret molecular observations (Juvela et al.,
1996a).

Specifically, in this study we attempt to answer the question of whether NMF
can be used to recover gas components from a set of blended line intensity
maps of transitions tracing different physical conditions. Here, unlike
previous approaches such as those presented in Boulais et al. (2021); Juvela
et al. (1996a), our interest is less about identifying co-moving parcels of
gas within observations and more about recovering components tracing a common
physical environment even if they may not have originated from the same
spatial location. Whereas previous approaches have operated on the spectral
datacube where frequency information is available and only made use of a small
number of line-transitions, here we operate on the derived line-intensity maps
for a larger number of lines and in doing so ignore the spectral shape of
line-transitions. This means that our factorization ignore the radial velocity
of spectral lines and only use the molecular composition when constructing
components. Our hope is then for NMF components to group together all emission
arising from gas under similar physical conditions leading to NMF components
capturing the archetypical signatures of physical environments in
observations.

##  2 Data model

We begin by detailing a model for the blending of components which relates
line intensities and spatial locations of components in a given region to the
intensities measured in line-intensity maps.

Let us represent a set of molecular line intensity-maps, observed over a
region of $N_{x}$ $\times$ $N_{y}$ pixels, by a non-negative matrix $X$ of
shape $N_{M}\times N_{p}$ whose entries contain the intensities for $N_{M}$
molecular line transitions measured at the $N_{p}=N_{x}\times N_{y}$ pixels.
We assume that multiple components contribute additively to the intensities
such that the intensity at every pixel is the sum of the intensities of the
separate components. We represent the strength of emission of the components
in the various line-transitions by a (non-negative) matrix $W_{true}$ of shape
$N_{M}\times N_{C}$ whose $N_{C}$ columns each represent the intensities
emitted in the $N_{M}$ molecular transitions for a component.

Similarly, we encode the spatial contributions of each of the $N_{C}$
components to the $N_{p}$ pixels into a (non-negative) matrix $V_{true}$ of
shape $N_{C}\times N_{p}$. Spatial contributions are factors between 0 and 1
indicating how much a component contributes to a pixel. The intensity of a
given line at a given pixel will simply be the sum of the emission in that
line for all components multiplied by their contribution to that same pixel.
That is to say the line-intensity maps can be written as

$X=W_{true}V_{true}+\Sigma$ (1)

where $\Sigma$ is an additional noise term chosen to be Gaussian.

So far, we have neglected the impact of convolution by a telescope beam factor
on observations. Provided that all line intensities are downsampled to a
common spatial resolution, convolution can be modeled through a matrix
multiplication by a convolution matrix K of shape $N_{p}\times N_{p}$. The
data-generating process then becomes $X=W_{true}V_{true}K+\Sigma$ which can be
written using the associativity property of matrix multiplication as
$X=W_{true}V_{conv}+\Sigma$ where $V_{conv}=V_{true}K$. Hence, including
convolution by a beam-filling factor does not change the original model beyond
the replacement of $V_{true}$ with $V_{conv}$.

A number of conditions must be met for our model of line-intensity maps to be
accurate. Firstly, molecular transitions must not be optically thick so as for
intensities to be (near) additive. Secondly, all components must be purely
emitting rather than absorbing as we assume all contributions are positive.
Finally, as explained in the previous paragraph, transitions must be
preliminarily downsampled to a common spatial resolution for the convolution
to be parameterizable by a matrix $K$.

##  3 Non-negative matrix factorization

In this study, we consider the problem of recovering a set of blended
components $W_{true}$ and their convolved spatial locations $V_{conv}$ from
intensity maps $X$, under the assumption that $X=W_{true}V_{conv}+\Sigma$.
Although many algorithms could be relevant for this task, here we focus on
investigating the non-negative matrix factorization algorithm (NMF) (Juvela et
al., 1996b). Given a matrix $X$, the NMF algorithm aims to find non-negative
matrices $W$ and $V$ such that $X\approx WV$, where $W$ and $V$ are matrices
of shape $N_{M}\times N_{NMF}$ and $N_{NMF}\times N_{p}$ and $N_{NMF}$ is the
number of components used. In the context of molecular observations, $N_{NMF}$
is the number of physically distinct emitting components assumed to exist
within telescope beam. It is a user-defined parameter for the algorithm which
should ideally be equal to $N_{C}$ the true number of components so as for the
dimensions of $W$ and $V$ to match those of $W_{true}$ and $V_{true}$. NMF
falls within the realm of blind-source separation algorithms which is to say
that it is an algorithm making minimal assumptions on the contents of matrices
$W$ and $V$.

Formally, the NMF algorithm proceeds by determining non-negative $W$ and $V$
minimizing a loss function

$L=\|X-WV\|$ (2)

where $\|A\|$ is a matrix-norm which is often chosen to be the Frobenius norm
$\|A\|_{\text{F}}^{2}=\sum_{i,j}A_{ij}^{2}$. Minimizing the loss function is a
non-convex optimization problem which is typically approached using iterative
solvers. For this paper, we use the scikit-learn implementation of NMF making
use of a coordinate-descent solver (Cichocki & Phan, 2009) with $W$ and $V$
initialized as random matrices.

Matrix factorization is an ill-posed problem. Indeed, if $X=WV$, then it would
also be the case that $X=WPP^{-1}V$ for any invertible matrix P (of shape
$N_{C}\times N_{C}$) and so if $W$ and $V$ are solutions to the matrix
factorization problem $WP$ and $P^{-1}V$ will also be solutions, provided they
are non-negative matrices. Because of this, there is no formal guarantees that
the NMF algorithm applied to line-intensity maps $X$ will converge towards
$W_{true}$ and $V_{conv}$, even in the limit of infinite data.

In NMF, to further constrain the matrix factorization problem towards
desirable solutions, it is common to add to the loss function additional
regularization terms on $W$ and $V$. In such cases, the loss function
minimized by NMF becomes:

$L=\|X-WV\|_{F}+\alpha_{1}\|W\|_{R}+\alpha_{2}\|V\|_{R}$ (3)

In addition to regularization, the initialization of matrices $W$ and $V$ can
impact the convergence of NMF and hence the matrices retrieved by the
algorithm. In this project, we initialize the algorithm with non-negative
random matrices as defined in the scikit-learn NMF documentation.

It is also worth mentioning that the NMF algorithm, like many other matrix
factorization algorithms, is most effective when there are significantly fewer
components than molecular lines; that is to say when $N_{C}<<N_{M}$. This is
because NMF relies on finding a set of components capturing the most
information about observations X which requires grouping together correlated
molecular lines. In the limit where the number of components matches the
number of molecular lines observed ($N_{C}=N_{M}$), the NMF loss is trivially
minimized by setting $W$ to the identity matrix (i.e. one component
reconstruct each line) and $V$ to be equal to $X$. Such a solution perfectly
reconstructs observations but brings no useful information about the
correlations amongst molecular lines.

##  4 Synthetic data generation

In order to assess the use of NMF in interpreting molecular line-intensity
maps, we test its effectiveness at retrieving components from synthetic line-
intensity maps. Synthetic data allow us to have exact knowledge of the ground-
truth components and therefore evaluate how closely they are recovered by NMF.
Our aim when constructing these synthetic observations is to approximately
reproduce components and observations found within real-world datasets. Our
aim is not to build extremely accurate simulations. As such, we have striven
for simplicity over complexity in our data-generating process.

We construct synthetic line-intensity maps of transitions of CO, CS, HCN and
HCO+ (see Section 5), using numpy (Harris et al., 2020), by first using
chemical and radiative transfer models to produce emission profiles
($W_{true}$) from chosen physical conditions. We then specify an arbitrary
spatial location of the emission of each component ($V_{true}$). Finally, we
select a noise level ($\Sigma$) and convolution due to the beam-filling factor
($K$) to produce the maps ($X$) following Eq 1 as discussed in Section 2.

The matrix $W_{true}$ contains the emission of each component for the
molecular lines considered in $X$. For our experiments, we select a
temperature and density for each component and the intensity of molecular
lines are then determined using the radiative transfer modelling code RADEX
(van der Tak et al., 2007). RADEX requires a gas density and temperature as
well as the column density of the emitting species. Since the emission is
proportional to the column density, and we choose the noise level, the
absolute value of the column density is completely arbitrary and has no effect
on the experiments. However, the differences in column density between
components is observationally important so we use an astrochemical model to
obtain reasonable values.

To do this, we use the astrochemical code UCLCHEM (Holdship et al., 2017) to
find the steady state abundances of our species for the chosen temperature and
density of each component, assuming standard Milky Way values for the cosmic-
ray ionization rate and UV field of, respectively, $1.3\times 10^{-17}\rm
s^{-1}$ and 1 Draine (equivalent to 2.74$\times$10-3 erg/s/cm2). We convert
abundances to molecular column densities by finding the factor that would
convert the component with the highest CO abundance to a CO column density of
$10^{19}\rm cm^{-2}$. We then scale the other components by the same factor so
that their column density ratios are equal to their abundance ratios.
Therefore, whilst our absolute column density values are arbitrary, their
ratios are informed by chemistry.

The matrix $V_{true}$ encodes how the emission of each component varies
spatially from the peak value given by RADEX. We model the spatial
contribution of a component at every pixel as a 2d Gaussian whose dimensions
correspond to the x and y coordinates of the pixel. That is to say, the
contribution at a pixel $v_{xy}$ is $v_{xy}=A\times N(\mu,\sigma)$ evaluated
at the x,y coordinates of the pixel center. In this data-generating process,
the mean $\mu$ of a component controls its location, the covariance $\sigma$
its extent and the amplitude $A$ its relative strength. For all components,
the amplitude A is set such that $v_{xy}=1$ at the peak location of the
component. The entries of $V_{true}$ are then the contribution $v_{xy}$ of the
$N_{p}$ pixels for the $N_{D}$ components.

The telescope beam-response is approximated by a Gaussian kernel K. After
creating noiseless line-intensity maps according to $W_{true}V_{true}K$, in
our final step gaussian noise N is added to line-intensity maps such that
every map has a signal-to-noise-ratio of 100 with gaussian noise
independently, identically, and uniformly added to all pixels.

##  5 Experiments

We now present our experiments on applying NMF to synthetic observations (with
results plotted using matplotlib (Hunter, 2007)). All experiments are carried
out using a common blueprint for the characteristics/properties of the line-
intensity map. We choose to produce synthetic observations of two
environments: protostellar cores and high redshift galaxies. We arbitrarily
pick the molecular transitions from the observational study of the nearby
galaxy NGC 1068 from Viti et al. (2014). While the choice of this particular
sample of lines was arbitrary, we note that it contains some of the most
observed molecules in both galactic as well as extragalactic environments. All
line-intensity maps are constructed on a 50 $\times$ 50 dimensional grid (i.e.
$50\times 50$ pixels), for the transitions observed in Viti et al. (2014), to
which Gaussian noise is added such that maps have a mean per-pixel SNR of 100
and to which a convolution with a Gaussian kernel of width 1/10 of the overall
width of maps has been applied. Although this blueprint is only one amongst
many possible parametrizations, we found the conclusions of our study to be
robust to changes to the blueprint.

In the following experiments, before application of NMF, we always apply a
series of preprocessing steps to the line intensity maps. These preprocessing
steps would likely need to be reproduced when working on real observations. As
a first step, in order to satisfy the non-negativity requirement of NMF, all
negative-entries in the line-intensity maps are zeroed-out. As a second step,
so as for the NMF loss to equally weight all line-intensity maps, we rescale
line-intensity maps prior to running NMF through multiplication by a scaling
factor such that, after multiplication, all line-intensity maps have an equal
average flux of unity. After the NMF algorithm has finished running, this
scaling factor is erased from the NMF components through dividing entries in
$W$ by the scaling factor associated to their transition.

As described in Section 3, there are two regularisation terms that can be used
to control how the algorithm converges. To investigate the sensitivity of our
approach to this regularisation, we try 3 different regularisation approaches
(regularise $W$ on its own, $V$ on its own, or both $W$ and $V$
simultaneously) and for each of these regularisation approaches we try 21
different values of $\alpha$.

Whilst testing the sensitivity to regularisation is important, there is a
further confounding factor. The result of NMF varies depending on the
initialisation of matrices $W$ and $V$. To capture this variation, we repeat
each regularisation set up 300 times with different random matrix
intialisations which is achieved by passing init=“random” in the scikit-learn
NMF implementation. This allows us to be sure that where one regularisation
set up performs better than another, it is not due to random fluctuations
because we can test whether it is consistent across 300 repeats.

Overall, each experiment is performed with 3 different regularisation
approaches at 21 different regularisation strengths repeated 300 times each
for a total of 18,900 runs. To analyse this population of NMF runs, we have
created a goodness-of-fit metric condensing all information about each
individual run into a single number quantifying the level of agreement between
recovered and true components. This metric is defined as the mean-absolute
error (MAE) between the true components emission profiles $W_{true}$ and those
retrieved by the NMF run ($W$). However, since components returned by NMF can
only be correct up to a scaling factor, we rescale the emission profiles to
have an average flux (averaged over all transitions) of unity before
calculating the MAE. This definition of MAE means that the line intensities
($W$) of components retrieved by a run with a MAE of 0.1 will on average
differ from the ground-truth intensities ($W_{true}$) by $10\%$ of the average
line-intensity of their closest matching ground-truth component. Finally,
since any of the three NMF components could correspond to any of the three
true components, we choose the permutation that gives the lowest MAE. We show
in Appendix 6.1 the pseudo-code summarising the overall procedure for
generating results as described in this section.

###  5.1 Protostellar environments

Figure 1: Synthetic line-intensity maps of molecular clouds ($X_{proto}$ :
bottom) alongside the emission profile ($W_{true}$ : top-right) and convolved
(unitless) spatial contribution ($V_{conv}$ : top-left) of each components
contributing to the maps.

In this section we present insights derived from applying NMF to the first of
two sets of mock line-intensity maps $X_{proto}$, which models emission from
blended protostellar environments. The maps in $X_{proto}$ are generated from
three components, each meant to approximate a gas component found in star
forming regions. Generated components are i) An extended “low density”
component (T=20 K and n=$10^{4}\rm cm^{-3}$), ii) a “medium density” component
(T=10 K and n=$10^{5}\rm cm^{-3}$), and iii) a “dense” warmer component (T=100
K and n=$10^{7}\rm cm^{-3}$). These three idealized components may, for
example, be representing the surrounding ISM, the envelope, and the more
central compact region of a protostellar core.

The synthetic-line intensity maps ($X_{proto}$) as well as emission profiles
and convolved contributions of emitting components are shown in Figure 1. From
this Figure, we see some moderate spatial overlap between components as well
as some components emitting more strongly than others. Both of these factors
make it difficult, at least visually, to recover components from the line-
intensity maps motivating the use of specialized retrieval algorithms.

We show in Figure 2 outputs obtained from running a population of NMF runs on
$X_{proto}$. Runs in the population use three NMF components, the same number
as the number of chemical components used in the modelling, which means that
the algorithm should have enough capacity to fit all components. Figure 2a (on
the left) characterises the goodness of fit of NMF across this population of
runs. Each panel of the figure represents a binned density plot of runs in
which the x-axis coordinates quantify the amplitude of the regularisation term
and the y-axis coordinate the goodness of fit of the retrieved components as
defined by the MAE metric introduced in Section 5. The binned density plot
captures the probability distribution of the goodness-of-fit metric at a given
regularisation value (i.e. at a fixed x-axis value). Each panel can be viewed
as a succession of vertically displayed histograms capturing the probability
distribution in the MAE metric for a regularisation level specified by the
axis. Here, the probability distributions capture the variability which arises
due to the random initialisation for the initial guesses of matrices $W$ and
$V$. Each panel considers a different type of regularisation: i) only on $W$
(top panel), ii) only on $V$ (central panel), or iii) on both $W$ and $V$
(bottom panel). Figure 2 b and c focus-in on two runs within the population
and compare the relative line strengths ($W$) and spatial contributions ($V$)
of these runs to those of the ground-truth components. The locations of these
runs within the density plots are shown by markers #1 and #2. These two runs
are selected for further study from the population of NMF runs because they
represent two extremes (in terms of agreement with $W_{true}$) amongst all
runs.

We find an excellent agreement between the components from run #1, one of the
better fitting runs in the population, and the true components. Although there
are minor mismatches between true and recovered components, such as most of
the HCN(1-0) originating from the diffuse (top-panel) true component being
mistakenly attributed to the middle-panel retrieved component, the recovered
components are by and large in excellent agreement with the true gas
components. However other runs in the population differ more strongly from the
true components. Run #2 is an extreme example of such a run with one of the
lowest y-values across the population. This run exhibits a much higher level
of disagreement with the true components, with in particular its top-panel
component not matching well any of the true components. The disagreement
observed stems from the run learning a different decomposition into components
than that used to generate the data. In this decomposition, the top panel
component captures emission from both the top and bottom panel of the true
components. Such a decomposition, while different from the data-generating
decomposition, was found to match the lower MAE runs (e.g. run#1) in terms of
its quality of fit to the line-intensity maps. The decomposition found in
run#2 is thus an equally valid decomposition that cannot be ruled-out from
studying the quality of the match between $X_{proto}$ and the factorized
approximation $WV$.

By comparing the distribution of the y-coordinates of runs in the models
population to the y-values of runs #1 and #2, we can understand how effective
and reliable NMF is at recovering components closely matching with the
underlying component. We observe that the vast majority of runs, irrespective
of regularisation, have y-values comparable to #1, indicating that most runs
retrieve components matching well with the underlying components. While poor-
match runs comparable to #2 do exist, these runs are relatively rare. This
seems to suggest that, at least for the considered dataset, NMF is an
effective, albeit imperfect, algorithm for recovering the underlying
components.

Figure 2: a) Binned density plots characterizing the performance of a
population of three-component NMF runs on $X_{proto}$ in which x-axis
coordinates quantify the amplitude of the regularization term and y-axis
coordinates the goodness of fit of the retrieved components as defined by the
MAE metric introduced in Section 5. Each panel considers a different type of
regularization: i) only on $W$ (top panel), ii) only on $V$ (central panel),
or iii) on both $W$ and $V$ (bottom panel). b) Emission profiles of the true
components alongside those of two runs in the population whose location in the
scatter-plot are represented by markers #1 & #2. c) Convolved spatial
contributions of the same two runs.

We can also study the impact of hyperparameters on NMF results. Regularization
appears to have a beneficial effect on the agreement between recovered and
ground-truth components with even small amounts of regularization on $W$ being
sufficient to steer the factorization away from the higher MAE solutions.
Whilst regularization on $V$ does not have as strong of an impact, it also
improves the match between true and recovered gas components as can be seen
from how the best performance is obtained with high-levels of regularization
on both $W$ and $V$. For now we do not discuss how the choice of the number of
NMF components affects results (but such a discussion can be found in Section
5.3). However, it is worth acknowledging that we expect performance to degrade
when using a sub-optimal number of components.

###  5.2 High-z galaxy (three-component fits)

Disentangling emission arising from blended components is of course a problem
that exacerbates the further away the object is. In this section, we therefore
study the performance of the NMF algorithm on a synthetic dataset
approximating line-intensity maps from a high-redshift galaxy $X_{gal}$.
Unlike in the previous experiment, where components had only moderate spatial
overlap, in this scenario all components peak in the same location - the
center of the galaxy - but have different spatial extents. This high-overlap
between components complicates the retrieval process as there are fewer pixels
which can clearly be attributed back to a unique component. This dataset thus
constitutes a more challenging test for NMF but an important one as the NMF
algorithm will only be practical if it can reliably produce useful retrievals
across a range of different scenarios.

The components used in this run are a less dense component (T=10 K and
n=$10^{4}\rm cm^{-3}$), a medium density moderately extended and warm
component (T=30 K and n=$10^{5}\rm cm^{-3}$), and a dense compact hotter
central component (T=50 K and n=$10^{7}\rm cm^{-3}$). We note that, for a
spatially unresolved high redshift galaxy, our three gas components are
definitely more arbitrary in nature than in our previous example.
Nevertheless, we may identify the first component with the average cold
molecular gas, typical of a Giant Molecular Cloud, the second component with
dense molecular gas heated by either nearby star formation or an AGN, and the
third component with compact very dense gas found in the proximity of the
nucleus of a galaxy.

Figure 3: Synthetic line-intensity maps of high-redshift galaxy ($X_{gal}$ :
bottom) alongside the emission profile ($W_{true}$ : top-right) and convolved
(unitless) spatial contribution ($V_{conv}$ : top-left) of each components
contributing to the maps. Figure 4: a) Binned density plots characterizing the
performance of a population of three-component NMF runs on $X_{gal}$ in which
x-axis coordinates quantify the amplitude of the regularization term and
y-axis coordinates the goodness of fit of the retrieved components as defined
by the MAE metric introduced in Section 5. Each panel considers a different
type of regularization: i) only on $W$ (top panel), ii) only on $V$ (central
panel), or iii) on both $W$ and $V$ (bottom panel). b) Emission profiles of
the true components alongside those of two runs in the population whose
location in the scatter-plot are represented by markers #1 & #2. c) Convolved
spatial contributions of the same two runs.

Once again, we validate the algorithm through plotting i) line-intensity maps
and associated components, ii) density plots showing the MAEs for a population
of three-component runs, and iii) the retrieved components for two runs whose
locations in the scatter plot are shown by markers #1 and #2. These plots are
shown in Figure 3 and 4. NMF runs use three NMF components to fit the line-
intensity maps which is a number equal to the number of components within
observations and as many as was used in the previous experiment.

We can see from these figures that there is a relatively good level of
agreement between the true components and those recovered in runs #1 and #2,
with runs having at most a MAE of 0.2 (which represents an average of 20%
difference in intensity between the true line intensities and those produced
by the NMF algorithm). This suggests that for this high-overlap scenario, NMF
can approximately recover the underlying components. However, the mismatch
between true and recovered components appears larger here (at least in terms
of contributions) than in the previous low-overlap scenario (Section 5.1). In
particular, the algorithm seems to struggle with decomposing the fully
overlapping central component, favouring a decomposition into spatially non-
overlapping components, despite the high overlap betweenthe underlying
components. This can be seen by how, in both runs, emission arising from more
extended components in the central portion of the observations is
misattributed to the less extended retrieved component (bottom panel). This in
turn then leads to a mismatch between the true and retrieved emission profile
for this less extended components as well as the creation of donut-like holes
in the spatial profiles of the other more extended components (at the location
of misattributed emission).

Overall, these results highlight the presence of degeneracies when doing
retrievals of observations containing highly overlapping components. When such
high-spatial overlap components exist, the NMF algorithm struggles to recover
the proper overlap between components. Furthermore, no choice of
regularisation appears capable of lifting these degeneracies. Nevertheless,
recovered components still match fairly well with observations.

###  5.3 High-z galaxy (two-component fits)

In real-world applications, the number of components used in the NMF retrieval
will need to be estimated from observations. However, because of differences
in physical conditions and chemical make-up, different spatial locations will
emit in different molecular transitions leading to the existence of a higher
number of components. In such cases, to maintain the regime of fewer NMF
components than observed line-transitions, NMF retrievals must be run with
fewer NMF components than required to fully capture observations. The aim is
then not to retrieve the true components but instead to retrieve useful
components. We examine here how NMF behaves in this regime.

To do this, we run NMF on our high-z line intensity maps but using only two
components instead of the three components actually required for reproducing
the line-intensity maps. This provides a test for the behaviour of the NMF
algorithm when data is more complex than the available number of components.
Figure 5 shows the output of such two-component runs. Here, we do not show
density plots since our goodness of fit metric is not applicable when the
number of retrieved and ground truth gas components differ.

Comparing these runs to their three-component counterparts (Figure 4) reveals
that the two-component runs seem to only differ from the three-component runs
in their treatment of the top and middle-panel components. Because the two-
component runs are incapable of fully capturing the emission from all three
components, they group together the top and middle-panel components and
instead return a component whose emission profile is approximately a mixture
of the emission from these two components. On the other hand, the lowest-panel
component of the two-component runs, characterizing the densest component, is
almost identical to the analogous component in the three-component runs. In
fact, similarly to the three-component runs, it (mistakenly) captures emission
from other gas-phases within the central region leading to a similar mismatch
between retrieved and true dense component and a similar (incorrect) donut
shape spatial emission in other components.

Conceptually, it appears that because the two-component NMF runs did not have
enough capacity to reproduce all three components, they selectively grouped
the two most similar components - the top and middle panel components - into a
single “archetypical component” capturing the average emission across both
phases whilst dedicating the remaining component to the densest component
because it emitted more uniquely. Such behaviour suggests that, when lacking
enough components to fit observations, the NMF algorithm will fall-back on
grouping together the most similar components into archetypical components.
This means that the algorithm could allow for synthesizing highly-complex
observations into a manageable number of components for human interpretation.

Figure 5: Emission profiles and (unitless) convolved spatial contributions of
the true components alongside those of two random runs in a population of two-
component runs fitted to $X_{gal}$.

###  5.4 High-z galaxy - dependence on the number of components

In the previous section we studied the fit to our high-z environment
$X_{gal}$, retrieved by a two-component NMF model. We found that the two-
component model did a good job at reproducing line-intensity maps despite
these technically containing three components. In this section we explicitly
study the effect of varying the number of components $N_{C}$ used in the NMF
algorithm on the quality of fit to line-intensity maps. It is worth noting
that whereas previous sections have evaluated the NMF algorithm by assessing
the quality of its retrieved components, here we evaluate NMF by assessing how
well it reproduces the line-intensity maps (irrespective of whether this
emission arises from components resembling the true gas phases). This better
matches the reality of fitting an unknown number of components where the only
data we have are the line-intensity maps.

To show this, in Figure 6 we plot the mean-squared error (MSE) between the NMF
fits to $X_{gal}$ (obtained as $WV$) and $X_{gal}$ for different values of
$N_{C}$, the number of components in the NMF (shown on the x-axis). Here, to
calculate the MSE, we first normalise each line-intensity map in order to
prevent any single line-intensity map from dominating over others in the
calculation of the MSE, then flatten the line-intensity maps into a 1D vector
and finally calculate $MSE=\frac{1}{N}\sum_{i=1}^{N}(x_{i}-\hat{x}_{i})^{2}$.
We also show, as a baseline, the MSE between $X_{gal}$ and the noiseless
version of $X_{gal}$ which represents some ideal MSE value for a model
reproducing the ground-truth gas phases and none of the noise.

Interestingly, we find that a two-component NMF fit, similar to that explored
in the previous subsection, does a good job at matching the emission in
observations $X_{gal}$ despite using fewer components than the number of gas
phases present in the observations. This is particularly interesting because
it may be conceptually similar to analysing an emission map as a set of
archetypical components (gas originating from a starburst ring; dense gas
surrounding an AGN; etc) as it is commonly done, rather than considering its
true physical gas components. In fact, we find that the two-component fit has
a lower MSE than the baseline representing the MSE solely due to noise. This
is due to the NMF decomposition overfitting the noise in the line-intensity
maps. From this figure, we further find that increasing the number of
components in the NMF algorithm leads to further overfitting to the noise in
observations, as evidenced by the seemingly linear decrease in MSE with
increasing dimensionality of $N_{C}$.

Since the noise level can be estimated from the map itself, this provides a
method for estimating the number of components one should use to analyse
observations. A similar plot can be produced for a real line-intensity map and
the optimal number of components is the number that causes the MSE of the NMF
algorithm to reach the true noise level. Crucially, this would be the number
of gas components that are sufficiently distinct for the algorithm to
separate, and not the true number of physical components. However, this would
be useful since, even when analysing archetypical components, there is often a
question of how many to include. This approach should be further tested by
analysing low resolution maps of objects where high resolution imaging has
uncovered distinct components.

Figure 6: Mean-squared error (MSE) between the NMF fits to $X_{gal}$ (obtained
as $WV$ and denoted as $\hat{X}(N_{C})$) and $X_{gal}$ for varying number of
NMF components $N_{C}$ (blue) where MSE is plotted on log-scale. As a baseline
we also plot the MSE between $X_{gal}$ and the noiseless version of $X_{gal}$
(green).

###  5.5 Discussion and conclusions

We have run a series of experiments aiming at recovering the underlying
components within observations of molecular line-intensity maps. These
experiments were run on a set of synthetic datasets designed to replicate
conditions expected from real observations while also capturing some
diversities of environments where the algorithm could be applicable. Because
NMF ignores the spatial locations of pixels and would perform equally well if
pixels were scrambled, we did not put a strong emphasises on reproducing the
geometrical structure of real molecular maps and instead focused on
reproducing the expected overlap between components. To capture the diversity
of environments we built two synthetic line-intensity maps on which to test
our observations: one with a moderate overlap between components and a second
dataset with a high overlap. We validated the effectiveness of NMF on these
synthetic environments through studying the algorithm’s ability to recover the
components used in the data-generating process. Because the convergence of the
NMF algorithm was found to be sensitive to hyperparameters, we have
characterised the performance across a range of hyperparameters rather than
for a single set of hyperparameters.

The results of our investigations were as follows. We found that the NMF
algorithm often recover components resembling those used in the data-
generating process. We found that when the NMF algorithm lacked the capacity
to reproduce all components, it would group together the most similarly
emitting components into one archetypical component. However, in general, the
matrix-factorization problem solved by NMF did not have a unique solution and
so its convergence to a set of components closely resembling the underlying
components was not guaranteed. Instead, initialisation and regularisation
played a crucial role towards influencing the type of solution the algorithm
converged towards, with sometimes very different solutions depending on the
hyperparameters. This sensitivity of NMF to model hyperparameters is obviously
not ideal as it means that the algorithm can sometimes recover components that
match poorly with observations. In general we also found that the NMF
algorithm struggled with recovering the proper spatial overlap of overlapping
components.

Our experiments provide a proof-of-concept towards the applicability of NMF
for analysing line-intensity maps. However, some differences do still exist
between real data and the data within our experiments. One such difference is
that our mock line-intensity maps shared a common beam size. When working with
real data, it will be crucial to down-sample observations so as to guarantee
that parcels of emitting gas have the same spatial extent across all line
transition maps. Another difference is that because of effects such as optical
depth and differences in local physical conditions, it is unlikely that real
data will be perfectly characterised as having arisen from a small number of
components. Nevertheless, as demonstrated by our experiment in which we used a
smaller number of NMF components than the number of components, even when the
data-generating process is not perfectly satisfied by observations, NMF still
retains the ability to uncover meaningful components in the form of
archetypical components.

A more important weakness of this work is that our experiments only considered
scenarios in which observations were created from components of similar column
densities and thus in which all components contributed more or less equally to
the line-intensity maps. If components have dissimilar column densities then
it becomes much more challenging for NMF to recover the fainter components.
This is because as the lower column density components emit noticeably more
weakly, correctly fitting them does not contribute much to the minimisation of
the loss function. For synthetic observations containing low-column density
components, in order to recover the weakest components it is crucial that the
signal-to-noise within observations be high enough for these weakest
components to not be drowned by noise.

We also wish to emphasise that our work provides only one, amongst potentially
many, approaches for incorporating matrix factorization constraints into the
component retrieval process. A potential extension of this work would be to
combine NMF with radiative-transfer codes to enforce on retrieved components
both the matrix factorization constraints of NMF and the physical realism
constraints of radiative-transfer codes.

Regardless of the exact approach taken, the large degeneracies associated with
radiative transfer currently make it difficult to constrain the properties of
the ISM especially at high-redshift. Combining the information available
across all spatial locations during inference, as is done in NMF and other
matrix factorization approaches, is, we believe, a promising approach to
address such degeneracies. However, whilst NMF is useful for understanding the
relationship between lines within observations, on its own it is not
sufficient for fully constraining components.

## Acknowledgements

The authors thank the anonymous referee for their constructive comments which
significantly improved the paper. DdM and J. Heyl are funded by an STFC
studentship in Data-Intensive Science (grant number ST/P006736/1). This work
was supported by the European Research Council (ERC) Advanced Grant MOPPEX
833460.

## References

  * Berné et al. (2007) Berné, O., Joblin, C., Deville, Y., et al. 2007, A&a, 469, 575, doi: 10.1051/0004-6361:20066282
  * Boulais et al. (2021) Boulais, Berné, O., Faury, G., & Deville, Y. 2021, A&a, 647, A105, doi: 10.1051/0004-6361/201936399
  * Cichocki & Phan (2009) Cichocki, A., & Phan, A.-H. 2009, IEICE Transactions on Fundamentals of Electronics Communications and Computer Sciences, 92, 708, doi: 10.1587/transfun.E92.A.708
  * Colombo et al. (2020) Colombo, D., Sanchez, S. F., Bolatto, A. D., et al. 2020, A&a, 644, A97, doi: 10.1051/0004-6361/202039005
  * de Mijolla et al. (2019) de Mijolla, D., Viti, S., Holdship, J., Manolopoulou, I., & Yates, J. 2019, A&a, 630, A117, doi: 10.1051/0004-6361/201935973
  * Foschino et al. (2019) Foschino, Berné, O., & Joblin, C. 2019, A&a, 632, A84, doi: 10.1051/0004-6361/201935085
  * Gratier, Pierre et al. (2017) Gratier, Pierre, Bron, Emeric, Gerin, Maryvonne, et al. 2017, A&a, 599, A100, doi: 10.1051/0004-6361/201629847
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Holdship et al. (2017) Holdship, J., Viti, S., Jiménez-Serra, I., Makrymallis, A., & Priestley, F. 2017, The Astronomical Journal, 154, 38, doi: 10.3847/1538-3881/aa773f
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * James et al. (2021) James, T. A., Viti, S., Yusef-Zadeh, F., Royster, M., & Wardle, M. 2021, The Astrophysical Journal, 916, 69, doi: 10.3847/1538-4357/abfd99
  * Juvela et al. (1996a) Juvela, M., Lehtinen, K., & Paatero, P. 1996a, MNRAS, 280, 616, doi: 10.1093/mnras/280.3.616
  * Juvela et al. (1996b) —. 1996b, MNRAS, 280, 616, doi: 10.1093/mnras/280.3.616
  * Lee & Seung (1999) Lee, D. D., & Seung, H. S. 1999, Nature, 401, 788
  * Lo et al. (2009) Lo, N., Cunningham, M. R., Jones, P. A., et al. 2009, MNRAS, 395, 1021, doi: 10.1111/j.1365-2966.2009.14594.x
  * Neufeld et al. (2007) Neufeld, D. A., Hollenbach, D. J., Kaufman, M. J., et al. 2007, ApJ, 664, 890, doi: 10.1086/518857
  * Tunnard & Greve (2016) Tunnard, R., & Greve, T. R. 2016, The Astrophysical Journal, 819, 161, doi: 10.3847/0004-637x/819/2/161
  * Tunnard et al. (2015) Tunnard, R., Greve, T. R., Garcia-Burillo, S., et al. 2015, The Astrophysical Journal, 815, 114, doi: 10.1088/0004-637x/815/2/114
  * van der Tak et al. (2007) van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, Astronomy & Astrophysics, 468, 627, doi: 10.1051/0004-6361:20066820
  * van Dishoeck & Blake (1998) van Dishoeck, E. F., & Blake, G. A. 1998, Annual Review of Astronomy and Astrophysics, 36, 317, doi: 10.1146/annurev.astro.36.1.317
  * Viti (2017) Viti, S. 2017, Astronomy & Astrophysics, 607, A118, doi: 10.1051/0004-6361/201628877
  * Viti et al. (2014) Viti, S., García-Burillo, S., Fuente, A., et al. 2014, Astronomy & Astrophysics, 570, A28, doi: 10.1051/0004-6361/201424116

##  6 Appendix

###  6.1 Pseudo-code