# Intent Detection and Slot Filling for Home Assistants: Dataset and Analysis
for Bangla and Sylheti

Fardin Ahsan Sakib, A H M Rezaul Karim, Saadat Hasan Khan, Md Mushfiqur Rahman
Department of Computer Science, George Mason University 111All authors have
made equal contributions
<EMAIL_ADDRESS>
###### Abstract

As voice assistants cement their place in our technologically advanced
society, there remains a need to cater to the diverse linguistic landscape,
including colloquial forms of low-resource languages. Our study introduces the
first-ever comprehensive dataset for intent detection and slot filling in
formal Bangla, colloquial Bangla, and Sylheti languages, totaling 984 samples
across 10 unique intents. Our analysis reveals the robustness of large
language models for tackling downstream tasks with inadequate data. The
GPT-3.5 model achieves an impressive F1 score of 0.94 in intent detection and
0.51 in slot filling for colloquial Bangla. 222The dataset and the analysis
code can be found in the following directory:
https://github.com/mushfiqur11/bangla-sylheti-snips.git

##  1 Introduction

Smart devices have become commonplace, establishing home assistants as
indispensable fixtures in contemporary households. These voice-activated
virtual companions adeptly manage an array of tasks, ranging from setting
reminders to controlling room temperatures. The efficacy of home assistants in
performing these tasks is closely intertwined with their underlying Natural
Language Understanding (NLU) models, which enable seamless interactions in
high-resource languages Chen et al. (2019); Stoica et al. (2021); Antoun et
al. (2020); Upadhyay et al. (2018). However, this advantage in NLU
capabilities is not extended to low-resource languages Stoica et al. (2019);
Schuster et al. (2018), presenting a notable discrepancy. This discrepancy
holds considerable significance, especially considering the global demand for
home assistants and the extensive usage of low-resource languages, which have
a substantial speaker base.

Bangla and Sylheti Ethnologue (2023), with 285 million native speakers
combined, have rich cultural and colloquial nuances. Specialized datasets are
needed to capture these intricacies as users prefer to interact with home
assistants in their native languages, highlighting the research need Bali et
al. (2019).

The language understanding of home assistants is dependent on two key NLU
tasks: intent detection and slot filling Weld et al. (2022); Louvan and
Magnini (2020). Intent detection determines user actions, like playing music
or checking the weather, while slot filling extracts specific details, such as
song titles or locations. These tasks enable seamless human-device
interactions, especially for home assistants.

Research on intent detection and slot filling primarily focuses on high-
resource languages Liu and Lane (2016); Qin et al. (2021); Niu et al. (2019);
Zhang et al. (2018). While there have been limited studies dedicated to the
Bangla language Bhattacharjee et al. (2021); Alam et al. (2021); Hossain et
al. (2020), none of them have addressed the tasks of intent detection and slot
filling in Bangla. Furthermore, these studies have not taken into account
colloquial variants or closely related languages like Sylheti. This gap in
research leaves a significant portion of the speaker base underserved.

This paper bridges this research gap with several notable contributions.
Firstly, we introduce a comprehensive dataset encompassing 328 entries for
intent detection and slot filling for each of the three languages – totaling
984 samples. These languages include formal Bangla, colloquial Bangla, and
colloquial Sylheti. We further show a comparative study between generative
LLMs and state-of-the-art language models for intent detection and slot
filling.

##  2 Dataset

At the core of our exploration stands a meticulously curated dataset that is
inspired by the SNIPS dataset Coucke et al. (2018), which caters to the broad
audience.

###  2.1 Dataset Size and Distribution

Originating from the 328 English samples present in the SNIPS dataset, our
dataset underwent a manual correction phase to ensure that the English samples
were of optimal quality. Then, we created three linguistically diverse
variants, maintaining the same distribution across intent classes and slots as
the original samples. These are:

  1. 1.

Formal Bangla: This represents the standard version of the Bangla language,
majorly used in contexts like official documents, news broadcasts, and
literature. Formal Bangla tends to adhere strictly to grammatical rules.

  2. 2.

Colloquial Bangla: An informal variant predominantly used in Bangladesh,
colloquial Bangla resonates with everyday conversations of its people. While
there are numerous dialects in different regions of Bangladesh, this form
remains more or less consistent across the country. Colloquial Bangla is more
flexible regarding syntax and incorporates a significant number of loanwords
from English, Arabic, Persian, and other languages.

  3. 3.

Colloquial Sylheti: A language with unique intricacies, Sylheti stands apart
from Bangla and is spoken in the Sylhet region of Bangladesh and among
diaspora communities. It’s rich in expressions, proverbs, and idiomatic
language that reflect the history and culture of the Sylhet region.

The curated dataset spans 10 distinctive intents. Each specific intent has a
distinct set of slot categories. Figure 1 shows the number of samples for each
intent and Figure 2 shows the fraction of slots that frequently occur for each
intent, with respect to infrequently occurring slots.

###  2.2 Data Generation Process

The generation of our dataset was methodical and rigorous to ensure
authenticity and accuracy.
Annotator Engagement
Four doctoral students were on board as annotators for our project. The
initial phase involving the rectification of English data from the SNIPS
dataset was a collaborative effort, with each annotator working on a distinct,
non-overlapping segment. Subsequent phases involved two individuals fluent in
Bangla for the Bangla datasets and two native Sylheti speakers for the
colloquial Sylheti dataset.
Base Creation
The base dataset was created using the Bangla-T5 model Bhattacharjee et al.
(2023), a state-of-the-art English-to-Bangla translation tool, following the
work of De bruyn et al.. The refined English samples served as the foundation
to produce the initial Bangla translations for each sample. An auto-generated
dataset comes with a myriad of issues. Therefore, these samples were manually
re-translated and annotated with the auto-translations as the base.
Inter-Annotator Agreement
An essential step in ensuring the reliability of our dataset was to gauge the
consistency between annotators. For each language variant, 28 randomly chosen
samples were annotated independently by both designated annotators, followed
by calculating their inter-annotator agreement (Table 1). This exercise helped
us discern the degree of concordance and areas of divergence.
Consensus Building
Post the initial agreement calculation, a meeting was convened where the
annotators discussed and reconciled their differences. This step was
instrumental in ironing out inconsistencies and ensuring a unified approach
going forward.
Blind Overlap
As the annotators progressed with data creation, a random 10% of the samples
were earmarked for blind overlap. These served as a secondary check on inter-
annotator agreement after dataset creation.
Independent Adjudication
After the final compilation of the dataset, each entry underwent a rigorous
review by an independent adjudicator who had not previously worked on that
particular language variant. This added an additional layer of scrutiny and
quality assurance.

Inter-annotator agreement
---
|  | Cohen’s
---
Kappa
| Average
---
BLEU
First 28 samples | 0.42 | 0.43
Blind overlap (10%) | 0.55 | 0.51
Table 1: There was an increase in annotator agreement before and after the
annotator’s meeting. This ensures the homogeneity of annotations in the
dataset.

###  2.3 Ensuring Quality

Our data generation process, featuring multiple checks, blind overlaps, third-
party reviews, and inter-annotator agreement stages, highlights our commitment
to quality. It minimizes biases and discrepancies that could result from a
single annotator’s viewpoint. The inclusion of an independent adjudicator in
the final review further bolsters the dataset’s integrity and reliability.
Using a well-established dataset as the baseline ensures proper distribution
of the data across different labels (Figure 1 and Figure 2).

Figure 1: The number of samples for each intent varies, but they are fairly
distributed, with 18 to 68 samples per intent. Figure 2: Slot categories
appearing in at least 30% of the instances are marked as "frequent," while
others are "infrequent." Despite varying slot categories per intent, frequent
ones are evenly distributed.

##  3 Methodology and Experimental Setup

Our experiments were divided into four phases. In our initial experiment, we
employed JointBERT Chen et al. (2019), the state-of-the-art model in this
domain, for both intent detection and slot-filling tasks. In our next
experiment, JointBERT was retained for intent detection, while we explored the
capabilities of GPT-3.5 (Generative Pre-trained Transformer) Brown et al.
(2020) model for slot filling. The third experiment fully utilized GPT-3.5 for
both tasks. For our concluding experiment, we provided GPT-3.5 with the
original intents and then analyzed its performance on the slot-filling task.
The final experiment gives the raw result of slot-filling for the GPT model.

JointBERT leverages the BERT Devlin et al. (2019) model to provide a unified
approach encompassing both intent classification and slot filling by utilizing
the representations from the pre-trained BERT model. We employed the default
BERT tokenizer and maintained consistent parameters for all three languages.
The utilization of these default settings and tokenization methods ensures an
equitable and consistent evaluation across the languages.

GPT-3.5 (Generative Pre-trained Transformer) Brown et al. (2020) model
operates on the Transformer architecture and is adept at generating text
resembling human language by predicting subsequent words or tokens in a
sequence. GPT-3.5’s deep contextual understanding is a result of extensive
pre-training on a diverse corpus of textual data, encompassing various
languages and linguistic intricacies enabling it to excel across a spectrum of
NLP tasks Goyal et al. (2022); Liu et al. (2021); Sakib et al. (2023); Kumar
et al. (2020). We used GPT in a few-shot setting, passing 5 training samples
along with the prompt. Rigorous prompt engineering was performed before
settling on the two prompts for the two tasks. Figure 3 and Figure 4 show the
final versions of the prompts used in the experimentations.

Figure 3: The figure illustrates how the input is formatted for the intent-
detection task. A base-prompt is passed on to the GPT model. A few samples (5)
from the training set are also passed as the context. From these sentence-
output pairs, the LLM understands how the task needs to be solved. Finally,
the current query is passed Figure 4: The input structure for the slot-filling
task is quite similar to the intent detection task. The major difference is
the prompt. For slot-filling, the set of possible slots is based on the intent
type of the query. The intent type is obtained from a separate model and then
from the train set, all possible slots for the given intent are fetched

###  3.1 Experimental Setup

We divided each of the three datasets into training, development, and test
sets using a standard 80-10-10 split. The JointBERT model was trained and
evaluated on an A100 GPU, using a batch size of 8. We closely followed the
setup provided by the original authors for this phase. For GPT, we used the
OPENAI API with the “GPT-3.5-turbo” engine and set the token limit to 50.

##  4 Results

Tables 2 and 3 present the performance of the models we evaluated on our
intent detection and slot-filling tasks. A clear pattern emerges: GPT-3.5
consistently outperforms JointBERT in both tasks.

While intent detection is generally more straightforward, JointBERT performs
reasonably well in this aspect, although it doesn’t quite match the
exceptional performance achieved by GPT-3.5. However, when it comes to the
more intricate task of slot-filling, JointBERT’s performance falls
significantly short, leaving ample room for improvement. In contrast, GPT-3.5
demonstrates its proficiency in handling the complexities of this task.

A significant reason behind GPT-3.5’s superior performance is its broader
exposure to diverse languages during training, including Bangla. JointBERT,
conversely, hasn’t been specifically trained on any Bangla dataset. This
linguistic familiarity gives GPT-3.5 a clear advantage, enabling it to process
and interpret Bangla’s nuances far more effectively than JointBERT. The
results underline the significance of using LLMs for low-resource languages,
especially in scenarios where obtaining high volumes of training data for a
particular downstream task is challenging.

Intent Detection (Accuracy and F1 Score)
---
Models |  | Formal
---
Bangla
| Colloquial
---
Bangla
| Colloquial
---
Sylheti
JointBERT | 0.57 | 0.56 | 0.63 | 0.61 | 0.45 | 0.46
GPT-3.5 | 0.94 | 0.94 | 0.94 | 0.94 | 0.87 | 0.89
Table 2: While the performance of JointBERT is noteworthy for Bangla and its
variants, the GPT-3.5 model excels across all metrics for all three datasets
Slot Filling (F1 Score)
---
| Slot Filling
---
Model
| Intent
---
From
| Formal
---
Bangla
| Colloquial
---
Bangla
| Colloquial
---
Sylheti
JointBERT | JointBERT | 0.14 | 0.11 | 0.07
GPT-3.5 | JointBERT | 0.43 | 0.45 | 0.52
GPT-3.5 | GPT-3.5 | 0.45 | 0.51 | 0.57
GPT-3.5 | Original | 0.54 | 0.53 | 0.57
Table 3:  The slot-filling task is separate from but dependent on the intent
detection task. Intent needs to be passed to the model for good performance.
In slot-filling tasks, GPT massively outperforms JointBERT

##  5 Conclusion

In the era of smart devices, a home assistant’s voice interfaces must resonate
with the authentic linguistic intricacies of its users. Our research presents
the first-ever dataset for intent detection and slot filling in Bangla and
Sylheti, emphasizing their colloquial forms. This focus on colloquial forms
bridges the often-overlooked gap between formal language models and the
nuances of everyday speech. By championing colloquial forms, we ensure a voice
interface that’s more natural and attuned to genuine communication habits.
Through rigorous data collection and validation, we have produced a high-
quality benchmark dataset, providing a solid foundation for subsequent
analyses and model evaluations. The comparative study between large language
models (LLM) like GPT-3.5 and non-LLMs underscores the remarkable capability
of LLMs to excel even with minimal datasets, marking a considerable stride for
underrepresented languages.

##  6 Limitations

While our research has made significant strides in understanding intent
detection and slot filling for Bangla and Sylheti, like any study, it has its
limitations. Our dataset, although carefully curated for the Bangla and
Sylheti variants, is on the smaller side compared to established benchmarks. A
precise and robust data generation process was prioritized, naturally limiting
our data volume. We confined our evaluations to the JointBERT model and
GPT-3.5. The pronounced difference in their performance deterred us from
testing a broader range of models. Moreover, the dearth of optimized Bangla
models for specific tasks posed challenges. An attempt with a Bangla BERT
tokenizer didn’t yield satisfactory outcomes, affecting the JointBERT’s
efficacy. As promising as our results are, they are tied to our specific
dataset and context. Extending our findings to diverse settings or other
languages requires further exploration, marking just the beginning of this
exciting journey.

## References

  * Akbari et al. (2023) Masoud Akbari, Amir Hossein Karimi, Tayyebeh Saeedi, Zeinab Saeidi, Kiana Ghezelbash, Fatemeh Shamsezat, Mohammad Akbari, and Ali Mohades. 2023.  A persian benchmark for joint intent detection and slot filling.
  * Alam et al. (2021) Firoj Alam, Arid Hasan, Tanvirul Alam, Akib Khan, Janntatul Tajrin, Naira Khan, and Shammur Absar Chowdhury. 2021.  A review of bangla natural language processing tasks and the utility of transformer models.  _arXiv preprint arXiv:2107.03844_.
  * Antoun et al. (2020) Wissam Antoun, Fady Baly, and Hazem Hajj. 2020.  Arabert: Transformer-based model for arabic language understanding. arxiv 2020.  _arXiv preprint arXiv:2003.00104_.
  * Bali et al. (2019) Kalika Bali, Monojit Choudhury, Sunaya Sitaram, and Vivek Seshadri. 2019.  Ellora: Enabling low resource languages with technology.  In _Proceedings of the 1st International Conference on Language Technologies for All_ , pages 160–163.
  * Bhattacharjee et al. (2023) Abhik Bhattacharjee, Tahmid Hasan, Wasi Ahmad, and Rifat Shahriyar. 2023.  Banglanlg and banglat5: Benchmarks and resources for evaluating low-resource natural language generation in bangla.  In _Findings of the Association for Computational Linguistics: EACL 2023_ , pages 714–723.
  * Bhattacharjee et al. (2021) Abhik Bhattacharjee, Tahmid Hasan, Wasi Uddin Ahmad, Kazi Samin, Md Saiful Islam, Anindya Iqbal, M Sohel Rahman, and Rifat Shahriyar. 2021.  Banglabert: Language model pretraining and benchmarks for low-resource language understanding evaluation in bangla.  _arXiv preprint arXiv:2101.00204_.
  * Brown et al. (2020) Tom B. Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, Sandhini Agarwal, Ariel Herbert-Voss, Gretchen Krueger, Tom Henighan, Rewon Child, Aditya Ramesh, Daniel M. Ziegler, Jeffrey Wu, Clemens Winter, Christopher Hesse, Mark Chen, Eric Sigler, Mateusz Litwin, Scott Gray, Benjamin Chess, Jack Clark, Christopher Berner, Sam McCandlish, Alec Radford, Ilya Sutskever, and Dario Amodei. 2020.  Language models are few-shot learners.
  * Chen et al. (2019) Qian Chen, Zhu Zhuo, and Wen Wang. 2019.  Bert for joint intent classification and slot filling.
  * Coucke et al. (2018) Alice Coucke, Alaa Saade, Adrien Ball, Théodore Bluche, Alexandre Caulier, David Leroy, Clément Doumouro, Thibault Gisselbrecht, Francesco Caltagirone, Thibaut Lavril, et al. 2018.  Snips voice platform: an embedded spoken language understanding system for private-by-design voice interfaces.  _arXiv preprint arXiv:1805.10190_.
  * Dao et al. (2021) Mai Hoang Dao, Thinh Hung Truong, and Dat Quoc Nguyen. 2021.  Intent detection and slot filling for vietnamese.  _arXiv preprint arXiv:2104.02021_.
  * De bruyn et al. (2022) Maxime De bruyn, Ehsan Lotfi, Jeska Buhmann, and Walter Daelemans. 2022.  Machine translation for multilingual intent detection and slots filling.  In _Proceedings of the Massively Multilingual Natural Language Understanding Workshop (MMNLU-22)_ , pages 69–82, Abu Dhabi, United Arab Emirates (Hybrid). Association for Computational Linguistics.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  Bert: Pre-training of deep bidirectional transformers for language understanding.
  * Ethnologue (2023) Ethnologue. 2023.  Ethnologue 200: Languages of the world.  Accessed on September 3, 2023.
  * Goo et al. (2018) Chih-Wen Goo, Guang Gao, Yun-Kai Hsu, Chih-Li Huo, Tsung-Chieh Chen, Keng-Wei Hsu, and Yun-Nung Chen. 2018.  Slot-gated modeling for joint slot filling and intent prediction.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 2 (Short Papers)_ , pages 753–757.
  * Goyal et al. (2022) Tanya Goyal, Junyi Jessy Li, and Greg Durrett. 2022.  News summarization and evaluation in the era of gpt-3.  _arXiv preprint arXiv:2209.12356_.
  * Hossain et al. (2020) Md Zobaer Hossain, Md Ashraful Rahman, Md Saiful Islam, and Sudipta Kar. 2020.  Banfakenews: A dataset for detecting fake news in bangla.
  * Kumar et al. (2020) Varun Kumar, Ashutosh Choudhary, and Eunah Cho. 2020.  Data augmentation using pre-trained transformer models.  _arXiv preprint arXiv:2003.02245_.
  * Li et al. (2018) Changliang Li, Liang Li, and Ji Qi. 2018.  A self-attentive model with gate mechanism for spoken language understanding.  In _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing_ , pages 3824–3833.
  * Liu and Lane (2016) Bing Liu and Ian Lane. 2016.  Attention-based recurrent neural network models for joint intent detection and slot filling.  _arXiv preprint arXiv:1609.01454_.
  * Liu et al. (2021) Jiachang Liu, Dinghan Shen, Yizhe Zhang, Bill Dolan, Lawrence Carin, and Weizhu Chen. 2021.  What makes good in-context examples for gpt-$3$?  _arXiv preprint arXiv:2101.06804_.
  * Louvan and Magnini (2020) Samuel Louvan and Bernardo Magnini. 2020.  Recent neural methods on slot filling and intent classification for task-oriented dialogue systems: A survey.  _arXiv preprint arXiv:2011.00564_.
  * Mesnil et al. (2013) Grégoire Mesnil, Xiaodong He, Li Deng, and Yoshua Bengio. 2013.  Investigation of recurrent-neural-network architectures and learning methods for spoken language understanding.  In _Interspeech_ , pages 3771–3775.
  * Niu et al. (2019) Peiqing Niu, Zhongfu Chen, Meina Song, et al. 2019.  A novel bi-directional interrelated model for joint intent detection and slot filling.  _arXiv preprint arXiv:1907.00390_.
  * Qin et al. (2019) Libo Qin, Wanxiang Che, Yangming Li, Haoyang Wen, and Ting Liu. 2019.  A stack-propagation framework with token-level intent detection for spoken language understanding.  _arXiv preprint arXiv:1909.02188_.
  * Qin et al. (2021) Libo Qin, Tailu Liu, Wanxiang Che, Bingbing Kang, Sendong Zhao, and Ting Liu. 2021\.  A co-interactive transformer for joint slot filling and intent detection.  In _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 8193–8197. IEEE.
  * Ravuri and Stolcke (2015) Suman Ravuri and Andreas Stolcke. 2015.  Recurrent neural network and lstm models for lexical utterance classification.  In _Sixteenth annual conference of the international speech communication association_.
  * Sakib et al. (2023) Fardin Ahsan Sakib, Saadat Hasan Khan, and AHM Karim. 2023.  Extending the frontier of chatgpt: Code generation and debugging.  _arXiv preprint arXiv:2307.08260_.
  * Schuster et al. (2018) Sebastian Schuster, Sonal Gupta, Rushin Shah, and Mike Lewis. 2018.  Cross-lingual transfer learning for multilingual task oriented dialog.  _arXiv preprint arXiv:1810.13327_.
  * Stoica et al. (2019) Anda Stoica, Tibor Kadar, Camelia Lemnaru, Rodica Potolea, and Mihaela Dînşoreanu. 2019.  The impact of data challenges on intent detection and slot filling for the home assistant scenario.  In _2019 IEEE 15th International Conference on Intelligent Computer Communication and Processing (ICCP)_ , pages 41–47. IEEE.
  * Stoica et al. (2021) Anda Stoica, Tibor Kadar, Camelia Lemnaru, Rodica Potolea, and Mihaela Dînşoreanu. 2021.  Intent detection and slot filling with capsule net architectures for a romanian home assistant.  _Sensors_ , 21(4):1230.
  * Tur et al. (2012) Gokhan Tur, Li Deng, Dilek Hakkani-Tür, and Xiaodong He. 2012.  Towards deeper understanding: Deep convex networks for semantic utterance classification.  In _2012 IEEE international conference on acoustics, speech and signal processing (ICASSP)_ , pages 5045–5048. IEEE.
  * Upadhyay et al. (2018) Shyam Upadhyay, Manaal Faruqui, Gokhan Tür, Hakkani-Tür Dilek, and Larry Heck. 2018.  (almost) zero-shot cross-lingual spoken language understanding.  In _2018 IEEE international conference on acoustics, speech and signal processing (ICASSP)_ , pages 6034–6038. IEEE.
  * Vu et al. (2016) Ngoc Thang Vu, Pankaj Gupta, Heike Adel, and Hinrich Schütze. 2016.  Bi-directional recurrent neural network with ranking loss for spoken language understanding.  In _2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 6060–6064. IEEE.
  * Wang et al. (2020) Jixuan Wang, Kai Wei, Martin Radfar, Weiwei Zhang, and Clement Chung. 2020.  Encoding syntactic knowledge in transformer encoder for intent detection and slot filling.
  * Weld et al. (2022) Henry Weld, Xiaoqi Huang, Siqu Long, Josiah Poon, and Soyeon Caren Han. 2022.  A survey of joint intent detection and slot filling models in natural language understanding.  _ACM Computing Surveys_ , 55(8):1–38.
  * Xu and Sarikaya (2013) Puyang Xu and Ruhi Sarikaya. 2013.  Convolutional neural network based triangular crf for joint intent detection and slot filling.  In _2013 ieee workshop on automatic speech recognition and understanding_ , pages 78–83. IEEE.
  * Zhang et al. (2018) Chenwei Zhang, Yaliang Li, Nan Du, Wei Fan, and Philip S Yu. 2018.  Joint slot filling and intent detection via capsule neural networks.  _arXiv preprint arXiv:1812.09471_.
  * Zhang and Wang (2016) Xiaodong Zhang and Houfeng Wang. 2016.  A joint model of intent determination and slot filling for spoken language understanding.  In _IJCAI_ , volume 16, pages 2993–2999.
  * Zhang et al. (2019) Zhichang Zhang, Zhenwen Zhang, Haoyuan Chen, and Zhiman Zhang. 2019.  A joint learning framework with bert for spoken language understanding.  _Ieee Access_ , 7:168849–168858.

##  Appendix A Appendix

###  A.1 Related Work

Efforts to enhance datasets for intent detection and slot-filling within low-
resource languages, such as Bangla and Sylheti in this context, commence with
the intricate process of translating individual English lexemes extracted from
established benchmarks like ATIS and SNIPS. Previous works in intent detection
and slot filling for low resource languagesDao et al. (2021); Akbari et al.
(2023), have translated each English utterance to their respective languages.
Recent works have shown that there are great performance achievements on
intent detection and slot-filling tasks on datasets that have been derived
from the SNIPS dataset Weld et al. (2022); Qin et al. (2019); Wang et al.
(2020), and this gives a reason to choose the SNIPS dataset over the ATIS
dataset as it is a good starting point for a work with a language that has
never been explored.

Spoken Language Understanding, a pivotal endeavor in the domain of task-
oriented dialogue systems, encompasses the tasks of intent detection and slot-
filling. Traditionally, these tasks were regarded as distinct domains in which
significant progress was made Tur et al. (2012); Ravuri and Stolcke (2015);
Mesnil et al. (2013); Vu et al. (2016). However, recent research has garnered
notable attention by achieving remarkable advancements in performance through
the concurrent learning of intent detection and slot-filling tasks Zhang et
al. (2018); Weld et al. (2022). In this section, we’re primarily looking at
how intent detection and slot-filling tasks are combined. We’ll focus on two
well-known strategies for this integration:

  * •

A strategy devised through parameter sharing and the exchange of hidden
states, utilizing a common BiLSTM/BERT encoder, along with two distinct
decoders dedicated to intent detection and slot filling, on top of the shared
encoder. Chen et al. (2019); Xu and Sarikaya (2013); Liu and Lane (2016);
Zhang and Wang (2016).

  * •

Another strategy, extending the initial approach to a more advanced level,
involves the model acquiring an understanding of the relationships between
slots and intent labels. This frontier has been explored in research in two
distinct ways. Some studies Goo et al. (2018); Li et al. (2018); Niu et al.
(2019) have demonstrated the use of attention mechanisms to discern the
correlation between the overarching intent context representation and the slot
vectors generated by the encoder. Alternatively, other works Qin et al.
(2019); Zhang et al. (2019) have approached this by initially learning the
representation of the utterance, which aligns with the representation of the
global intent context, utilizing a self-attention mechanism. Subsequently,
they join this representation with the encoder’s vector outputs before feeding
the combined vectors into the slot-filling decoder.

###  A.2 Examples from the dataset

Here we include a few examples from each of the datasets.

Figure 5: Few examples from the Formal Bangla dataset. (Input sentence -
Intent - Expected slots) Figure 6: Few examples from the Colloquial Bangla
dataset. (Input sentence - Intent - Expected slots) Figure 7: Few examples
from the Sylheti dataset. (Input sentence - Intent - Expected slots)

###  A.3 Prompts used for GPT

For the intent detection task we used the following prompt: "You are a
language model for classifying the intent of the given text. There are 10
intent classes. These are: BookRestaurant, ShareETA, GetPlaceDetails,
ShareCurrentLocation, ComparePlaces, GetDirections, GetTrafficInformation,
RequestRide, SearchPlace, GetWeather. You are given a text in Bangla and you
have to classify it into one of these classes. Only return the class name."

In this approach, we clearly outlined the potential intent classes, specified
the input language as Bangla, and directed the model to solely return the
class name. Such structuring was essential to elicit precise responses from
the model.

For our slot-filling task, we utilized the following prompt: "You are a
language model for slot filling. You will be given a Bangla sentence. You need
to classify each word in the given sentence according to the appropriate slot.
The possible slots are as follows: list of possible slots extracted from the
train set (based on the training intent)"

We equipped the model with both the potential slots and their associated
intent. Notably, the performance fluctuated depending on the source of the
intent— GPT-3.5, JointBERT, or the Original dataset.

###  A.4 Inter-annotator metrics

In order to assess inter-annotator agreement, this study utilized two primary
evaluation metrics: Cohen’s Kappa and Average BLEU.

Cohen’s Kappa provides a statistical measure of agreement between two
annotators, while accounting for the possibility of chance agreement.
Specifically, it involves calculating the actual observed agreement between
the annotators and comparing that to the level of agreement that would be
expected by random chance. Cohen’s Kappa expresses the ratio between these two
values as a score ranging from 0 to 1, with higher scores indicating greater
reliability.

Average BLEU (Bilingual Evaluation Understudy) is a commonly employed metric
for evaluating machine translation outputs by comparing them against one or
more reference translations. It analyzes the co-occurrence of n-grams between
the translated text and human reference texts to produce a score reflecting
the quality and fluency of the translation. Taking the average BLEU score
across multiple translations provides an overall indicator of the fidelity of
the translations with respect to the reference materials.

Together, these two metrics enable analysis of both the reliability of
individual annotators via Cohen’s Kappa and the accuracy and fluency of
translations via Average BLEU in relation to trusted references. The
combination provides a robust means of evaluating key aspects of annotation
quality for this study.


Original Article Journal Section Philipp Berens, Institute for Ophthalmic
Research, University of Tübingen, Tübingen, Germany philipp.berens@uni-
tuebingen.de The German Research Foundation: BE5601/4-1 and BE5601/8-1; the
Collaborative Research Center 1233 “Robust Vision”: 276693517; Individual
research grants: EU 42/10-1 and BE5601/6-1; the German Ministry of Education
and Research: 01GQ1601. ALD, automatic locality determination; ASD, automatic
smoothness determination; LG, Linear-Gaussian model; LNLN, Linear-Nonlinear
cascade model; LNP, Linear-Nonlinear-Poisson model; MAP, maximum a posteriori;
MSE, mean squared error; MLE, maximum-likelihood estimate; NMF, nonnegative
matrix factorization; RGC, retinal ganglion cell; STA, spike-triggered
average; STRF, spatio-temporal receptive field; wSTA, whitened-STA.

# Estimating smooth and sparse neural receptive fields with a flexible spline
basis

Ziwei Huang  Institute for Ophthalmic Research, University of Tübingen,
Tübingen, Germany  Yanli Ran  Institute for Ophthalmic Research, University of
Tübingen, Tübingen, Germany  Centre for Integrative Neuroscience, University
of Tübingen, Tübingen, Germany  Jonathan Oesterle  Institute for Ophthalmic
Research, University of Tübingen, Tübingen, Germany  Thomas Euler  Institute
for Ophthalmic Research, University of Tübingen, Tübingen, Germany  Centre for
Integrative Neuroscience, University of Tübingen, Tübingen, Germany  Philipp
Berens

###### Abstract

Spatio-temporal receptive field (STRF) models are frequently used to
approximate the computation implemented by a sensory neuron. Typically, such
STRFs are assumed to be smooth and sparse. Current state-of-the-art approaches
for estimating STRFs based on empirical Bayes are often not computationally
efficient in high-dimensional settings, as encountered in sensory
neuroscience. Here we pursued an alternative approach and encode prior
knowledge for estimation of STRFs by choosing a set of basis functions with
the desired properties: natural cubic splines. Our method is computationally
efficient and can be easily applied to a wide range of existing models. We
compared the performance of spline-based methods to non-spline ones on
simulated and experimental data, showing that spline-based methods
consistently outperform the non-spline versions.

###### keywords:

Spatio-temporal receptive field estimation, natural cubic spline

##  1 Introduction

Spatio-temporal receptive fields (STRFs) are frequently used in neuroscience
to approximate the computation implemented by a sensory neuron. Such models
consist typically of one or more linear filters summing up sensory inputs
across time and space, followed by a static nonlinearity and a probabilistic
output process [1]. For this, different distributions can be used depending on
the type of recorded neuronal responses.

The simplest way to estimate a STRF for a given data set is arguably to
compute the spike-triggered average (STA), the average over all stimuli
preceding a spike [2]. However, the results of the spike-triggered analysis
are often noisy due to limited data and prone to be distorted when the neuron
is stimulated with correlated noise or natural stimuli [3]. The maximum
likelihood estimate (MLE), also known as whitened-STA (wSTA) [4], does not
suffer from such shortcomings, but it requires even more data to converge as
the whitening procedure tends to amplify noise at frequencies with low power,
making it a less preferred option in many experimental settings. A modified
wSTA that removes low power frequencies can result in smoother STRF, while the
threshold for such a frequency cutoff is normally chosen by cross-validation
[4].

It has been suggested that these shortcomings can be overcome by computing the
maximum a posteriori (MAP) estimates within the framework of Bayesian
inference [5]. Under this framework, prior knowledge about STRFs such as
sparsity, smoothness [6], and locality[7] can be encoded into a prior
covariance matrix, whose hyperparameters are learnt automatically from the
data via evidence optimization (also known as empirical Bayes, or Type II
maximum likelihood in the literature[8]). These algorithms thus provide STRF
estimates with the desired properties even with little or noisy data. Although
elegant in theory, the evidence optimization step of this framework is
computationally costly: the cubic complexity in the dimensionality of the
parameter space renders it difficult to scale to high-dimensional settings as
often encountered in sensory neuroscience [5]. Also, current popular
implementations of MAP estimates work mostly under the Linear-Gaussian
encoding model, which assumes additive Gaussian noise to the neural response,
and are restricted to simple one filter models.

An alternative approach to encoding prior knowledge for estimating STRFs is to
parameterize the STRF with a set of basis functions with the desired
properties. In this way, the model fitting becomes much more efficient as the
number of parameters to be optimized is significantly reduced, and some degree
of smoothness is automatically enforced. Previous studies used a raised-cosine
basis (Figure 1A) for 1D retinal ganglion cell STRFs in macaque monkey [9] and
Gaussian spectral or Morlet wavelet basis for 2D auditory STRFs in ferrets
[10]. Those basis functions are chosen to capture certain STRF properties in
their corresponding studies, thus they are relatively inflexible to generalize
to other settings. Besides, they are typically controlled by multiple
hyperparameters, which creates a burdensome and time-consuming model selection
process.

Figure 1: Illustration of the raised cosine and natural cubic spline basis. A:
A typical raised-cosine bases as in [9] controlled by several parameters, such
as number of basis functions, the degree of nonlinearity and the placement of
endpoints. The choice of these parameters determines the final goodness of
fit. B: An equally spaced 1D natural cubic spline, controlled by only one
parameter: the number of basis functions (aka the degree of freedom). C:
Fitted 1D filters using different bases and different sets of parameters. D:
An example of 2D natural cubic splines. Left: Two selected bases of 2D natural
cubic splines. Middle: Weighted bases after fitted to a simulated 2D STRF.
Right: An example of a fitted 2D STRF overlaid on the ground truth.

Here, we propose to use an alternative basis for receptive field estimation:
natural cubic splines (Figure 1B). These are known as the smoothest possible
interpolant and can be easily extended to high-dimensional settings with
minimal assumptions and few hyperparameters [11]. This spline basis can be
incorporated into a wide range of existing receptive field models (Figure 2),
such as the Linear-Gaussian (LG) and the Linear-Nonlinear-Poisson (LNP)
model[12, 9] for single STRF estimation as well as Linear-Nonlinear-Linear-
Nonlinear (LNLN) cascade models[13, 14] and spike-triggered clustering methods
[14, 15] for subunit estimation. When combined with proper regularization,
smooth and sparse STRFs can be retrieved with very little data and in a
computationally efficient manner, while their predictive performance reaches
state-of-the-art level. Furthermore, we introduce diagnostic tools for
assessing the quality of an estimated STRF: (1) a significance test based on
the analytical solution of the STRF covariance matrix and (2) a permutation
test based on model prediction. We provide all developed methods as part of
the RFEst Python toolbox.

##  2 Results

We developed the RFEst Python toolbox for spline-based STRF estimation in
various models, including the single-filter LG and LNP models as well as the
multi-filter LNLN cascade model (Figure 2). We measured the performance of the
spline-based methods and compared it with their non-spline versions or other
previous methods, using simulated data as well as neural recordings from the
retina and the visual cortex. To quantify the goodness of fits for the
simulated benchmark data, we measured the mean squared error (MSE) between the
ground truth STRF and the estimated STRFs for various methods and studied
their dependence on the amount of available data. For the experimental data,
we computed the correlation between the predicted and measured responses on a
held-out test set as a performance measure.

Figure 2: Overview of the different spatio-temporal receptive field (STRF)
models used in this study. A: Linear-Gaussian (LG) encoding model.
$y=Xw+\epsilon,\epsilon\sim N(0,\sigma^{2})$, where $X$ is the stimulus design
matrix, $y$ the response, $w$ the STRF. B: Linear-Nonlinear-Poisson (LNP)
model. $y=f(Xw),y\sim Poisson$, where $f(x)$ is the filter nonlinearity. C:
Linear-Nonlinear-Linear-Nonlinear-Poisson (LNLN) model. $y=g(\sum f(Xw)),y\sim
Poisson$, where $g$ is the output nonlinearity. Optionally, an autoregressive
term convolving previous response with a response-history filter can be added
before the nonlinearity to improve the model predictive performance.

###  2.1 Estimating receptive fields from simulated LG, LNP and LNLN models

First, we simulated the responses of a sensory neuron $y$ with additive
Gaussian noise using a STRF with 30 time frames and 40 pixels in space in
response to white and pink noise stimuli for various stimulus lengths (see
Methods 4.1 for details). First, we fit LG models (Figure 2A) to these
responses using various STRF estimation techniques.

We measured the MSEs between the ground truth STRF and the different STRF
estimates with respect to the ratio of the number of samples $n$ to the number
of parameters $d=30\times 40$, averaged over 10 different random seeds (Figure
3A). As expected, the STA and wSTA approached the ground truth only in high
data-to-parameters regimes. For limited data (e.g. n/d=4), the estimated STRF
showed substantial background noise (Figure 3B). In the case of pink noise,
the STA was in fact never able to recover the STRF due to the high correlation
between pixels in the stimulus [3] (Figure 3C-D).

We also applied more advanced techniques to estimate the STRF to the simulated
data, such as automatic smoothness determination (ASD, aka MAP with smoothness
prior [6]) and automatic locality determination (ALD, aka MAP with locality
prior [7]). These estimates were much closer to the ground truth already with
a relatively small amount of data, with ALD outperforming ASD for both white
and pink noise (Figure 3). However, their computation time for estimating the
optimized prior was long due to the high computational complexity.

We then applied a spline-based version of the wSTA (SPL wSTA) with and without
L1 regularization on the basis function coefficients (see Methods). This
simple change of basis functions led to a STRF estimator which performed very
well for all n-d-ratios tested, even better than ASD and ALD, especially when
using L1 regularization (Figure 3). At the same time, this estimator was also
computationally very efficient, as the number of coefficients to be estimated
was reduced to the number of basis functions $b$ (for this example, the STRF
had $d=30\times 40=1200$ parameters, whereas the number of basis functions was
only $b=9\times 12=108$)

Figure 3: Linear-Gaussian (LG) Model. 2D simulated RF stimulated by white and
pink noise in various lengths. A and C: Similarity between ground truth and
estimated RFs measured as mean squared error (MSE) as a function of the amount
of simulated data. Error bars are computed by averaging over 10 trials with
different random seeds. B and D: Ground truth and examples of estimated RFs by
different methods at n/d=4, indicated by vertical grey dashed line in A and C.

We next simulated responses $y$ from an LNP model (Figure 2B) with Poisson
noise similar to before using white and pink noise for stimulation. For
simplicity, we used an exponential nonlinearity and treated it as known during
the estimation procedure. We estimated the linear filter from different
amounts of data using the STA, wSTA, ASD and ALD, as well as spline and non-
spline LNP fitted via gradient descent with L1 regularization (referred to as
LNP+SPL and LNP, respectively). Note that the original ASD and ALD were not
designed for the LNP but the LG model (though some progress has been made for
the former [16]). Despite this model mismatch, they still recovered good
estimates of the STRF (as has been demonstrated in the original papers [6,
7]).

Again, we found that different methods resulted in estimates of different
quality depending on the available amount of data (Figure 4). Here, in
contrast to the LG model, we changed the x-axis of Figure 4A and C to time in
minutes for better presentation, as the model performance was dependent on the
number of effective samples (frames that contain at least one spike) instead
of the total number of samples, and the stochastic nature of the Poisson spike
generator yielded different number of spikes for each trial. For both stimulus
conditions, we found that the spline-based LNP yielded extremely data-
efficient estimates of the linear filter compared to all other methods. In
contrast, gradient-based optimization of an LNP model without spline basis,
even with L1 regularization, performed only slightly better than the STA and
wSTA computed under the mismatched LG model assumptions. The gap in
performance of the LNP+SPL and MAP estimators of ASD and ALD were larger here
compared to the LG case. This was partially due to the fact that the number of
effective samples was smaller in Poisson spiking data, but the mentioned model
mismatch between the cost function of the LG model optimized by the ASD and
ALD methods and the data likely also contributed.

Figure 4: Linear-Nonlinear-Poisson (LNP) Model. A simulated 2D RF stimulated
by white and pink noise for various stimulus lengths. A and C: Details are the
same as in Figure 3. B and D: Example estimated RFs by different methods at
data=4 minutes, indicated by vertical grey dash line in A and C.

We finally investigated how recently proposed subunit models (Figure 2C), in
particular an LNLN cascade model [13, 17] and spike-triggered clustering or
non-negative matrix factorization methods[15, 14], can be made more data-
efficient using a spline basis. We compared their performance with and without
spline augmentation (see Methods 4.6). The assumption behind the spike-
triggered clustering methods is that each spike can be attributed to a
specific subunit of a STRF, so that each subunit can be viewed as a STA of a
cluster of spikes. If this assumption holds, then a spline-approximated STA
should be able to serve as a better centroid for clustering, as the smoothed
STA is closer to the ground truth subunit STRF. We used k-means clustering as
an in-place of the previous published soft-clustering method [15] and used a
multiplicative update rule[18] for semi-nonnegative matrix factorization
(semi-NMF) [14]. Moreover, we modified the previous published semi-NMF [14]
and imposed the nonnegative constraint to the weight matrix instead of the
subunit matrix, which allowed us to retrieve antagonistic STRFs.

Figure 5: Subunit models. Two antagonistic difference-of-Gaussian 2D RFs
stimulated by white noise. The subunit STRFs were estimated by k-means
clustering, semi-NMF and LNLN, without and with splines. Examples shown are
estimated subunits by different methods at data = 4 and data = 32 minutes,
indicated by vertical light and dark grey dash line in the bottom panel,
respectively. Similarity between ground truth and estimated RFs measured as
mean squared error (MSE) as a function of the amount of simulated data.

We simulated responses $y$ with Poisson noise from the model by stimulating
two 20x20 pixels antagonistic center-surround subunits with different lengths
of white noise stimuli. Both filter outputs first went through an exponential
nonlinearity, were then summed up, and finally went through a softplus output
nonlinearity.

We estimated subunit STRFs using k-means clustering or semi-NMF of the spike-
triggered ensemble, and maximum-likelihood estimation of the LNLN model
directly, with or without a spline basis (Figure 5). Similar to the results
from previous simulations, retrieved subunit STRFs from spline-based models
were generally much more similar to the ground truth than models without
spline, and required a smaller amount of data to achieve a similar level of
performance. For example, when using a spline basis, we were able to retrieve
the two subunit filters at high quality using only 4 minutes of data. With
this amount of data, all other methods still showed essentially only noise
with a hint of signal in the estimate of the linear filter.

It is worth noting that the computational efficiency of LNLN estimation was
improved a lot by using a spline basis, as the number of coefficients to be
estimated was significantly reduced. In contrast, the computation time for
both spike-triggered clustering methods was actually increased by
incorporating splines, because more computational steps were added to the
original algorithms in those modified versions (see Methods 4.6).

###  2.2 Comparing computation time between different methods and model
selection

It was not straightforward to fairly compare the computation time between MAP-
derived and spline-based GLMs, as the most time-consuming part for MAP was to
estimate the prior hyperparameters via evidence optimization, while spline-
based GLMs need to resort to cross-validation for selecting the number of
basis functions ($b$, also known as the degrees of freedom) for each
dimension. Yet, we can still get an impression of the algorithmic efficiency
of each method by profiling their most time-consuming step.

Therefore, we compared the time taken by the computation of the STA (see
equation 1 in Methods), wSTA (equation 2), MAP (equation 3) and SPL wSTA
(equation 12) for different amounts of data by separately profiling the single
most time consuming computation (Figure 6). Not surprisingly, STA was the
fastest, as it had a linear time complexity ($O(n)$) and its computation time
depended mainly on the amount of data ($n$), while wSTA and MAP were the
slowest, as they were both limited by the calculation of the inverse matrix,
and their computation time grew in a cubic manner with respect to the number
of coefficients ($w$) in the STRF due to the cubic time complexity
($O(w^{3})$) of the matrix inversion. SPL on the other hand, was less time-
consuming than wSTA and MAP, even though the time complexity was still cubic.
However, the improvement was due to the reduced number of model parameters by
parameterizing the STRF coefficients with the spline basis ($b$ coefficients),
which resulted in a much smaller covariance matrix to invert ($O(b^{3})$),
thereby speeding up the most time-consuming step of the algorithm
substantially. Thus, spline-based STRF estimators offer the best available
compromise between STRF estimation quality and computation time.

Figure 6: Comparison of computation time between different methods with
different number of STRF coefficients, amount of data and number of spline
basis. MAP refers to both ASD and ALD methods, as the size of the covariance
matrix was the same for both techniques. Measurements of computation time were
averaged over 10 repetitions and performed on a 2019 16-inch MacBook Pro with
a 2.3 GHz 8-core Intel Core i9 and 16 GB of RAM. See Method 4.10 for more
information.

The total amount of time spent on finding an optimal STRF depends on many
factors, such as how the model is initialized, the choice of learning rate and
stopping criteria for gradient descent optimization. Importantly, the way
hyperparameters are selected also contributes significantly: most commonly,
the model prediction performance is monitored for on a held-out data set over
a discrete grid of hyperparameters. Even for evidence optimization methods, in
which the main selling point is to learn the optimal hyperparameters
automatically, cross-validation can still be used to assess model goodness of
fit and generalization performance [6, 19]. For spline-based GLM, this
requires refitting the same model as often as the number of hyperparameter
sets. Given that splines are equally spaced, a good strategy to accelerate the
process would be to measure the performance of SPL wSTA over a grid of
possible candidates from very few to many, choose the number when the
performance is starting to plateau, then fit SPL+L1 with gradient descent to
retrieve a sparse STRF.

###  2.3 Diagnostic tools

Figure 7: Examples of diagnosing STRF quality with confidence interval,
p-value, prediction performance. A: Ground truth STRF with the dimensions
25x25x25. The 3D STRF is the Kronecker product of the respective temporal and
spatial RF. B: A well fitted STRF with the sample size equals to $d$. (i) the
estimated spatial RF at the minimum of temporal RF, indicated by the gray
vertical dotted line in (ii); (ii) the estimated temporal RF; (iii) the change
of correlation between the data and predicted response over all gradient
descent iteration. (iv) the predicted performance on the permuted validation
set. C: Same as B but with half amount of data. D: Same as B but with stronger
additive Gaussian noise. E: Same as B but the simulated response were
permuted.

In addition, even if a well-regularized model is chosen according to
predictive performance, additional diagnostic tools can be useful to decide
whether a certain STRF is reasonable and shows a clear structure (Figure 7).
We devised three diagnostic tools that can be employed to this end: (1) the
confidence interval (CI) of the fitted STRF and its corresponding p-values
based on Wald statistics (with the null hypothesis that all STRF coefficients
are zero, see also Method 4.7); (2) the difference between the metrics of the
training and the validation set over all iterations of the optimization; and
(3) a permutation test on model prediction for validation and/or test data set
(Method 4.8).

To demonstrate their usefulness, we investigated a ground truth 3D STRF under
a Linear-Gaussian model. We showcased four different scenarios (Figure 7):
First, we looked at a well-fitted STRF (Figure 7A). Here, the CI was narrow
and significantly different from an all-zero STRF (p<0.001), the gap between
the training and validation prediction metrics was small and the prediction on
the validation set was better than chance (p<0.001). The confidence interval
and its corresponding p-value were determined by the training data. Next, we
looked at scenarios where the amount of data was small or the signal-to-noise
ratio (SNR) was low (Figure 7 C and D). Here, the STRF did not pass the
significance test due to a lack of power, but it is still possible that the
model predictions on the validation set were better than chance. Additionally,
if the gap between the performance metric on the training and validation sets
gets too large, it would be a good idea to refit the model with a stronger
regularization. Finally, in the case of the data with extremely low SNR
(Figure 7E), which we generated by randomly permuting the simulated response,
the prediction performance of the fitted STRF was just at the chance level
regarding all indicators.

###  2.4 Application to experimental data

We next documented the complete process of using gradient-based estimation for
the parameters of spline-based LNP and LNLN models in a previously published
experimental data set: extracellularly recorded spikes from tiger salamander
retinal ganglion cells (RGC), stimulated by 3D non-repeat checkerboard white
noise [14] (Figure 9A). In addition, we applied the same techniques to
2-photon calcium imaging data recorded from a mouse retinal ganglion cell soma
stimulated by checkerboard white noise (Figure 9B) (using the positive
gradient of the calcium trace as the response). We compared the results of
spline-based models with non-spline versions.

For both data sets, only the first 10 minutes of the whole recording were used
for fitting the models, another 2 minutes were used for model validation and
another 4 sets of 2 minutes for testing. The predictive performance was
measured by the average of the four correlation coefficients of the predicted
and measured responses. We set the number of subunits $k=4$ throughout. The
dimensionality of the 3D STRFs were [20, 25, 25] and [25, 15, 15],
respectively. We added diagnostic information as in Figure 7. Moreover, a
contour was drawn on the spatial filter to illustrate the spatial extent of
the STRF.

We fitted the non-spline LNP and LNLN with different L1 regularization weights
from 0 and 10, with steps of 1.

For spline LNP and LNLN, the number of basis functions (aka degree of freedom,
df) for temporal and spatial dimensions (assuming x and y dimensions here used
the same df) were first selected based on the predicted performance of the
maximum likelihood estimates on the validation set (from 1/3 to 2/3 of the
pixels in each dimension, with steps of 1). Finally, for data set (1), the
grid search took 11.2 minutes on the same hardware setup used for Figure 6 and
yielded the optimal df=[9, 9, 9] ; for data set (2), it took 2.4 minutes and
yielded the optimal df=[12, 9, 9] (Figure 8). See also how different amounts
of training data affects the optimal df in Supplementary Section 1.

Figure 8: Choosing the number of basis functions with cross-validation. The
predictive performance of the SPL wSTA on the validation set for a grid of
different degrees of freedom (df). The optimal df are indicated by red
crosses. A: Spike recordings from a salamander RGC. B: 2-photon calcium
imaging recordings from a mouse RGC (see Methods).

Then, we fitted spline LNP and LNLN with L1 regularization weights ranged from
0.5 and 1.5, with steps of 0.1 (as we observed, the spline model required a
smaller value for regularization due to the reduction of model parameters).
The ordered grid search for the L1 weight, starting with the smallest value,
was interrupted if the current predictive performance was lower than the
previous one. All models were initialized with maximum likelihood estimates
(plus small additive Gaussian noise) and fitted by gradient descent with 1500
iterations, but an early stop was also deployed with two triggers: either the
training cost changed less than 1e-5 for the last 10 iterations, or the
validation cost monotonically increased for the last 10 iterations. Sometimes
the optimization didn’t stop after reaching the lowest validation cost due to
the adaptive learning rate resulting in the cost changing in a zigzag manner.
But regardless of the stopping trigger, the model at the lowest validation
cost was returned. The computation time of each fit for all four models was
summarized by box plots in Figure 9Aiv and Biv.

Fitting a 3D STRF with gradient descent directly was considered infeasible
unless certain regularization was imposed [17], and the situation would become
even worse for subunit models, as the number of model parameters increases
linearly with the number of subunits. Previously methods [14, 15] (as shown in
Figure 5) relied on dimensionality reduction techniques to retrieve subunit
STRFs but those methods could not be applied to 2-photon calcium imaging data
directly without an extra step of spike inference. Spline GLMs not only
overcome both of these problems but also allow us to assess the quality of the
estimated STRF with CI and their corresponding p-values (In principle, the
same statistics can also be applied to non-spline GLMs, yet a larger number of
model coefficients would result in high computational cost for computing the
STRF posterior covariance).

We found that for both example neurons from the two data sets, STRFs estimated
by non-spline models were noisy and prediction performance was only moderate
(Figure 9Ai, Bi, Aiii, Biii). This was especially prominent for LNLN
estimates, as here the number of parameters to estimate exceeded the amount of
available data to fit the models (Figure 9Aii and Bii). In contrast, the
spline-based LNP and LNLN estimated in general were smooth and sparse even
with a severely limited amount of data, and they consistently outperformed
their non-spline versions (Figure 9A and B).

The number of subunits that could be retrieved by the model was limited by the
amount of available data (total length of the recording, and for spike data,
the total number of spikes) as well as the signal-to-noise ratio of the data.
With the limited data in our examples, we found three and two non-overlapping
subunits for the two neurons, respectively (Figure 9Aii and Bii). Additional
subunits available to the model were pushed to near zero by the L1
regularization and the significance test failed to reject the null hypothesis
that all coefficients in those subunit STRFs were close to zero. But as noted
in the previous section, the p-value should not be the only judge for the
quality of the estimated STRF. As shown in 9B, the STRF from the spline LNP
model was also non-significant but looked reasonable from the snapshot
visualization. In cases like this, other tools such as the permutation test
shown in the previous section would be required to determine whether the STRF
in question was good or not.

We provided two more examples of visual neurons (RGC and V1) and one example
of non-visual neuron (entorhinal grid cell) in Supplementary Section 2.

Figure 9: Spline-based LNP and LNLN consistently outperformed their non-spline
versions on both example neurons from two data sets: A: Spike recordings from
a salamander RGC, B: 2-photon calcium imaging data from a mouse RGC. (see
Methods 4.1 Experimental data for details). i and ii: STRFs estimated from LNP
and LNLN with or without spline, where number of subunits $k=4$ for all LNLN
models; Shading for temporal filters show 95% confidence interval (see
Methods, 4.7), which is too small in A to be visible. iii: predictive
performance comparison of LNP and LNLN with and without spline, measured by
the averaged correlation coefficients of the test set response and the
predictive ones. Error bars were computed from 4 sets of test data. iv:
Computation time of each fit for the same models were summarized in the box
plots.

##  3 Discussion

Here, we showed that natural cubic splines form a good basis for STRF
estimation by tackling one of the main challenges in sensory neuroscience: to
efficiently estimate high-dimensional STRFs with limited data. Using this
basis function set, we showed that inference for single filter LN models with
different noise models as well as LNLN cascade models can be efficiently
performed from both spiking as well as two-photon imaging data, even if the
signal-to-noise ratio is low (and see also how spline-based methods compare to
previous methods in Table 1). We provide an implementation of the proposed
methods in the RFEst-package for Python.

Methods | Stimuli Type | Noise Distribution | Data needed | Complexity | Multiple Filters
---|---|---|---|---|---
STA | White Noise only | Gaussian | a lot | $O(n)$ | No
wSTA | Any | Gaussian | a lot | $O(w^{3})$ | No
MAP | Any | Gaussian | little | $O(w^{3})$ | No
LNLN | Any | Poisson | a lot | $O((wk)^{3})$ | Yes
SPL LG | Any | Gaussian | little | $O(b^{3})$ | No
SPL LNP | Any | Poisson | little | $O(b^{3})$ | No
SPL LNLN | Any | Poisson | little | $O((bk)^{3})$ | Yes
Table 1: Model comparison. Here, $n$ is the number of samples, $w$ the number
of STRF coefficients, $b$ the number of spline basis coefficients and $k$ the
number of subunits. In most cases, $b\ll w$.

###  3.1 Choice of basis functions

We chose the natural cubic spline basis to represent the linear filters as
default in our toolbox, because among many other spline bases such as cyclic
cubic splines, B-splines, or thin-plate splines, natural cubic splines have
been shown to be the smoothest possible interpolant through any set of data
and theoretically near-optimal to approximate any true underlying functions
closely [11]. Conveniently, assuming all knots are spaced equally in each
dimension, natural cubic splines only require setting one hyperparameter: the
degree of freedom, or the number of basis functions. This is unlike the
raised-cosine basis used previously[9], which is governed by multiple
hyperparameters. Also, natural cubic splines can be easily extended to high-
dimensional settings, also known as tensor product bases, by simply taking the
row-wise Kronecker product of the basis matrix along each dimension (see
Method 4.4), making them uniquely suitable for estimating STRFs. There is also
an interesting connection between smoothing splines and Gaussian processes. As
shown by [20, 21, 21], spline regression converges to the MAP solution of a
Gaussian process regression under certain constraints.

Through extensive validation with simulated and real experimental data sets,
we showed that STRFs estimated with this basis can not only be more accurate
than previous state-of-the-art methods, but the estimation procedure was also
computationally much more efficient. The estimated STRFs were smoothed
automatically, as they were the linear combination of a set of smooth basis
functions. Further, sparsity and localization could often be achieved by
adding an L1 regularization term to the loss function, which pushed the
coefficient for the less relevant basis functions to zeros. Moreover, by
parameterizing STRFs with basis functions, the number of coefficients could be
significantly reduced, hence model fitting became much more computationally
efficient, and the amount of data needed for the optimization to converge was
also significantly less.

###  3.2 Potential pitfalls of using splines in non-ideal situations

One caveat of using splines to fit STRFs is that some STRFs might not need
smoothing, thus using splines would lead to suboptimal results. For example,
if a STRF occupies only one pixel in space or the amplitude changes sharply in
the temporal profile, spline-based model would likely yield suboptimal results
by yielding overly smoothed filters. However, one can identify these cases by
monitoring the cross-validation performance on held-out data while increasing
the number of basis functions in each dimension. If the performance keeps
increasing as the number of basis functions reaches the number of pixels in
each dimension, a spline-free model could be preferable.

Another concern of using splines is that it would be difficult to tell a poor
fit from a good one in some cases, as a poorly-fitted spline-based STRF would
also look like a smooth STRF, unlike pixel-based methods that yield a noisy
STRF when the fit is poor. In these cases, the diagnostic tools we introduced
in Section 2.3 and Figure 7 would be helpful to assess the quality of the
STRFs.

###  3.3 Other approaches to efficient estimation of receptive fields

Similar to our approach, fast automatic smoothness determination (fASD) [22]
has been suggested to speed up the process of evidence optimization in ASD.
Because the smoothness prior covariance can be represented in frequency
coordinates and many high-frequency coefficients of the Fourier-transformed
filter tend to be small, the size of the prior covariance matrix can be
reduced by cutting off those small value frequencies (an idea similar to the
modified wSTA [4]), leading to a speed-up in the optimization process.
However, as the efficiency of fASD is controlled completely by the initial
frequency truncation, and the level of frequency to be truncated has to be
fixed before the evidence optimization starts, one needs to resort back to
grid search and cross-validation for choosing a good starting point, making
the automatic hyperparameters selection via evidence optimization less
automatic and comparable to the need to find the optimal degrees of freedom
for the spline basis. Also, as fASD (and original ASD) rely on only a few
hyperparameters, it might not even require to do evidence optimization at all
but just to find the optimal hyperparameters through cross-validation
completely. Importantly, fASD is so far only applicable within the framework
of linear models with Gaussian additive noise, while our approach is ready to
use in nonlinear models with non-Gaussian noise, as shown for LNP and LNLN
models.

Regarding the inference for multiple STRFs in LNLN-cascade models, previous
work has attempted to provide efficient algorithms by first retrieving subunit
STRFs through spike-triggered clustering methods such as soft-clustering [15]
or semi-NMF [14]. These were subsequently used to estimate other model
parameters via gradient descent while holding the subunit STRFs fixed.
Utilizing only the collection of stimuli that elicit a spike, these approaches
have advantages in terms of the computational cost in time and memory compared
to the direct fitting of the complete model using gradient descent [13, 17].

The ease of using these algorithms relying on the spike-triggered ensemble
under white noise stimulation is also the disadvantage of these approaches –
they can only be applied to spike data under white noise stimulation. In many
cases, experiments today record more diverse measurements from neurons, such
as two-photon calcium imaging [23, 24] or synaptic current recording [25]
under more diverse types of stimulation, such as correlated noise and natural
stimuli [26]. In those cases, spike-triggered analysis can not be used
directly without modification, as the spike-triggered average is no longer the
maximum likelihood estimate [3]. In contrast, our method does not rely on
Gaussianity of the stimulus or the discreteness of the response, and in
principle can be generalized to other stimuli distribution and other response
types.

In addition, current spike-triggered clustering methods rely on temporally
collapsing the 3D stimuli to further reduce the model complexity, assuming all
subunits share the same temporal response kernel. This limits the usability of
these methods to some very specific cell types and will not account for the
spatio-temporal inseparable computation implemented by neurons receiving
inputs from multiple sources, such as bi-stratified On-Off direction-selective
cell in the retina [27]. In this case, a more general approach such as our
method should be preferred.

###  3.4 Further extensions

Here, we only showcased how to estimate STRFs using a natural cubic spline in
a set of popular models with standard gradient descent with L1 regularization.
This simple spline add-on can be directly incorporated into other previous
models [13, 17] with different constraints, regularization and optimization
techniques, as it only changed the linear step of the models without modifying
other parts.

Furthermore, to further improve model prediction performance, a response-
history filter (as shown in Figure 2) and flexible nonlinearity can be added
to the model [28, 17], as both of which can also be parameterized by a 1D
spline basis. With such extra flexibility, a model can better adapt to the
cell-specific firing threshold and adaptation to the stimulus ensemble. Still,
even with the response-history filter, a STRF model may still not be able to
account for the total variability of the neuronal responses, if the response
is influenced by other sources. Future extension of the current GLM for STRFs
can consider incorporating more diverse sources of behavioral data [29] of the
experimental animal under experimental or naturalistic settings as 1D time-
dependent filters.

##  4 Material and Methods

###  4.1 Data sets

#### Simulated data

For the simulated data shown in Figure 3 and 4, we used a 2D rank-2 STRF with
30 time frames and 40 pixels in space ($d=1200$) as shown as "ground truth" in
both Figures and generated responses with two types of flicker bar stimuli:
Gaussian white noise and pink noise. For LNP models, we used a fixed
exponential nonlinearity.

For Figure 5, we used two 20x20 pixels antagonistic center-surround STRF and
with Gaussian white noise as stimulus. We used a fixed exponential function
for the filter nonlinearity and a softplus function for the output
nonlinearity. We controlled the mean firing rate of all simulated spike trains
to about 21 Hz by adjusting the intercept of the corresponding models (for
Figure 4A and C, the intercepts were 3.5 and 2.5, respectively. For Figure 5,
the intercept was 0) and also multiplying a constant $R=10$ to $\lambda$.

For all simulations, we measured the similarity of the ground truth and the
model estimators by computing the MSE between the respective normalized STRFs.
The STRFs were normalized by element-wise division by the matrix-norm.

For Figure 3, the similarity was measured for various stimulus length
($n_{\mathrm{samples}}$), which were scaled by different factors
($\mathrm{factor}\in{0.5,1,2,4,8,16,32,64}$) with respect to the number of
coefficients of the STRF. For Figure 4 and 5, the similarity was measured at
various lengths of time ($\mathrm{time}\in{0.5,1,2,4,8,16,32,64}$ minutes),
which is equal to setting $\mathrm{time}=n_{\mathrm{samples}}\times\Delta$,
where $\Delta$ is the time bin size and took the value of 0.033 s.

Error bars in Figure 3-5 were obtained by averaging over 10 trials of
simulations with different random seeds for the stimulus.

For the simulated data shown in Figure 7, we used a 3D STRF with 25 time
frames and 25x25 pixels in space ($d=25\times 25\times 25=15625$), where the
spatial profile was a 2D Gaussian field and the temporal profile had a sharp
dip close to zero time lag, and generated responses with Gaussian white noise
(sample size $n=d$ unless noted otherwise). Additive Gaussian noise was added
to the responses ($\sigma=1$ except (c), for which $\sigma=3$).

#### Experimental data

We used a previously published data set for Figure 9A: extracellularly
recorded spikes from tiger salamander retinal ganglion cells [14], available
to download at
https://gin.g-node.org/gollischlab/Liu_etal_2017_RGC_spiketrains_for_STNMF.

Another two data sets of visual neurons and one data set of a non-visual
neuron were used for the supplementary Figure 2 and 3:

  1. 1.

Extracellularly recorded spikes from a tiger salamander retinal ganglion cell
[17], available to download at https://github.com/baccuslab/inferring-hidden-
structure-retinal-circuits;

  2. 2.

Extracellular recorded spikes from a V1 complex cell of macaque monkey [30],
available on request from the authors;

  3. 3.

Extracellular recorded spikes from a grid cell in rat entorhinal cortex [31],
available to download at https://www.ntnu.edu/kavli/research/grid-cell-data;

The details of the second data set for Figure 9B are described below.

###  4.2 Experimental procedures for 2-photon calcium imaging

#### Animals and tissue preparation

One wild-type mouse (C57Bl/6J, JAX 000664) housed under a standard 12 hour
day/night rhythm was used in this study. Before experiments, the mouse was
dark adapted for at least 2 hours. Then, the animal was anaesthetized with
isoflurane (Baxter) and killed with cervical dislocation. Eyes were quickly
enucleated in carboxygenated (95% O2, 5% CO2) artificial cerebral spinal fluid
(ACSF) solution containing (in mM): 125 NaCl, 2.5 KCl, 2 CaCl2, 1 MgCl2, 1.25
NaH2PO4, 26 NaHCO3, 20 glucose, and 0.5 L-glutamine (pH 7.4). After carefully
removed the cornea, sclera and vitreous body, the retina was flattened on an
Anodisc (0.2 $\mu$M pore size, GE Healthcare) with the RGC side facing up and
electroporated (see below). Then, the retina was transferred to the recording
chamber of the microscope. To visualize the blood vessels and the damaged
cells, we added 5 $\mu$M sulforhodamine-101 (SR101, Invitrogen) into 2l of
ACSF solution. All experimental procedures were carried out under very dim red
light. All animal procedures were approved by the governmental review board
(Regierungspräsidium Tübingen, Baden-Württemberg, Konrad-Adenauer-Str. 20,
72072 Tübingen, Germany) and performed according to the laws governing animal
experimentation issued by the German Government.

#### Electroporation

For two-photon $\mathrm{Ca}^{2+}$ imaging, Oregon-Green BAPTA-1 (OGB-1,
hexapotassium salt; Life Technologies) was bulk electroporated [32]. After the
retina was flattened on an Anodisc and placed between two 4-mm horizontal
plate electrodes (CUY700P4E/L, Nepagene/Xceltis), 10 $\mu$l OGB-1 (5 mM)
diluted in ACSF solution was suspended on the upper electrode. Then the upper
electrode was lowered onto the retina and applied 10 pulses ( 9.2 V, 100 ms
pulse width, at 1 Hz) from a pulse generator/wide-band amplifier combination
(TGP110 and WA301, Thurlby handar/Farnell). For complete recovery of the
electroporated retina, we started recordings 30 min after electroporation.

#### Two photon imaging and light stimulation

A MOM-type two-photon microscope (designed by W. Denk, MPI, Martinsried;
purchased from Sutter Instruments/Science Products) as described previously
[33, 34] was used for this study. Briefly, the system was equipped with a
mode-locked Ti:Sapphire laser (MaiTai-HP DeepSee, Newport Spectra-Physics,
Darmstadt, Germany), green and red fluorescence detection channels for OGB-1
(HQ 510/84, AHF, Tübingen, Germany) and SR101 (HQ 630/60, AHF), and a water
immersion objective (16x/0.80W, DIC N2, $\infty$/0 WD 3.0, Nikon). For
$\mathrm{Ca}^{2+}$ imaging, we tuned the laser to 927 nm, and used a custom-
made software (ScanM, by M. Müller, MPI, Martinsried, and T.E.) running under
IGOR Pro 6.3 for Windows (Wavemetrics). Time-elapsed $\mathrm{Ca}^{2+}$
signals were recorded with 64$\times$16 pixel image sequences (31.25 Hz).

Light stimuli were projected through the objective lens [34]. A digital light
processing (DLP) LightCrafter 4500 (Texas Instruments, Dallas, TX) coupled
with external unit with light-emitting diodes (LEDs) – "green" (576 nm) and UV
(387 nm) that match the spectral sensitivities of mouse M- and S-opsins (for
details, see [23, 35]) was used in this study. Both LEDs were intensity-
calibrated to range from $\mathrm{0.1\cdot 10}^{3}$ ("black" background) to
$\mathrm{20.0\cdot 10}^{3}$ ("white" full field) photoisomerizations
$\mathrm{P^{*}s^{-1}}$ per cone. The light stimulus was centered before every
experiment, ensuring that its center corresponded to the center of the
microscope’s scan field. For all experiments, the tissue was kept at a
constant mean stimulator intensity level for $\geq$ 15 s after the laser
scanning started and before light stimuli were presented.

Binary dense noise (20$\times$15 matrix of 30$\times$30 $\mu$M pixels; each
pixel displayed an independent, balanced random sequence at 5 Hz for 20
minutes) were generated and presented using the Python-based software package
QDSpy (https://github.com/eulerlab/QDSpy).

#### Data pre-processing

$\mathrm{Ca}^{2+}$ imaging data were pre-processed using custom-written
scripts in IGOR Pro. Regions-of-interest (ROIs) were manually defined
depending on the outline of the soma. For each ROI, $\mathrm{Ca}^{2+}$ traces
were extracted using the image analysis toolbox SARFIA for IGOR Pro [36]. A
stimulus time marker embedded in the recorded data served to align the traces
relative to the light stimulus at a temporal resolution of 2 ms. All stimulus-
aligned traces together with the relative ROI position were exported for
further analysis.

###  4.3 Previous methods

#### Spike-triggered Average (STA)

The STA was computed as the average of all stimuli that proceeded a spike [1]:

$w_{STA}=\frac{1}{n}X^{T}y$ (1)

where $X$ is the stimulus design matrix, $y$ the response, $n$ can be the
number of spikes or the number of samples (rows) in $X$ and $y$ for non-spike
data.

#### Whitened spike-triggered average (wSTA)

The maximum-likelihood estimator under a Gaussian noise model, or whitened-STA
(wSTA), was computed by multiplying the STA with the inverse of the stimulus
auto-correlation matrix:

$w_{wSTA}=(X^{T}X)^{-1}X^{T}y$ (2)

In our code base, we did not compute the matrix inverse explicitly but
computed the least squared solution instead, which is common practice.

#### Maximum-a-posteriori estimator (MAP)

Similar, the MAP under a Gaussian noise model was computed by adding the
inverse of the prior covariance to the stimulus auto-correlation matrix:

$w_{MAP}=\left(X^{T}X+C^{-1}\right)^{-1}X^{T}y$ (3)

where $C$ is a prior covariance matrix of choice, e.g. corresponding to a
smoothness or sparseness assumption [6, 7]. Its hyperparameters can be
optimized by minimizing the negative log-evidence[6, 7]:

$\mathrm{Cost}_{MAP}=\frac{n}{2}\mathrm{log}|2\pi\sigma^{2}|+\frac{1}{n}\mathrm{log}|C\Lambda^{-1}|-\frac{1}{n}\mu^{T}\Lambda\mu+\frac{1}{2\sigma^{2}}Y^{T}Y$
(4)

where $\sigma$ is variance of the additive Gaussian noise, $\mu$ and $\Lambda$
are posterior mean and covariance, respectively:

$\mu=\frac{1}{\sigma^{2}}\Lambda X^{T}y$ (5)

and

$\Lambda=\left(\frac{1}{\sigma^{2}}X^{T}X+C^{-1}\right)^{-1}$ (6)

#### Smoothness prior covariance matrix for ASD

The 1D ASD prior covariance matrix is formulated as following [6]:

$C_{ij}=\exp\left(-\rho-\frac{\Delta_{ij}}{2\delta^{2}}\right)$ (7)

where $\Delta_{ij}$ is the squared distance between STRF coefficients $w_{i}$
and $w_{j}$, $\rho$ controls the scale and $\delta$ controls the smoothness.

#### Locality prior covariance matrix for ALD

The 1D ALD prior covariance matrix is formulated as following [7]:

$C=C_{s}^{\frac{1}{2}}(B^{T}C_{f}B)C_{s}^{\frac{1}{2}}$ (8)

where $B$ is an orthogonal basis matrix for the 1D discrete Fourier transform,
$C_{s}$ and $C_{f}$ are the diagonal ALD prior covariance in space and
frequency, respectively, which can be constructed as:

$C_{s,ii}=\exp\left(-\frac{1}{2\tau_{s}^{2}}(\chi_{i}-\nu_{s})^{2})\right)$
(9)

where $\tau_{s}$ is scalar values control the shape and extent of the local
region, $\chi$ is the STRF coefficient coordinates in space-time and $\nu_{s}$
is center coordinate of the STRF. And

$C_{f,ii}=\exp\left(-\frac{1}{2}(|\tau_{f}\omega_{i}|-\nu_{f})^{2})\right)$
(10)

where $\tau_{s}$ is scalar values control the scale of smoothness in
frequency, $\omega$ is of STRF coefficient coordinates in frequency and
$\nu_{f}$ is center coordinate of the STRF in Fourier domain.

High-dimensional prior covariance matrix can be approximated by taking the
Kronnecker product ($\otimes$) of the covariance matrix in each dimension.

###  4.4 Estimating STRFs with Natural Cubic Regression Splines

To generate a natural cubic spline matrix ($S$), we used the implementation
from Patsy (0.5.1), a Python equivalent of the original implementation in R
package mcgv [11]. The only hyperparameter needed to generate a 1D spline
basis matrix was the number of basis functions, or the degrees of freedom
(df), as we set all knots to be equally spaced.

$S=a^{-}_{j}(x)+a^{+}_{j}(x)+c^{-}_{j}(x)F_{j}+c^{+}_{j}(x)F_{j+1}\;\mathrm{if}\;x_{j}<x\leq
x<x_{j+1}$ (11)

where $S$ is the spline basis matrix, $x$ is knot locations, $F^{-1}=B^{-1}D$,
$F=(0,F^{-1},0)$ and $a^{-}_{j}$, $a^{+}_{j}$, $c^{-}_{j}$, $c^{+}_{j}$, $B$,
$D$ are defined in Table 2.

$h_{j}=x_{j+1}-x_{j}$ |  |
---|---|---
$a^{-}_{j}(x)=(x_{j+1}-x)/h_{j}$ | $c^{-}_{j}(x)=[(x_{j+1}-x)^{3}/h_{j}-h_{j}(x_{j+1}-x)]/6$ |
$a^{+}_{j}(x)=(x-x_{j})/h_{j}$ | $c^{+}_{j}(x)=[(x-x_{j})^{3}/h_{j}-h_{j}(x-x_{j})]/6$ |
$D_{i,i}=1/h_{i}$ | $D_{i,i+1}=-1/h_{i}-1/h_{i+1}$ | $D_{i,i+2}=1/h_{i+1}$
$B_{i,i}=(h_{i}+h_{i+1})/3$ |  | i=1…k-2
$B_{i,i+1}=h_{i+1}/6$ | $B_{i+1,i}=h_{i+1}/6$ | i=1…k-3
Table 2: Basis functions for natural cubic spline.

An illustration of 1D and 2D natural cubic splines can be seen in Figure 1B. A
spline-based STRF under a Gaussian noise model can be fitted in the similar
manner as the wSTA:

$b_{SPL,wSTA}=(S^{T}X^{T}XS)^{-1}S^{T}X^{T}y$ (12)

where $b$ is the spline coefficients, and the corresponding STRF is:

$w_{SPL,wSTA}=Sb_{SPL,wSTA}=S(S^{T}X^{T}XS)^{-1}S^{T}X^{T}y$ (13)

The spline matrix can be extended to high dimensions, which is also called
tensor product smooth [11], by simply taking the Kronnecker product
($\otimes$) of the basis functions in each dimension:

$S=S_{t}\otimes S_{x}\otimes S_{y}$ (14)

###  4.5 Spline-based GLMs with gradient descent

The spline coefficients can also be fitted by gradient descent with respect to
the cost function defined in each model with L1 regularization.

For the LG model, the cost function is the sum square error:

$\mathrm{Cost}_{LG}=\frac{1}{n}\sum(y-XSb)^{2}+\alpha\|b\|_{1}$ (15)

where $\alpha$ is the weight of L1 regularization.

For LNP and LNLN model, the cost function is the negative log-likelihood:

$\begin{split}\lambda=f(XSb),\\\
\mathrm{Cost}_{LNP/LNLN}=-ylog(\lambda)+\sum\Delta\lambda+\alpha\|b\|_{1}\end{split}$
(16)

where $\lambda$ is the conditional intensity or the instantaneous firing rate,
$\Delta$ is the size of time bin, and $f(.)$ is a nonlinearity applied to the
filter output, which can be fixed as a softplus or exponential nonlinearity.

We used JAX [37] for automatic differentiation, and Adam optimizer [38]for
parameters optimization.

###  4.6 Spline-based spike-triggered clustering methods

Both k-means clustering and semi-NMF can be formulated as a matrix
factorization problem:

$\mathrm{min}\|V_{\pm}-W_{\pm}H_{+}^{T}\|^{2}$ (17)

where $V$ is the spike-triggered ensemble, $W$ is clustering centroids matrix
(or subunit matrix), and $H$ is the clustering label (or subunit weight
matrix). A spline basis can be easily incorporated into k-Means clustering by
assuming that each centroid can be approximated by the same spline basis
matrix: $W_{SPL,centroid}=SB$, where $B$ is the spline coefficient matrix for
each subunit, and can be computed as the least-squares solution to the linear
matrix equation $B=S^{\dagger}W_{centroid}$, where $\dagger$ is pseudo-
inverse. For semi-NMF, splines can also be combined with the update rules
[18]:

$W=SS^{\dagger}VH(H^{T}H)^{-1}$ (18a)
$H=H*\sqrt{\frac{(V^{T}W)^{+}+[H(W^{T}W)^{-}]}{(V^{T}W)^{-}+[H(W^{T}W)^{+}]}}$
(18b)

where $*$ is point-wise multiplication, $A^{+}=(\mid A\mid+A)/2$ and
$A^{-}=(\mid A\mid-A)/2$.

###  4.7 Confidence interval and statistical significance of the estimated
STRFs

The asymptotic posterior probability of the estimated STRF coefficients can be
computed analytically [11]: $b|y,X\sim N(b,V_{b})$, where $b$ is the fitted
STRF coefficients, $V_{b}$ is the posterior covariance. For LG model:

$V_{b}=(S^{T}X^{T}XS)^{-1}\sigma^{2}$ (19)

where $\sigma^{2}$ is the residual sum of squares of the data and predicted
response from the training set. And for LNP model:

$V_{b}=(S^{T}X^{T}WXS)^{-1}$ (20)

where $W$ is a diagonal matrix with entries $w_{ii}=\frac{1}{f(XSb)^{2}}$, and
$f(.)$ is the nonlinearity of the model. Thus, the standard error of $b$ is
the square root of the diagonal of $V_{b}$, and the corresponding confidence
interval (CI) of the estimated STRF is as following:

$b-1.96\times b_{se}\leq CI\leq b+1.96\times b_{se}$ (21)

Given the computed confidence interval, Wald Chi-Squared Test can be used to
determine whether a STRF is statistically significant (against the null
hypothesis that all STRF coefficients are zero), with the Wald statistic
$T=bV_{b}^{-1}b^{T}$.

###  4.8 Permutation Test

To evaluate the predictive performance of the estimated STRF, we permuted the
validation test stimulus in time and made prediction over 100 repetitions. A
one-tailed one-sample Student’s T-test was used to evaluate if the model
prediction was better than chance level.

###  4.9 Visualization of 3D STRF

For visualization of 3D STRF in Figure 7 and 9, we first performed singular
value decomposition (SVD) to separate the spatial and temporal components of
the estimated 3D STRF. The coordinates of the extremum pixel in the spatial
component were selected, and the temporal filter of the original STRF is shown
for those coordinates. The confidence interval of the same extremum point was
plotted as a gray shaded region on the temporal filter. For the spatial
filter, we show the frame of the original STRF at the extremum point of the
temporal filter (indicated by a vertical dotted line). The maximum value of
the color map was set to the maximum absolute value of the STRF.

###  4.10 Comparison of Computation Time

For comparison of computation time in Figure 6, we used simulated 2D RFs of
different sizes (15x15, 20x20 and 30x30 pixels, respectively) and white noise
stimuli. Simulations were carried out at different lengths of time in the same
manner as Figure 4 and 5. The number of basis functions for SPL wSTA were
fixed as 5 and 10 along each dimension in all three cases. Measurements of
computation time were averaged over 10 repetitions and performed on a 2019
16-inch MacBook Pro with a 2.3 GHz 8-core Intel Core i9 and 16 GB of RAM.

###  4.11 Data availability

All data used are available in GitHub
(https://github.com/huangziwei/data_RFEst).

###  4.12 Code availability

All methods (spline-based GLMs, evidence optimization and spike-clustering
methods) were implemented in Python3 and available in GitHub
(https://github.com/berenslab/RFEst). And Jupyter notebooks for all figures
are available also in GitHub (https://github.com/huangziwei/notebooks_RFEst)

## Acknowledgements

We thank the authors of [17, 14, 30] and [31] for making their data available.
Research was funded by the Deutsche Forschungsgemeinschaft through a
Heisenberg Professorship (BE5601/4-1 and BE5601/8-1, PB), the Collaborative
Research Center 1233 “Robust Vision” (reference number 276693517) as well as
individual research grants (EU 42/10-1, BE5601/6-1), the German Ministry of
Education and Research through the Bernstein Award (01GQ1601, PB). PB is a
member of the Excellence Cluster 2064 “Machine Learning — New Perspectives for
Science” (reference number 390727645) and the Tübingen AI Center (FKZ:
01IS18039A).

## References

  * Aljadeff et al. [2016] Aljadeff J, Lansdell BJ, Fairhall AL, Kleinfeld D.  Analysis of Neuronal Spike Trains, Deconstructed.  Neuron 2016 Jul;91(2):221–259.  http://www.sciencedirect.com/science/article/pii/S0896627316302501.
  * Schwartz et al. [2006] Schwartz O, Pillow JW, Rust NC, Simoncelli EP.  Spike-triggered neural characterization.  Journal of Vision 2006 Feb;6(4):13–13.  http://jov.arvojournals.org/article.aspx?articleid=2192881, publisher: The Association for Research in Vision and Ophthalmology.
  * Paninski [2003] Paninski L.  Convergence properties of three spike-triggered analysis techniques.  Network: Computation in Neural Systems 2003;14(3):437–464.  Publisher: Taylor & Francis.
  * Theunissen et al. [2001] Theunissen FE, David SV, Singh NC, Hsu A, Vinje WE, Gallant JL.  Estimating spatio-temporal receptive fields of auditory and visual neurons from their responses to natural stimuli.  Network (Bristol, England) 2001 Aug;12(3):289–316.
  * Wu et al. [2006] Wu MCK, David SV, Gallant JL.  Complete functional characterization of sensory neurons by system identification.  Annu Rev Neurosci 2006;29:477–505.  Publisher: Annual Reviews.
  * Sahani and Linden [2002] Sahani M, Linden J.  Evidence optimization techniques for estimating stimulus-response functions.  Advances in neural information processing systems 2002;15:317–324.
  * Park and Pillow [2011] Park M, Pillow JW.  Receptive field inference with localized priors.  PLoS Comput Biol 2011;7(10):e1002219.  Publisher: Public Library of Science.
  * Casella [1985] Casella G.  An Introduction to Empirical Bayes Data Analysis.  The American Statistician 1985 May;39(2):83–87.  https://www.tandfonline.com/doi/abs/10.1080/00031305.1985.10479400, publisher: Taylor & Francis _eprint: https://www.tandfonline.com/doi/pdf/10.1080/00031305.1985.10479400.
  * Pillow et al. [2008] Pillow JW, Shlens J, Paninski L, Sher A, Litke AM, Chichilnisky EJ, et al.  Spatio-temporal correlations and visual signalling in a complete neuronal population.  Nature 2008 Aug;454(7207):995–999.  https://www.nature.com/articles/nature07140/, number: 7207 Publisher: Nature Publishing Group.
  * Thorson et al. [2015] Thorson IL, Liénard J, David SV.  The Essential Complexity of Auditory Receptive Fields.  PLOS Computational Biology 2015 Dec;11(12):e1004628.  https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1004628, publisher: Public Library of Science.
  * Wood [2017] Wood SN.  Generalized additive models: an introduction with R.  CRC press; 2017.
  * Pillow et al. [2005] Pillow JW, Paninski L, Uzzell VJ, Simoncelli EP, Chichilnisky EJ.  Prediction and Decoding of Retinal Ganglion Cell Responses with a Probabilistic Spiking Model.  Journal of Neuroscience 2005 Nov;25(47):11003–11013.  https://www.jneurosci.org/content/25/47/11003, publisher: Society for Neuroscience Section: Behavioral/Systems/Cognitive.
  * McFarland et al. [2013] McFarland JM, Cui Y, Butts DA.  Inferring Nonlinear Neuronal Computation Based on Physiologically Plausible Inputs.  PLOS Computational Biology 2013 Jul;9(7):e1003143.  https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003143, publisher: Public Library of Science.
  * Liu et al. [2017] Liu JK, Schreyer HM, Onken A, Rozenblit F, Khani MH, Krishnamoorthy V, et al.  Inference of neuronal functional circuitry with spike-triggered non-negative matrix factorization.  Nature Communications 2017 Jul;8(1):149.  https://www.nature.com/articles/s41467-017-00156-9, number: 1 Publisher: Nature Publishing Group.
  * Shah et al. [2020] Shah NP, Brackbill N, Rhoades C, Kling A, Goetz G, Litke AM, et al.  Inference of nonlinear receptive field subunits with spike-triggered clustering.  eLife 2020 Mar;9:e45743.  https://doi.org/10.7554/eLife.45743, publisher: eLife Sciences Publications, Ltd.
  * Park and Pillow [2013] Park M, Pillow JW.  Bayesian inference for low rank spatiotemporal neural receptive fields.  In: Advances in Neural Information Processing Systems; 2013. p. 2688–2696.
  * Maheswaranathan et al. [2018] Maheswaranathan N, Kastner DB, Baccus SA, Ganguli S.  Inferring hidden structure in multilayered neural circuits.  PLOS Computational Biology 2018 Aug;14(8):e1006291.  https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1006291, publisher: Public Library of Science.
  * Ding et al. [2010] Ding CHQ, Li T, Jordan MI.  Convex and Semi-Nonnegative Matrix Factorizations.  IEEE Transactions on Pattern Analysis and Machine Intelligence 2010 Jan;32(1):45–55.  Conference Name: IEEE Transactions on Pattern Analysis and Machine Intelligence.
  * Park and Pillow [2011] Park IM, Pillow JW.  Bayesian spike-triggered covariance analysis.  In: Advances in neural information processing systems; 2011. p. 1692–1700.
  * Kimeldorf and Wahba [1970] Kimeldorf GS, Wahba G.  A correspondence between Bayesian estimation on stochastic processes and smoothing by splines.  The Annals of Mathematical Statistics 1970;41(2):495–502.  Publisher: JSTOR.
  * Bay et al. [2016] Bay X, Grammont L, Maatouk H.  Generalization of the Kimeldorf-Wahba correspondence for constrained interpolation.  Electronic journal of statistics 2016;10(1):1580–1595.  Publisher: Institute of Mathematical Statistics and Bernoulli Society.
  * Aoi and Pillow [2017] Aoi MC, Pillow JW.  Scalable Bayesian inference for high-dimensional neural receptive fields.  bioRxiv 2017 Nov;p. 212217.  https://www.biorxiv.org/content/10.1101/212217v1, publisher: Cold Spring Harbor Laboratory Section: New Results.
  * Baden et al. [2016] Baden T, Berens P, Franke K, Román Rosón M, Bethge M, Euler T.  The functional diversity of retinal ganglion cells in the mouse.  Nature 2016 Jan;529(7586):345–350.  https://www.nature.com/articles/nature16468, number: 7586 Publisher: Nature Publishing Group.
  * Ran et al. [2020] Ran Y, Huang Z, Baden T, Schubert T, Baayen H, Berens P, et al.  Type-specific dendritic integration in mouse retinal ganglion cells.  Nature Communications 2020 Apr;11(1):2101.  https://www.nature.com/articles/s41467-020-15867-9, number: 1 Publisher: Nature Publishing Group.
  * Trong and Rieke [2008] Trong PK, Rieke F.  Origin of correlated activity between parasol retinal ganglion cells.  Nature Neuroscience 2008 Nov;11(11):1343–1351.  https://www.nature.com/articles/nn.2199, number: 11 Publisher: Nature Publishing Group.
  * Qiu et al. [2020] Qiu Y, Zhao Z, Klindt D, Kautzky M, Szatko KP, Schaeffel F, et al.  Mouse retinal specializations reflect knowledge of natural environment statistics.  bioRxiv 2020 Dec;p. 2020.12.08.416172.  https://www.biorxiv.org/content/10.1101/2020.12.08.416172v1, publisher: Cold Spring Harbor Laboratory Section: New Results.
  * Taylor and Vaney [2002] Taylor WR, Vaney DI.  Diverse Synaptic Mechanisms Generate Direction Selectivity in the Rabbit Retina.  Journal of Neuroscience 2002 Sep;22(17):7712–7720.  https://www.jneurosci.org/content/22/17/7712, publisher: Society for Neuroscience Section: ARTICLE.
  * Williamson et al. [2015] Williamson RS, Sahani M, Pillow JW.  The Equivalence of Information-Theoretic and Likelihood-Based Methods for Neural Dimensionality Reduction.  PLOS Computational Biology 2015 Apr;11(4):e1004141.  https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1004141, publisher: Public Library of Science.
  * Balzani et al. [2020] Balzani E, Lakshminarasimhan K, Angelaki D, Savin C.  Efficient estimation of neural tuning during naturalistic behavior.  Advances in Neural Information Processing Systems 2020;33.  https://papers.nips.cc/paper/2020/hash/94d2a3c6dd19337f2511cdf8b4bf907e-Abstract.html.
  * Rust et al. [2005] Rust NC, Schwartz O, Movshon JA, Simoncelli EP.  Spatiotemporal elements of macaque v1 receptive fields.  Neuron 2005;46(6):945–956.  Publisher: Elsevier.
  * Hafting et al. [2005] Hafting T, Fyhn M, Molden S, Moser MB, Moser EI.  Microstructure of a spatial map in the entorhinal cortex.  Nature 2005 Aug;436(7052):801–806.  https://www.nature.com/articles/nature03721, bandiera_abtest: a Cg_type: Nature Research Journals Number: 7052 Primary_atype: Research Publisher: Nature Publishing Group.
  * Briggman and Euler [2011] Briggman KL, Euler T.  Bulk electroporation and population calcium imaging in the adult mammalian retina.  Journal of Neurophysiology 2011 Feb;105(5):2601–2609.  https://journals.physiology.org/doi/full/10.1152/jn.00722.2010, publisher: American Physiological Society.
  * Euler et al. [2009] Euler T, Hausselt SE, Margolis DJ, Breuninger T, Castell X, Detwiler PB, et al.  Eyecup scope—optical recordings of light stimulus-evoked fluorescence signals in the retina.  Pflügers Archiv - European Journal of Physiology 2009 Apr;457(6):1393–1414.  https://doi.org/10.1007/s00424-008-0603-5.
  * Euler et al. [2019] Euler T, Franke K, Baden T.  Studying a Light Sensor with Light: Multiphoton Imaging in the Retina.  In: Hartveit E, editor. Multiphoton Microscopy Neuromethods, New York, NY: Springer; 2019.p. 225–250.  https://doi.org/10.1007/978-1-4939-9702-2_10.
  * Franke et al. [2019] Franke K, Maia Chagas A, Zhao Z, Zimmermann MJ, Bartel P, Qiu Y, et al.  An arbitrary-spectrum spatial visual stimulator for vision research.  eLife 2019 Sep;8:e48779.  https://doi.org/10.7554/eLife.48779, publisher: eLife Sciences Publications, Ltd.
  * Dorostkar et al. [2010] Dorostkar MM, Dreosti E, Odermatt B, Lagnado L.  Computational processing of optical measurements of neuronal and synaptic activity in networks.  Journal of Neuroscience Methods 2010 Apr;188(1):141–150.  https://www.sciencedirect.com/science/article/pii/S0165027010000725.
  * Bradbury et al. [2018] Bradbury J, Frostig R, Hawkins P, Johnson MJ, Leary C, Maclaurin D, et al.  JAX: composable transformations of Python+ NumPy programs.  Version 01 2018;55.
  * Kingma and Ba [2017] Kingma DP, Ba J.  Adam: A Method for Stochastic Optimization.  arXiv:14126980 [cs] 2017 Jan;http://arxiv.org/abs/1412.6980, arXiv: 1412.6980.