Modelling multivariate spatio-temporal data with identifiable variational autoencoders

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Modelling multivariate spatio-temporal data with identifiable variational autoencoders © 2024 The Authors. Published by Elsevier Ltd. Published version Sipilä, Mika; Cappello, Claudia; De Iaco, Sandra; Nordhausen, Klaus; Taskinen, Sara Sipilä, M., Cappello, C., De Iaco, S., Nordhausen, K., & Taskinen, S. (2025). Modelling multivariate spatio-temporal data with identifiable variational autoencoders. Neural Networks, 181, Article 106774. https://doi.org/10.1016/j.neunet.2024.106774 2025 Contents lists available at ScienceDirect Neural Networks journal homepage: www.elsevier.com/locate/neunet Full Length Article Modelling multivariate spatio-temporal data with identifiable variational autoencoders Mika Sipilä a,∗, Claudia Cappello b, Sandra De Iacob, Klaus Nordhausen a, Sara Taskinen a aDepartment of Mathematics and Statistics, University of Jyväskylä, Finland bDSE - Section of Mathematics and Statistics, University of Salento, Italy ARTICLE INFO Keywords: Blind source separation Dimension estimation Kriging Meteorological data Shapley values ABSTRACT Modelling multivariate spatio-temporal data with complex dependency structures is a challenging task but can be simplified by assuming that the original variables are generated from independent latent components. If these components are found, they can be modelled univariately. Blind source separation aims to recover the latent components by estimating the unknown linear or nonlinear unmixing transformation based on the observed data only. In this paper, we extend recently introduced identifiable variational autoencoder to the nonlinear nonstationary spatio-temporal blind source separation setting and demonstrate its performance using comprehensive simulation studies. Additionally, we introduce two alternative methods for the latent dimension estimation, which is a crucial task in order to obtain the correct latent representation. Finally, we illustrate the proposed methods using a meteorological application, where we estimate the latent dimension and the latent components, interpret the components, and show how nonstationarity can be accounted and prediction accuracy can be improved by using the proposed nonlinear blind source separation method as a preprocessing method. 1. Introduction Many real world phenomena, such as weather, epidemiological patterns and ecosystem dynamics, are multivariate spatio-temporal, meaning that multivariate observations 𝒙(𝒔, 𝑡) ∶= 𝒙∈R𝑆are observed in a spatial location 𝒔∈⊂R𝐷at time 𝑡∈⊂R, where is called a spatial domain, is called a temporal domain and 𝐷is a spatial dimension. Without loss of generality, we assume from now on that 𝐷= 2. A multivariate observation 𝒙contains measurements of multiple, usually dependent, random variables describing the phenomenon of interest. When modelling such multivariate spatio-temporal data, one has to account not only the dependence between the variables, but also the dependences in space and in time. The dependence structure is often described through spatio-temporal covariance function 𝑪(𝒙(𝒔, 𝑡),𝒙(𝒔′, 𝑡′)), where (𝒔, 𝑡)and (𝒔′, 𝑡′)are two spatio-temporal locations. The covariance 𝑪is a 𝑆×𝑆matrix valued functional with elements 𝐶𝑖𝑗 ,𝑖, 𝑗 = 1,…, 𝑆, defined as 𝐶𝑖𝑗 =𝐶(𝑥𝑖(𝒔, 𝑡), 𝑥𝑗(𝒔′, 𝑡′)) = 𝐸(𝑥𝑖(𝒔, 𝑡)𝑥𝑗(𝒔′, 𝑡′))−𝐸(𝑥𝑖(𝒔, 𝑡))𝐸(𝑥𝑗(𝒔′, 𝑡′)). Modelling the covariance function 𝑪is usually a highly demanding task, and often, in order to make the modelling feasible, some severely restricting assumptions, such as stationarity or separability, are made. When the spatio-temporal field is assumed to be stationary, the covariance function can be simplified ∗Corresponding author. E-mail address: [email protected] (M. Sipilä). to 𝑪(𝒙(𝒔, 𝑡),𝒙(𝒔′, 𝑡′)) = 𝑪(‖𝒔−𝒔′‖,|𝑡−𝑡′|),(1) meaning that the value of the function depends only on the distance between the spatial locations and the distance between temporal locations. If (1) does not hold, the data are nonstationary, meaning that the covariance function may differ when spatial or temporal location is altered. When separability is assumed, the spatio-temporal covariance function can be written as a product of spatial and temporal covariance functions as 𝑪(𝒙(𝒔, 𝑡),𝒙(𝒔′, 𝑡′)) = 𝑪(𝒔,𝒔′)𝑪(𝑡, 𝑡′), meaning that the spatial and temporal covariance models can be fitted independently and that the spatio-temporal interaction is not considered. The assumptions of stationarity or separability often lead to unrealistically simple models that hence produce nonoptimal results under nonstationary or nonseparable data. Accounting complex nonseparable and nonstationary correlation structures is complicated already in the univariate case, for which an overview can be found in Chen, Genton, and Sun (2021). For multivariate data, the task is even more demanding and computationally challenging as the crossdependencies between the variables have to be taken into account. https://doi.org/10.1016/j.neunet.2024.106774 Received 6 June 2024; Received in revised form 11 September 2024; Accepted 29 September 2024 Neural Networks 181 (2025) 106774 Available online 9 October 2024 0893-6080/© 2024 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). M. Sipilä et al. For more details of complexity of nonstationary covariance functions for multivariate spatio-temporal data, see Porcu, Furrer, and Nychka (2021), Salvana and Genton (2020). Another approach to simplify the modelling is to assume that the observations are composed of 𝑃latent, mutually independent components (ICs) 𝒛(𝒔, 𝑡) ∶= 𝒛∈R𝑃through some mixing environment. The main motivation for assuming the ICs is, that if the latent components 𝒛are recovered, they can be modelled univariately making for example nonstationarity much easier to account for. Being able to model components univariately is especially desirable in spatiotemporal settings, where multivariate modelling is highly demanding and computationally challenging as discussed previously. Additionally, the ICs may reveal some meaningful patterns and structures in the observed data that can lead to new insights of the phenomenon of interest. A linear blind source separation (BSS) (Comon & Jutten,2010) is a popular approach to recover the latent components 𝒛. In linear BSS, it is assumed that the mixing environment is linear and usually also that 𝑆=𝑃meaning that a 𝑃-variate observable random vector 𝒙= (𝑥1,…, 𝑥𝑃)⊤is generated as 𝒙=𝑨𝒛,(2) where 𝑨is an invertible 𝑃×𝑃mixing matrix and 𝒛= (𝑧1,…, 𝑧𝑃)⊤are the 𝑃-variate latent components. The objective is to recover 𝑨and 𝒛 using only 𝒙and varying assumptions on 𝒛depending on the method used. For example, spatial BSS (SBSS) (Bachoc, Genton, Nordhausen, Ruiz-Gazen, & Virta,2020;Nordhausen, Oja, Filzmoser, & Reimann, 2015), which is a method for multivariate stationary spatial data, assumes spatially stationary 𝒛, and a nonstationary extension of SBSS, spatial nonstationary source separation (SNSS) (Muehlmann, Bachoc, & Nordhausen,2022), assumes 𝒛to have nonstationary spatial covariance function. SBSS and SNSS recover the latent components by jointly diagonalizing two or more moment-based matrices. Recently, SBSS was also extended for stationary spatio-temporal data yielding spatiotemporal BSS (STBSS) (Muehlmann, De Iaco, & Nordhausen,2023). A drawback of STBSS and linear BSS methods in general is that they assume linear mixing (2) which may be too restrictive assumption for many real life applications. Similarly, the assumption that there are as many latent ‘‘signal’’ components as observed variables is in many applications undesirable and it is often hoped that there are significantly fewer signals. This assumption is often needed simply due to the lack of tools for estimating the correct number of signals. Finally, STBSS is developed only for stationary data, and to our knowledge, there are no spatio-temporal alternatives available for nonstationary data cases. Recent advancements in unsupervised deep learning, such as variational autoencoders (VAEs) (Kingma & Welling,2013) and generative adversarial networks (GANs) (Goodfellow et al.,2020), have increased interest for developing nonlinear BSS methods, where the mixing function is not restricted to be linear, but can be any injective function 𝒇∶R𝑃→R𝑆, which generates the observed data 𝒙as 𝒙=𝒇(𝒛).(3) The objective is then to identify an unmixing transformation 𝒒∶R𝑆→ R𝑃, which returns the latent components 𝒛as 𝒛=𝒒(𝒙) based on the observations 𝒙only. Without any additional assumptions on the mixing transformations 𝒇or on the latent components 𝒛, the model is unidentifiable as there exists infinite nonlinear transformations to generate mutually independent components from the observations (Hyvärinen & Pajunen,1999). For this reason both VAEs and GANs, in general, suffer from the unidentifiability issue. However, in many recent studies (Hälvä & Hyvärinen,2020;Hyvärinen & Morioka, 2016,2017;Hyvärinen, Sasaki, & Turner,2019;Khemakhem, Kingma, Monti, & Hyvärinen,2020) the identifiability have been achieved by introducing some constraints on the distribution of the latent components 𝒛. The main assumption leading to identifiability is that the components 𝑧1,…, 𝑧𝑃are statistically dependent on a 𝑚-dimensional auxiliary variable 𝒖, and that the components are conditionally independent yielding the joint distribution 𝑝(𝒛|𝒖) = ∏𝑃 𝑖=1 𝑝(𝑧𝑖|𝒖). In previous studies, the main focus has been in time series data for which several algorithms and examples for auxiliary variables exist in the literature. In case of stationary time series data, permutation contrastive learning (PCL) (Hyvärinen & Morioka,2017) can be used, for which 𝒖is usually given by one or more previous observations in time. For nonstationary time series data, the available methods are time contrastive learning (TCL) (Hyvärinen & Morioka,2016), hidden Markov nonlinear ICA (HM-NICA) (Hälvä & Hyvärinen,2020) and temporal identifiable VAE (iVAE) (Khemakhem et al.,2020), all of which use the time segment of the observation as 𝒖. Generalized contrastive learning (Hyvärinen et al., 2019) and nonlinear ICA with switching linear dynamical systems (𝛥SNICA) (Hälvä, Le Corff, Lehéricy, So, Zhu, Gassiat, & Hyvärinen,2021) can account both stationary and nonstationary time series. In HM-NICA and 𝛥-SNICA, the auxiliary variables 𝒖are not explicitly provided by the user, but they are instead assumed to be hidden states that are modelled simultaneously by the algorithms. In Sipilä, Nordhausen and Taskinen (2024), iVAE was studied further and extended to nonstationary spatial setting, where spatial segmentation was used as 𝒖.Hälvä et al. (2021) also introduced a structured nonlinear ICA framework which could be used for spatial process, but did not provide any algorithm for the method. In addition to these more general identifiable nonlinear BSS methods, many other deep learning based BSS methods (Ansari, Alatrany, Alnajjar, Khater, Mahmoud, Al-Jumeily, & Hussain,2023) have been introduced for mainly acoustic signal specific settings, where only serial dependence is present. However, the spatio-temporal data as discussed in this paper is special in sense that in the temporal domain there is natural direction of dependence (past-future) while in the spatial domain such direction is missing and the dependence is usually considered as a function of the distance between two points. Hence, none of the previous methods are directly applicable or optimal for such spatio-temporal data. Note that regularly spaced spatio-temporal data is often represented as tensor data, and BSS methods developed for such specific cases, like those in Virta and Nordhausen (2017), are generally not applicable to broader spatio-temporal settings. In particular, we are interested in iVAE, which utilizes the auxiliary variable to make VAE identifiable. iVAE is capable of estimating nonlinear injective mixing function, meaning that it allows the latent dimension 𝑃to be less or equal to the observed dimension 𝑆. However, the latent dimension 𝑃has to be estimated beforehand, and currently the nonlinear BSS framework lacks methods for the latent dimension estimation. In this paper, iVAE is extended to nonstationary spatio-temporal setting by introducing three novel approaches to construct the auxiliary variables. The proposed methods address two key limitations of previous STBSS approaches: they accommodate nonlinear mixing functions and allow for more observed variables than latent components. Moreover, the developed methods are suitable for nonstationary data, unlike earlier STBSS methods, which rely on the assumption of stationarity. The three developed methods, coordinate based, segmentation based and radial basis function based iVAE algorithms, are studied using comprehensive simulation studies to find how various types of nonstationarity affect the performance of the methods. The best performing method, radial basis function based iVAE, is illustrated in real life meteorological application where the recovered latent components are interpreted, and a new procedure to account for nonstationarity in modelling and predicting multivariate data is demonstrated. Moreover, nonlinear BSS framework currently lacks methods for estimating the latent dimension 𝑃, which is a crucial task in order to recover the true latent components and to obtain as low dimensional representation of the data as possible without losing much information. Therefore, two alternative procedures for latent dimension estimation are introduced. To conclude, the main contributions of this paper are: Neural Networks 181 (2025) 106774 2 M. Sipilä et al. Fig. 1. Schematic representations of VAE (left) and iVAE (right) models. For VAE, the lower bound of the data log likelihood (ELBO) is formed of 𝒙,𝒙′,𝒛′,𝝁𝒛|𝒙and 𝝈𝒛|𝒙. In iVAE, ELBO has in addition 𝝁𝒛|𝒖and 𝝈𝒛|𝒖which are provided by the auxiliary function. The latent components are obtained as 𝝁𝒛|𝒙,𝒖. 1. Extending iVAE to the nonstationary spatio-temporal setting by proposing three novel approaches for constructing auxiliary variables. 2. Introducing two alternative procedures for latent dimension estimation. 3. Developing a new iVAE-based method for addressing nonstationarity in the modelling and prediction of spatio-temporal data. The rest of the paper is organized as follows. In Section 2we review basic theory behind VAE and iVAE, and discuss the identifiability, after which the spatio-temporal iVAE extensions are introduced in Section 3. In Section 4, the introduced methods are compared using simulation studies, and two alternative latent dimension estimation methods are studied. Finally, Section 5shows a real data example and Section 6 concludes the paper. 2. Variational autoencoders and identifiability Let 𝒙∈R𝑆be an observable random vector and 𝒛∈R𝑃,𝑃≤𝑆, be a latent random vector, i.e., a source vector. Variational autoencoders (VAE) (Kingma & Welling,2013) assume that the observed data are generated from a deep latent variable model with the structure 𝑝∗(𝒙,𝒛) = 𝑝∗(𝒙|𝒛)𝑝∗(𝒛), where 𝑝∗is a true, unknown generative distribution, 𝒛∼𝑝∗(𝒛)and 𝒙∼𝑝∗(𝒙|𝒛). The distribution of the observed data is then obtained as 𝑝∗(𝒙) = ∫𝑝∗(𝒙,𝒛)𝑑𝒛. VAE consists of an encoder 𝒈(𝒙)and a decoder 𝒉(𝒛), which are parameterized by deep neural networks with parameters 𝜽= (𝜽⊤ 𝒈,𝜽⊤ 𝒉)⊤. The encoder maps the observed data to mean vector 𝝁𝒛|𝒙∈R𝑃and variance vector 𝝈𝒛|𝒙∈R𝑃, which are used to sample a new latent representation 𝒛′by applying the reparametrization trick (Kingma & Welling,2013). The decoder transforms the latent representation 𝒛′ back to the observable data 𝒙′. The VAE framework allows effective optimization of the parameters 𝜽so that after optimization we have that 𝑝𝜽(𝒙) ≈ 𝑝∗(𝒙). The VAE framework learns the full generative model 𝑝𝜽(𝒙,𝒛) = 𝑝𝜽(𝒙|𝒛)𝑝𝜽(𝒛)and a variational approximation 𝑞𝜽(𝒛|𝒙)of the posterior distribution 𝑝𝜽(𝒛|𝒙)by maximizing the lower bound of the data loglikelihood, or evidence lower bound (ELBO), defined as (𝜽|𝒙)≥𝐸𝑞𝜽(𝒛|𝒙)(log 𝑝𝜽(𝒙|𝒛) + log 𝑝𝜽(𝒛) − log 𝑞𝜽(𝒛|𝒙)) with respect to the parameter vector 𝜽. The problem however is that the model is not identifiable, meaning that even though we have a good estimate of the marginal distribution 𝑝∗(𝒙), there is no guarantee that 𝑝𝜽(𝒙,𝒛) ≈ 𝑝∗(𝒙,𝒛). More formally, the model is identifiable if for all (𝒙,𝒛)it holds that ∀(𝜽,𝜽′) ∶ 𝑝𝜃(𝒙) = 𝑝𝜃′(𝒙)⟹𝑝𝜽(𝒙,𝒛) = 𝑝𝜽′(𝒙,𝒛). This means that if we find a parameter vector 𝜽for which 𝑝𝜃(𝒙) = 𝑝∗(𝒙), we also have that 𝑝𝜽(𝒙,𝒛) = 𝑝∗(𝒙,𝒛). This leads to the fact that we have found the correct source density distribution 𝑝𝜽(𝒛) = 𝑝∗(𝒛)and correct conditional distributions 𝑝𝜽(𝒙|𝒛) = 𝑝∗(𝒙|𝒛)and 𝑝𝜽(𝒛|𝒙) = 𝑝∗(𝒛|𝒙). The whole VAE model is illustrated in Fig. 1. In the nonlinear BSS framework, the identifiability has been recently achieved by assuming that the latent sources 𝒛have a conditional distribution 𝑝(𝒛|𝒖), where 𝒖∈R𝑚is an auxiliary variable. The auxiliary variable can for example be previous observations in time (Hyvärinen & Morioka,2017) or current time index (Hälvä & Hyvärinen,2020; Hyvärinen & Morioka,2016;Hyvärinen et al.,2019). Similarly, by assuming that the true latent generating model has the form 𝑝∗(𝒙,𝒛|𝒖) = 𝑝∗(𝒙|𝒛)𝑝∗(𝒛|𝒖),(4) the identifiability can be achieved in the VAE framework, yielding identifiable VAE (iVAE) (Khemakhem et al.,2020). In iVAE, the distribution 𝑝∗(𝒙|𝒛)is defined as 𝑝∗(𝒙|𝒛) = 𝑝∗ 𝝐(𝒙−𝒇(𝒛)), which means that 𝒙can be decomposed into 𝒙=𝒇(𝒛) + 𝝐, where 𝝐is an independent noise vector with density 𝑝𝝐. Assuming the nonnoisy nonlinear BSS model (3), the distribution 𝑝𝝐can be modelled with Gaussian distribution with infinitesimal variance. The function 𝒇∶R𝑃→R𝑆is an injective, but possibly nonlinear function. The conditional distribution of latent sources 𝒛is assumed to be a part of the exponential family, that is, 𝑝𝑻,𝝀(𝒛|𝒖) = 𝑃 ∏ 𝑖=1 𝑄𝑖(𝑧𝑖) 𝑍𝑖(𝒖)exp [𝑘 ∑ 𝑗=1 𝑇𝑖,𝑗 (𝑧𝑖)𝜆𝑖,𝑗 (𝒖)],(5) where 𝑄𝑖(𝑧𝑖)is a base measure, 𝑍𝑖(𝒖)is a normalizing constant, 𝑻𝑖(𝑧𝑖) = (𝑇𝑖,1(𝑧𝑖),…, 𝑇𝑖,𝑘(𝑧𝑖))⊤contains sufficient statistics, and 𝝀𝑖(𝒖)=(𝜆𝑖,1(𝒖),…, 𝜆𝑖,𝑘(𝒖))⊤contains the parameters depending on 𝒖. The dimension 𝑘of each sufficient statistic 𝑻𝑖(𝑧𝑖)and 𝝀𝑖(𝒖)is assumed to be fixed. The latent components 𝒛are identifiable up to permutation and signed scaling under some generally mild conditions on the mixing function 𝒇, the sufficient statistics 𝑻and the auxiliary variable 𝒖. In this study, we construct iVAE assuming Gaussian latent components. Then, for Neural Networks 181 (2025) 106774 3 M. Sipilä et al. identifiability, the variances of the latent components are required to vary enough based on the auxiliary variable 𝒖and the mixing function 𝒇is required to have continuous partial derivatives. The exact identifiability conditions can be found in Khemakhem et al. (2020). The iVAE model is similar to the regular VAE model with the exception that iVAE has an additional auxiliary function 𝒘(𝒖)and its parameters 𝜽𝒘to be estimated, and the encoder 𝒈(𝒙,𝒖)takes both, observations 𝒙and the auxiliary variables 𝒖as an input. The auxiliary function 𝒘maps 𝒖into 𝝁𝒛|𝒖and 𝝈𝒛|𝒖, which are used to calculate the loss. For iVAE model, ELBO is obtained as (𝜽|𝒙,𝒖)≥𝐸𝑞𝜽(𝒛|𝒙,𝒖)(log 𝑝𝜽𝒉(𝒙|𝒛) + log 𝑝𝜽𝒘(𝒛|𝒖) − log 𝑞𝜽𝒈(𝒛|𝒙,𝒖)), where log 𝑝𝜽𝒉(𝒙|𝒛)controls the reconstruction accuracy and log 𝑝𝜽𝒘(𝒛|𝒖) − log 𝑞𝜽𝒈(𝒛|𝒙,𝒖)is a Kullback–Leibler (KL) divergence between 𝑝𝜽𝒘(𝒛|𝒖)and 𝑞𝜽𝒈keeping the distributions 𝑝𝜽𝒘and 𝑞𝜽𝒈as similar as possible. ELBO is maximized to obtain the estimated parameters 𝜽= (𝜽⊤ 𝒈,𝜽⊤ 𝒉,𝜽⊤ 𝒘)⊤. The distributions 𝑝𝜽𝒉,𝑝𝜽𝒘and 𝑞𝜽𝒈are typically Gaussian distributions, where the functions 𝒉,𝒘and 𝒈give the mean and variance parameters of the distributions. The distributions can also be other than Gaussian as long as the resampling can be done using the reparametrization trick to allow the backpropagation go through the resampling node. Then, the functions 𝒉,𝒘and 𝒈do not give mean and variance, but the parameters according the chosen distributions. As we assume Gaussian latent data in this paper, we have 𝑝𝜽𝒘=𝑁(𝒛|𝝁𝒛|𝒖,diag(𝝈𝒛|𝒖)),𝑞𝜽𝒈=𝑁(𝒛|𝝁𝒛|𝒙,𝒖,diag(𝝈𝒛|𝒙,𝒖)) and 𝑝𝜽𝒉=𝑁(𝒙|𝒙′, 𝛽𝑰), where 𝛽 > 0is a small constant as 𝑝𝜽𝒉(𝒙|𝒛)estimates the true distribution 𝑝∗(𝒙|𝒛)with infinitesimal variance. By increasing 𝛽, the weight of the reconstruction accuracy in the ELBO decreases. Based on our empirical investigations, we use 𝛽= 0.02 which provides a good balance between the reconstruction error and KL divergence in ELBO, and leads to good performance. Fig. 1 has representations of both VAE and iVAE models and illustrates the difference between the models. 3. iVAE for STBSS To perform nonlinear spatio-temporal blind source separation using iVAE, the auxiliary variables 𝒖must be selected appropriately. The main assumption for identifiability in spatio-temporal setting is that the variances of the latent components are varying in space and/or in time. This assumption is met by assuming that the latent components are second-order nonhomogeneous, meaning that the second moment of the marginal distribution 𝑝(𝑧𝑖)is not invariant with respect to the location shift in space and/or in time. In addition, the latent components are allowed to be first-order nonhomogeneous, meaning that the components can have nonconstant spatio-temporal trend. The auxiliary variables are constructed in a way that the auxiliary function 𝒘is capable to learn and estimate the mean and the variance vectors of the location of the corresponding multivariate observation. We propose three spatio-temporal iVAE methods; a naive coordinate based method, a segmentation based method extended from spatial iVAE (Sipilä et al.,2024) and a radial basis function based method utilizing ideas of Nag, Sun, and Reich (2023). Each of the three methods construct the auxiliary data differently based on the spatio-temporal location of the observation. Notice that although in many options below the auxiliary variables are constructed separately for spatial coordinates and temporal coordinates, the auxiliary functions can still learn complex spatio-temporal interactions in the mean and in the variance as the auxiliary functions are modelled by deep neural networks. Furthermore, even though the methods for constructing the auxiliary data are defined here for spatial dimension 𝐷= 2, the same ideas apply also for higher 𝐷. The approaches presented here are well scalable in terms of computation time. The computation time grows sublinearly with respect to the sample size 𝑛(as fewer training epochs are typically needed with larger datasets) and linearly with respect to the dimensions of observed data, latent data and auxiliary data. However, if the dimension of the auxiliary data and sample sizes are large, the memory usage may grow unless the auxiliary data are formed batch-wise. An in-depth analysis of computational complexity is provided in Appendix A. 3.1. Coordinate based algorithm In coordinate based iVAE, the preprosessed coordinates are used directly as auxiliary variable. The preprosessed coordinates are obtained by applying min–max normalization to each dimension. The preprosessed coordinates are then 𝑠1=𝑠1−𝑠min 1 𝑠max 1−𝑠min 1 , 𝑠2=𝑠2−𝑠min 2 𝑠max 2−𝑠min 2 and  𝑡=𝑡−𝑡min 𝑡max −𝑡min , where 𝑠min 1, 𝑠min 2and 𝑡min are the minimum coordinates of the locations of the observations, and 𝑠max 1, 𝑠max 2and 𝑡max are the maximum coordinates of the locations of the observations. The algorithm with the preprosessed coordinates, 𝒖(𝒔, 𝑡)=(𝑠1, 𝑠2, 𝑡)⊤, as auxiliary variable is denoted by iVAEc. 3.2. Segmentation based algorithm In segmentation based iVAE, a spatio-temporal segmentation is used as an auxiliary variable. The spatio-temporal segmentation means that the domain ×is divided into 𝑚nonintersecting segments 𝑖∈× so that 𝑖∩𝑗= ∅ for all 𝑖≠𝑗,𝑖, 𝑗 = 1,…, 𝑚, and ∪𝑚 𝑖=1𝑖=×. By using an indicator function 1, the auxiliary variable for the observation 𝒙(𝒔, 𝑡)can be written as 𝒖(𝒔, 𝑡) = (1((𝒔, 𝑡) ∈ 1),…,1((𝒔, 𝑡) ∈ 𝑚)))⊤, where 1((𝒔, 𝑡) ∈ 𝑖) = 1, if the location (𝒔, 𝑡)is within the segment 𝑖, and otherwise 1((𝒔, 𝑡) ∈ 𝑖) = 0. This results into 𝑚-dimensional standard basis vector, where the value 1 gives the spatio-temporal segment in which the location of the observation belongs. If the spatio-temporal domain is large and small segments are used, the dimension 𝑚of the auxiliary variable becomes very large. To lower the dimension, the spatial and temporal segmentations can be considered separately. This means that the auxiliary data is composed of 𝑚𝑆spatial segments 𝑖∈and 𝑚𝑇temporal segments 𝑖∈so that 𝑖∩𝑗= ∅ for all 𝑖≠𝑗,𝑖, 𝑗 = 1,…, 𝑚𝑆,∪𝑚𝑆 𝑖=1𝑖=,𝑖∩𝑗= ∅ for all 𝑖≠𝑗,𝑖, 𝑗 = 1,…, 𝑚𝑇, and ∪𝑚𝑇 𝑖=1𝑖=. Then, the auxiliary variable for the observation 𝒙(𝒔, 𝑡)is 𝒖(𝒔, 𝑡) = (1(𝒔∈1),…,1(𝒔∈𝑚𝑆),1(𝑡∈ 1),…,1(𝑡∈𝑚𝑇)))⊤. The auxiliary variable is (𝑚𝑆+𝑚𝑇)-dimensional and has two nonzero entries for each observation. The dimension can be reduced even further by considering also the 𝑥-axis and 𝑦-axis of the spatial domain separately. Segmentation based auxiliary variables are illustrated in Fig. 2, in which spatial and temporal segmentations are considered separately. We denote the algorithm with all dimensions segmented separately as iVAEs1, with space and time segmented separately as iVAEs2, and with spatio-temporal segmentation as iVAEs3, respectively. 3.3. Radial basis function based algorithm In radial basis function based iVAE, the auxiliary variable is defined using radial basis functions (see e.g. Hastie, Tibshirani, Friedman, and Friedman (2009)). The idea is that with large number of appropriate radial basis functions, the model incorporates much more spatio-temporal information than by using the coordinates only. Similar ideas have been used recently in Chen, Li, Reich, and Sun (2024), Nag et al. (2023) to perform deep learning based spatial and spatio-temporal predicting by using the spatial and spatio-temporal locations transformed into radial basis functions as input for deep neural networks. Following Nag et al. (2023), we transform spatial and temporal locations separately into radial basis functions. Let {𝒐 𝑖},𝑖= 1,…, 𝐾, where 𝒐 𝑖∈, be a set of spatial node points, and let {𝑜 𝑖},𝑖= 1,…, 𝐾, where 𝑜 𝑖∈, be a Neural Networks 181 (2025) 106774 4 M. Sipilä et al. Fig. 2. Illustrations of auxiliary variables of segmentation based iVAE (iVAEs2) (a) and radial basis function based iVAE (b). The top figure of (a) illustrates spatial segmentation, where each segment has size 20 ×20 producing 25 spatial segments, and the bottom figure illustrates temporal segmentation, where each segment has 5 time points, producing 20 temporal segments. In (b), the black lines in the top figure are the normalized 𝑥and 𝑦values at 1∕(𝐻+ 2) = 1∕4 and 1∕(𝐻+ 2) + 1∕𝐻= 3∕4 formed by resolution level 𝐻= 2, and the red points represent the produced spatial node points. The blue points represent temporal node points for resolution level 𝐺= 5. Spatial and temporal Gaussian radial basis functions are illustrated in (c) and (d), respectively. The radial basis functions are functions of distance between node point and a spatial or temporal location. set of temporal node points. The parameter 𝜁is a scale parameter. The spatial and temporal radial basis functions are given as 𝑣(𝒔;𝜁, 𝒐 𝑖) = 𝑣(‖𝒔−𝒐 𝑖‖∕𝜁)and 𝑣(𝑡;𝜁, 𝑜 𝑖) = 𝑣(|𝑡−𝑜 𝑖|∕𝜁), where 𝑣is a kernel function such as the Gaussian kernel 𝑣𝐺(𝑑) = 𝑒−𝑑2, or one of the Wendland kernels (Wendland,1995) such as 𝑣𝑊(𝑑) = {(1 − 𝑑)6(35𝑑2+ 18𝑑+ 3)∕3, 𝑑 ∈ [0,1] 0,otherwise. Following Nychka, Bandyopadhyay, Hammerling, Lindgren, and Sain (2015), we use a multi-resolution approach to form the spatial and temporal radial basis functions. Each resolution level is composed of its own number of evenly spaced node points and own scaling parameter. A low level resolution with small number of node points and large value of the scaling parameter aims to capture large-scale spatial or temporal dependencies, while a high level resolution with many node points and small scaling parameter aims to find finer details of the dependence structure. To form the radial basis functions, we first preprocess the spatial and temporal locations to range [0,1] using min–max normalization. A 𝐻level spatial resolution is formed of evenly spaced grid of node points {𝒐 𝑖}with spacing 1∕𝐻and an offset 1∕(𝐻+ 2) before the first node point, meaning that 𝐻-level resolution has the node points {(𝑖, 𝑗) ∶ 𝑖, 𝑗 ∈ {1 𝐻+2 ,1 𝐻+2 +1 𝐻,…,1 − 1 𝐻+2 }}. For example 2-level spatial resolution is then composed of the node points {(𝑖, 𝑗) ∶ 𝑖, 𝑗 ∈ {0.25,0.75}}. Similarly, a 𝐺-level temporal resolution is formed of evenly spaced one dimensional node points {𝑜 𝑖}with spacing 1∕𝐺and an offset 1∕(𝐺+2), meaning that 𝐺-level temporal resolution is composed of the node points {1 𝐺+2 ,1 𝐺+2 +1 𝐺,…,1 − 1 𝐺+2 }. As the scaling parameters 𝜁𝐻and 𝜁𝐺, for spatial and temporal radial basis functions, we use 𝜁𝐻=1 2.5𝐻 following Nychka et al. (2015) and 𝜁𝐺=|𝑜 1−𝑜 2| √2following Nag et al. (2023). Spatial and temporal node points and radial basis functions are illustrated in Fig. 2 for 𝐻= 2 spatial resolution, producing 4 spatial radial basis functions, and 𝐺= 5, producing 5 temporal radial basis functions. In practice, multiple spatial and temporal resolution levels, such as 𝐻= (𝐻1, 𝐻2) = (2,9), and 𝐺= (𝐺1, 𝐺2, 𝐺3) = (9,17,37), should be used to capture both large scale and finer dependencies. An advantage of using radial basis functions as auxiliary variables instead of spatio-temporal segments is that by using radial basis functions, iVAE’s auxiliary function provides a smooth spatio-temporal trend and variance functions, which can be used later for further analysis such as for prediction purposes. The radial basis function based iVAE is denoted as iVAEr in the rest of the paper. 4. Simulation studies The aim of this section is to demonstrate and compare the performances of iVAE methods using simulation studies and to discover how various types of nonstationarity in variance affect the performance. The section begins with a short review of some common procedures for generating spatio-temporal data and ways to introduce nonstationarity in it. The remainder of the section contains a large simulation study showing the unmixing performances of the iVAE methods under different types of nonstationarity scenarios, and then introduces two methods to estimate the number of latent signals. Finally, the latent dimension estimation methods are illustrated in a small simulation study. All simulations can be reproduced using R 4.3.0 (R Core Team, 2023) together with R packages fastICA (Marchini, Heaton, & Ripley, 2021), SpaceTimeBSS (Muehlmann, Piccolotto, Cappello, De Iaco, & Nordhausen,2022) and NonlinearBSS. NonlinearBSS package contains R implementations of all proposed spatio-temporal iVAE variants, and is available in https://github.com/mikasip/NonlinearBSS. The simulations were executed on the CSC Puhti cluster, a high-performance computing environment. Neural Networks 181 (2025) 106774 5 M. Sipilä et al. 4.1. Nonstationary spatio-temporal data generation Spatio-temporal data are typically composed of 𝑛𝑠spatial locations and 𝑛𝑡temporal points, making the total number of observations 𝑛= 𝑛𝑠𝑛𝑡usually very high. The observations are often collected regularly, for example daily or hourly, by some monitoring stations in different locations. This makes the observed data quickly very dense in time but more sparse in space. To study the properties of the models under the nature of real life spatio-temporal data, generating large datasets with various spatio-temporal covariance models is required. In the following simulations, we exploit a computationally efficient vector autoregressive process, see for example (Papalexiou & Serinaldi,2020;Sigrist, Künsch, & Stahel,2012;Xu & Gardoni,2018;Yan, Huang, & Genton, 2021), and a simplified version of improved latent space approach (ILSA) (Xu & Gardoni,2018) to generate nonstationary spatio-temporal data. Assume spatial field at time 𝑡= 1,…, 𝑛𝑡to be 𝜹(𝑡) = (𝛿(𝒔1, 𝑡),…, 𝛿(𝒔𝑛𝑠, 𝑡))⊤, where 𝒔1,…𝒔𝑛𝑠are the spatial locations in the spatio-temporal field. The vector autoregressive process can be written as 𝜹(𝑡) = 𝑅 ∑ 𝑟=1 𝜌𝑟𝑲𝑟𝜹(𝑡−𝑟) + 𝝐𝜹(𝑡),(6) where 𝑟= 1,…, 𝑅 is an autoregressive order, 𝜌𝑟is 𝑟th baseline autoregressive coefficient, 𝑲𝑟is a 𝑛𝑠×𝑛𝑠spatial kernel matrix determining the change of temporal correlation with spatial locations, and 𝜖𝜹(𝑡)is a𝑛𝑠-dimensional noise vector with covariance 𝐶(𝜖𝜹(𝒔, 𝑡), 𝜖𝜹(𝒔′, 𝑡)) with 𝒔,𝒔′∈ {𝒔1,…,𝒔𝑛𝑠}. With a simplified version of ILSA, one can generate nonstationary spatio-temporal data by using vector autoregressive process (6) as formulated next. Let  𝒔(𝒔) = [𝑠1,…, 𝑠𝑑]be a 𝑑-dimensional transformation of the original coordinate 𝒔, where the transformed coordinates 𝑠𝑖, 𝑖= 1,…, 𝑑, are called regressors or latent coordinates. Let 𝑑𝒔𝑖𝒔𝑗= [‖𝑠1 1−𝑠2 1‖,‖𝑠1 2−𝑠2 2‖]⊤,𝑑 𝒔𝑖 𝒔𝑗= [‖𝑠1 1−𝑠2 1‖,…,‖𝑠1 𝑑−𝑠2 𝑑‖]⊤and 𝑉to be any stationary covariance function. Simplified ILSA has the formulations 𝐾𝑟{𝑖,𝑗}=|||| 𝜽𝒔,𝑟 𝜽 𝒔,𝑟|||| −1 2exp (−[𝒅𝒔𝑖𝒔𝑗𝒅 𝒔𝑖 𝒔𝑗][𝜽𝒔,𝑟 𝜽 𝒔,𝑟][𝒅𝒔𝑖𝒔𝑗 𝒅 𝒔𝑖 𝒔𝑗]), 𝐶(𝜖𝜹(𝒔, 𝑡), 𝜖𝜹(𝒔′, 𝑡)) = 𝜎[ 𝒔(𝒔),𝒔, 𝑡]𝜎[ 𝒔(𝒔′),𝒔′, 𝑡]𝑉(𝑄𝑡)(7) where 𝑄𝑡=([𝒅𝒔𝑖𝒔𝑗𝒅 𝒔𝑖 𝒔𝑗][𝜽𝒔𝜽 𝒔][𝒅𝒔𝑖𝒔𝑗 𝒅 𝒔𝑖 𝒔𝑗])1 2 and 𝜽𝒔,𝑟 =diag(𝜃𝑠1,𝑟, 𝜃𝑠2,𝑟),𝜽 𝒔,𝑟 =diag(𝜃𝑠1,𝑟,…, 𝜃𝑠𝑑,𝑟),𝜽𝒔=diag(𝜃𝑠1, 𝜃𝑠2) and 𝜽 𝒔=diag(𝜃𝑠1,𝑟,…, 𝜃𝑠𝑑,𝑟)are diagonal matrices giving scaling parameters for the spatial coordinates and for the latent coordinates. The function 𝑄𝑡transforms the original coordinates based on the scaling parameters and the latent coordinates. By using this approach, one can easily introduce complex, nonstationary and nonseparable covariance structures through latent coordinates  𝒔, time varying spatial kernel matrices 𝑲𝑟and nonstationary variance function 𝜎. In the following simulations, we are mainly interested in having nonstationarity in the variance as that is required for the identifiability of the latent components. 4.2. Finite sample efficiencies In this section, four different iVAE configurations – regular VAE, symmetric FastICA (FICA) (Hyvärinen,1999) with hyperbolic tangent nonlinearity, and STBSS – are compared using simulated spatiotemporal data. Although FICA is not designed for spatio-temporal data or nonlinear mixing, it is included as a linear baseline for data with nonstationary variances. STBSS, developed for stationary spatio-temporal data and linear mixing, serves as a spatio-temporal baseline. While there are several deep learning-based approaches for nonlinear BSS in the literature, see Ansari et al. (2023) for a recent review, most lack identifiability and focus on acoustic data, which primarily exhibits serial dependence, making them suboptimal for spatio-temporal data. Nonetheless, we include regular VAE as an unidentifiable deep learning baseline. The aim of the simulations is to identify how the proposed iVAE methods perform as compared to other existing methods under various types of nonstationary spatio-temporal data, and how the type of nonstationary affects the performance. To identify how the reduction of either temporal or spatial observations affect the performance of the algorithms, we consider three sample sizes composed of 𝑛𝑠spatial locations and 𝑛𝑡temporal observations for each spatial location. The sample dimensions considered are (𝑛𝑠, 𝑛𝑡) = (150,300),(𝑛𝑠, 𝑛𝑡) = (50,300) and (𝑛𝑠, 𝑛𝑡) = (150,75) yielding 𝑛= 45000,𝑛= 15000 and 𝑛= 11250 observations, respectively. We generate the latent data 𝒛according to six different simulation settings. In some settings the nonstationarity is introduced only in time, in some settings only in space, and in some settings both in space and in time. In each simulation setting, 𝑛𝑠spatial locations 𝒔𝑖are sampled uniformly in spatial domain = [0,1] × [0,1], and for each spatial location 𝒔𝑖,𝑛𝑡observations 𝒙(𝒔𝑖, 𝑡)are generated. The true latent dimension is 𝑃= 5 and the dimension of the observations is 𝑆= 8. Every setting is repeated 500 times for each sample size and for each algorithm. Finally, each trial is repeated using three increasingly nonlinear mixing functions as described hereafter. The first three simulation settings are more simple ones, followed by three more complex ones which utilize the ILSA framework. The simulation settings and the mixing procedure are defined in the following. After introducing the data generation of the settings, the motivation behind each setting is carefully explained. Setting 1. The latent spatio-temporal field consists of three clusters in space and five segments in time yielding 15 spatio-temporal clusters, each of which has their own unique diagonal covariance matrix and unique mean vector. For 𝑘th cluster, 𝑘= 1,…,15, the covariance matrix is given as 𝑪𝑘=diag(𝜎1,𝑘,…, 𝜎5,𝑘), where 𝜎𝑖,𝑘,∼Unif(0.1,5) and unique mean vector is given as 𝝁𝑘= (𝜇1,𝑘,…, 𝜇5,𝑘)⊤, where 𝜇𝑖,𝑘 ∼Unif(−5,5), 𝑖= 1,…,5. Setting 2. The latent spatio-temporal field consists of 10 segments in time. The latent components are simulated by generating first Gaussian spatial data with Matern covariance function using unique parameters (𝜈𝑖, 𝜙𝑖)for each component 𝑖= 1,…,5, and then adding Gaussian iid data with unique covariance matrix and mean vector for each time segment. The Matern parameters are (𝜈1, 𝜙1) = (0.5,0.30), (𝜈2, 𝜙2) = (0.1,0.25),(𝜈3, 𝜙3) = (1,0.35),(𝜈4, 𝜙4) = (2,0.20),(𝜈5, 𝜙5) = (0.25,0.15) and the parameters for the time segment 𝑘= 1,…,10 are 𝝁𝑘= (𝜇1,𝑘,…, 𝜇5,𝑘)⊤, where 𝜇𝑖,𝑘 ∼Unif(−0.3,0.3), and 𝜮𝑘= diag(𝜎1,𝑘,…, 𝜎5,𝑘), where 𝜎𝑖,𝑘,∼Unif(0,0.4),𝑖= 1,…,5. Setting 3. The latent spatio-temporal field consists of five clusters in space and follows AR1 model. In 𝑘th cluster, 𝑘= 1,…,5, the latent components 𝑧𝑖,𝑖= 1,…,5, are generated as 𝑧𝑖(𝒔, 𝑡 + 1) = 𝜌𝑖,𝑘𝑧𝑖(𝒔, 𝑡) + 𝜖𝑖,𝑘,𝑡, 𝜖𝑖,𝑘,𝑡 ∼𝑁(𝜇𝑖,𝑘, 𝜎𝑖,𝑘), where 𝑡= 1,…, 𝑛𝑡− 1 and 𝑧𝑖(𝒔,1) ∼ 𝑁(𝜇𝑖,𝑘, 𝜎𝑖,𝑘). Each cluster has unique parameters 𝝆𝑘= (𝜌1,𝑘,…𝜌5,𝑘)⊤,𝝁𝑘= (𝜇1,𝑘,…𝜇5,𝑘)⊤and 𝝈𝑘= (𝜎1,𝑘,…𝜎5,𝑘)⊤generated as 𝜌𝑖,𝑘 ∼Unif(0.05,0.95),𝜇𝑖,𝑘 ∼Unif(−1,1) and 𝜎𝑖,𝑘 ∼Unif(0.1,5). Settings 4–6. The latent spatio-temporal field is generated using ILSA framework. Each setting has the same highly nonstationary covariance structure. In addition, Setting 4 has a variance 𝜎changing in space, Setting 5 has a variance changing in time, and Setting 6 has a variance changing both in space and in time. The latent coordinates  𝒔= (𝑠1, 𝑠2)⊤are transformed from the spatial coordinates 𝒔= (𝑠1, 𝑠2)⊤by using a swirl-like coordinate transformation according to Papalexiou, Serinaldi, and Porcu (2021) given by 𝑠1= (𝑠1−𝑠∗ 1)cos(𝜂exp (−( ℎ∗ 𝑏swirl )2)) − (𝑠2−𝑠∗ 2)sin(𝜂exp (−( ℎ∗ 𝑏swirl )2))+𝑠∗ 1, Neural Networks 181 (2025) 106774 6 M. Sipilä et al. Table 1 The ILSA parameters and coordinate deformation parameters for Settings 4–6. 𝜽𝒔,𝑙 𝜽 𝒔,𝑙 𝜽𝒔𝜽 𝒔𝜌1𝒔∗𝑏swirl 𝜂 𝜈 𝜙 IC1 (6,4) (7,7) (0.2,0.7) (0.7,0.2) 0.9 (0.5,0.5) 0.7 1.8𝜋0.25 0.5 IC2 (3,6) (4,7) (0.7,0.2) (0.25,0.5) 0.8 (0.7,0.7) 0.4 1.2𝜋0.2 0.9 IC3 (3,3) (6,3) (0.5,0.5) (0.7,0) 0.7 (0.3,0.3) 0.2 2𝜋0.05 1.5 IC4 (7,3) (2,6) (0.2,0.4) (0.3,0.7) 0.6 (0.7,0.3) 10.5𝜋0.1 0.25 IC5 (2,1) (6,2) (0.3,0.3) (0,0.7) 0.5 (0.3,0.7) 0.9 0.9𝜋0.15 1 𝑠2= (𝑠1−𝑠∗ 1)sin(𝜂exp (−( ℎ∗ 𝑏swirl )2)) − (𝑠2−𝑠∗ 2)cos(𝜂exp (−( ℎ∗ 𝑏swirl )2))+𝑠∗ 2, where 𝒔∗= (𝑠∗ 1, 𝑠∗ 2)is the centre point of the deformation, ℎ∗=‖𝒔−𝒔∗‖ is Euclidean distance between the original location and the centre point, 𝜂is a rotation angle, and 𝑏swirl is a scaling parameter controlling the magnitude of the swirl. Each latent component has their own set of deformation parameters. The stationary covariance function 𝑉 in (7) is the Matern covariance function with parameters (𝜈𝑖, 𝜙𝑖)for all Settings 4–6. The deformation parameters, ILSA parameters and Matern parameters mutual for Settings 4–6 are given in Table 1. The autoregressive order is 𝑅= 1 for all settings. In Setting 4, we have 𝜎[ 𝒔(𝒔),𝒔, 𝑡] = exp(𝜃𝑖 𝜎𝑠(𝑠1−0.5)), where 𝜃𝑖 𝜎𝑠is the scaling parameter of variance in space for 𝑖th latent component. This means that the variance of the latent fields vary in space based on the first latent coordinate. The variance scaling parameters for the latent components 𝑧𝑖,𝑖= 1,…,5, are 𝜃1 𝜎𝑠= 1,𝜃2 𝜎𝑠= 2,𝜃3 𝜎𝑠= 3,𝜃4 𝜎𝑠= −1 and 𝜃5 𝜎𝑠= −2. In Setting 5, the variances of the latent fields are changing in time. We set 𝜎[ 𝒔(𝒔),𝒔, 𝑡] = exp(sin((𝑡+𝜃𝑖 𝜎𝑡1) + 𝜃𝑖 𝜎𝑡2)∕2), where 𝜃𝑖 𝜎𝑡1and 𝜃𝑖 𝜎𝑡2 are variance coefficient and variance scaling parameter in time for 𝑖th latent component. The parameters (𝜃𝑖 𝜎𝑡1, 𝜃𝑖 𝜎𝑡2)for the latent components 𝑧𝑖,𝑖= 1,…,5, are (𝜃1 𝜎𝑡1, 𝜃1 𝜎𝑡2) = (50,0.1),(𝜃2 𝜎𝑡1, 𝜃2 𝜎𝑡2) = (0,0.05),(𝜃3 𝜎𝑡1, 𝜃3 𝜎𝑡2) = (100,0.005),(𝜃4 𝜎𝑡1, 𝜃4 𝜎𝑡2) = (20,0.01) and (𝜃5 𝜎𝑡1, 𝜃5 𝜎𝑡2) = (10,0.03). In Setting 6, the variances of the latent fields are changing in space and in time. We set 𝜎[ 𝒔(𝒔),𝒔, 𝑡] = exp(𝜃𝑖 𝜎𝑠(𝑠1−0.5)+sin((𝑡+𝜃𝑖 𝜎𝑡1)+𝜃𝑖 𝜎𝑡2)∕2). The parameters 𝜃𝑖 𝜎𝑠, 𝜃𝑖 𝜎𝑡1, 𝜃𝑖 𝜎𝑡2for the latent fields 𝑧𝑖,𝑖= 1,…,5, are identical as in Settings 4 and Setting 5. Setting 1 has the simplest latent fields by having a diagonal spatiotemporal covariance for each latent field. The variance and mean are changing explicitly between the spatio-temporal clusters as is assumed for segmentation based iVAE. This setting is a spatio-temporal variant of the simulation setting used in time series context in Hyvärinen and Morioka (2016), Khemakhem et al. (2020), where the latent components had multiple temporal segments with unique mean and/or variance parameters. Settings 2 and 3 are still relatively simple with no spatio-temporal interaction in the latent fields. Setting 2 is used to compare performances in cases where latent fields are stationary in space, but variance is changing over time. Setting 3 illustrates a scenario where the latent fields are stationary in time, but the variance is changing over the clusters in space. By having less variability in the variance, Settings 2 and 3 should be less optimal for iVAE. Settings 4– 6 have latent fields with a complex spatio-temporal covariance model and strong spatio-temporal interaction. The variance is not changing explicitly over segments, but instead through a nonstationary covariance function. In Setting 4, the latent fields have smoothly changing nonstationary variance in space, but the variance is stationary in time, and in Setting 5, the variance is nonstationary in time, but stationary in space. Setting 6 introduces nonstationarity in variance both in space and in time. With Settings 4–6 the aim is thus to find out how the presence of nonstationarity in variance affects the performances of iVAE methods in settings with more realistic and more complex spatio-temporal structures. Nonlinear mixing functions. The mixing function 𝑓𝐿is generated using multilayer perceptron (MLP) following Hyvärinen and Morioka (2016), Hyvärinen et al. (2019), Khemakhem et al. (2020). Here 𝐿 denotes the number of mixing layers used in MLP. To obtain an injective and differentiable mixing function, each layer of MLP has 𝑆= 8 hidden units with the activation function 𝜔𝑖being either linear or exponential linear unit (ELU). The matrices 𝑩𝑖,𝑖= 1,…, 𝐿, in the mixing procedure are normalized to have unit length row and column vectors to guarantee that none of the independent components vanish in the mixing process. The mixing function 𝑓𝐿is defined as 𝒇𝐿(𝒛) = {𝜔𝐿(𝑩𝐿𝒛), 𝐿 = 1, 𝜔𝐿(𝑩𝐿𝒇𝐿−1(𝒛)), 𝐿 ∈ {2,3,…}, where 𝑩1is a 8 ×5 matrix and the other matrices, 𝑩𝑖,𝑖≠1, are 8 ×8 matrices. In simulations we use linear activation 𝜔𝐿(𝑥) = 𝑥for the last layer and ELU activation 𝜔𝑖(𝑥) = {𝑥, 𝑥 ≥0, exp(𝑥)−1, 𝑥 < 0, 𝑖= 1,…, 𝐿 − 1, for the other layers. By this procedure, with the number of mixing layers 𝐿= 1, we obtain 𝑆= 8 linear mixtures of the independent components. When the number of mixing layers increase, the mixtures become increasingly nonlinear. In simulations, we consider three different mixing functions with the number of mixing layers 𝐿= 1,3,5. Model specifications. All iVAE models are set up with encoder, decoder and auxiliary function with three hidden layers in each. The hidden layers consist of 128 neurons and leaky rectified linear unit (leaky ReLU) activation (Maas, Hannun, Ng, et al.,2013). iVAEs1, iVAEs2 and iVAEs3 use 4 ×4 spatial segmentation, resulting a grid of 𝑚𝑆= 16 equally sized squares. The temporal segmentation is done by dividing the temporal domain to equally sized segments of length 5. This results the number of temporal segments 𝑚𝑇= 60 when 𝑛𝑡= 300 and 𝑚𝑇= 15 when 𝑛𝑡= 75. For iVAEr, we use spatial resolution levels 𝐻= (2,9), and temporal resolution levels 𝐺= (9,17,37). The iVAE models are trained for 60 epochs when (𝑛𝑠, 𝑛𝑡) = (150,300), for 120 epochs when (𝑛𝑠, 𝑛𝑡) = (50,300) and for 150 epochs when (𝑛𝑠, 𝑛𝑡) = (150,75). The number of epochs is increased when the sample size is decreased, as the number of training steps in each epoch is lower for the smaller sample size. For all sample sizes, the number of epochs are selected large enough to guarantee that the training converges. All iVAE models use learning rate of 0.001 with polynomial decay of secondorder over 10000 training steps, where the learning rate after the first 10000 training steps is 0.0001. VAE uses similar parameters as iVAE, but it does not use any auxiliary data or have an auxiliary function. The STBSS model is fitted with multiple different kernel settings, and the best one is selected, which is having two spatial ring kernels (0,0.15) and (0.15,0.3) and time lag of 1. For more about STBSS and its kernel settings, see Muehlmann et al. (2023). Performance index. To measure the performance of the methods, the mean correlation coefficient (MCC) is used following the previous studies, e.g.,Hälvä and Hyvärinen (2020), Hyvärinen and Morioka (2017), Hyvärinen et al. (2019), Sipilä et al. (2024). MCC is a function of the correlation matrix 𝜴=𝐶𝑜𝑟(𝒛, 𝒛)between the true latent components 𝒛and the estimated ones  𝒛. MCC is calculated as MCC(𝜴) = 1 𝑃sup 𝑷∈tr(𝑷abs(𝜴)),(8) where is a set of all possible 𝑃×𝑃permutation matrices, tr(⋅)is the trace of a matrix and abs(⋅)denotes taking elementwise absolute values of a matrix. MCC gets values in range [0,1], where the optimal value 1 is obtained when the estimated sources are correlated perfectly up to their signs with the true sources. Results. The simulation results are provided in Fig. 3 for (𝑛𝑠, 𝑛𝑡) = (150,300) and in Figs. B.8 and B.9 in the Appendix B for (𝑛𝑠, 𝑛𝑡) = Neural Networks 181 (2025) 106774 7 M. Sipilä et al. Fig. 3. Mean correlation coefficients of 500 trials for Settings 1–6 for sample size with the number of spatial locations 𝑛𝑠= 150 and the number of temporal observations 𝑛𝑡= 300. (50,300) and (𝑛𝑠, 𝑛𝑡) = (150,75), respectively. Based on the results, it is clear that only the iVAE methods are capable of recovering sources through nonlinear unmixing environment. The performances of iVAE methods are better, when the sample size grows, although the differences are small in some settings. The performance of iVAEc is slightly worse than the performances of the other iVAE methods in every setting. FICA performs well in the linear mixing settings, when the latent fields do not contain trend in mean. In nonlinear settings, its performance drops dramatically in all simulation settings. Similar behaviour is present for STBSS, although the performance in zero mean settings does not reach FICA. This is not surprising as STBSS is developed for stationary spatio-temporal random fields. VAE performs poorly in almost all settings, which is expected as the model is not identifiable. In Settings 1–3 all iVAE methods outperform FICA, STBSS and VAE. In these settings, the best performing method is iVAEr, followed by iVAEs1, iVAEs2 and iVAEs3, respectively. They all perform very well under the linear mixing, but when the number of mixing layers is increased, iVAEr outperforms the three other methods. iVAEs1 and iVAEs2 have very similar performance, and they perform slightly better than iVAEs3. The performance of iVAEc is worse than performances of other iVAE methods, especially in Setting 2. FICA performs relatively well in Setting 3 under the linear mixing, but the performance is poor in other settings. STBSS fails to recover the latent fields in Settings 1–3 even under the linear setting. VAE fails in Settings 1 and 2, but performs moderately in Setting 3 under linear mixing. In Settings 4–6, the best method under linear mixing is FICA, but its performance drops considerably in nonlinear settings. In case of nonlinear mixing, the best methods are iVAEs1 and iVAEs2 and iVAEr followed by iVAEs3 and iVAEc, in the order from best to worst. iVAEs1, iVAEs2 and iVAEr perform nearly as well as FICA in linear setting and keep up their good performance also in nonlinear settings. STBSS performs rather well under linear mixing, but is still worse than FICA and iVAE variants. VAE has slightly lower performance than STBSS under linear mixing, but it also fails under nonlinear mixing. In Setting 5, when 𝑛𝑡= 75, the performances of all the methods drop considerably. This is probably due to the fact that in Setting 5, the variance is varying less over the whole temporal domain when the number of time points is lower. The best performing methods, when 𝑛𝑡= 75, are iVAEs1 and iVAEs2 in both linear and nonlinear cases. In Setting 6, the performance of iVAE methods drop only slightly when the number of mixing layers is increased, especially when the sample size is high. The best performing methods for nonlinear mixing environment are iVAEr, iVAEs1 and iVAEs2. Neural Networks 181 (2025) 106774 8 M. Sipilä et al. Fig. B.9. Mean correlation coefficients of 500 trials for settings 1–6 for sample size with the number of spatial locations 𝑛𝑠= 150 and the number of temporal observations 𝑛𝑡= 75. CRediT authorship contribution statement Mika Sipilä: Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Conceptualization. Claudia Cappello: Writing – review & editing, Validation, Investigation, Formal analysis, Data curation. Sandra De Iaco: Writing – review & editing, Validation, Supervision, Investigation. Klaus Nordhausen: Writing – review & editing, Supervision, Resources, Project administration. Sara Taskinen: Writing – review & editing, Supervision, Resources, Project administration. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability Data will be made available on request. Acknowledgements We acknowledge the support from Vilho, Yrjö and Kalle Väisälä foundation for MS, the support from the Research Council of Finland (453691) to ST, the support from the Research Council of Finland (363261) to KN and the support from the HiTEc COST Action (CA21163) to KN and ST. Appendix A. Computational complexity analysis The computational complexities of the proposed algorithms are composed of two parts, forming the auxiliary data and training iVAE. We use Big O notation to represent the worst case time and space complexities, where 𝑂(𝑛)denotes linear growth in computation time or memory usage with respect to the input size 𝑛. First, let us address iVAE’s computational complexity which is very similar to any feed forward neural network such as regular VAE. Using the Big O notation, the computational time complexity for training the model is 𝑂(𝑛×𝑛𝑤×𝑛𝑒), where 𝑛is the sample size, 𝑛𝑤is the number of weights in the model and 𝑛𝑒is the number of epochs. When the sample size 𝑛grows, less epochs are typically needed for training, which makes the model well scaleable Neural Networks 181 (2025) 106774 15 M. Sipilä et al. Fig. B.10. Mean correlation coefficients for different radial basis function parameter settings of iVAEr. The boxplots present 500 trials for Setting 6 with the number of spatial locations 𝑛𝑠= 150 and the number of temporal observations 𝑛𝑡= 300. in terms of sample size. The memory usage, i.e. space complexity, is 𝑂(𝑛𝑤)for storing the weights of the model. The number of weights 𝑛𝑤can be broken down to number of weights 𝑛𝑤1in encoder–decoder part and the number of weights 𝑛𝑤2in auxiliary function. 𝑛𝑤1and 𝑛𝑤2 are heavily dependent on the number and size of the hidden layers. Typically fairly small neural networks (e.g. 3 layers with 128 units) in encoder, decoder and auxiliary functions are sufficient. In addition, 𝑛𝑤1depends linearly on the dimension of the input data and latent dimension whereas 𝑛𝑤2depends linearly on the dimension of auxiliary data. Also, if the input dimension is very large (e.g. more than 100), a larger encoder–decoder network might be needed. Hence, the time complexity grows more when input dimension, latent dimension or auxiliary data dimension grow. In iVAEc, the time and space complexities are the lowest as the coordinates are only scaled to form the auxiliary variables, meaning that using Big O notation, both time and space complexities are 𝑂(𝑛)for forming and storing the auxiliary data. The auxiliary data is only two dimensional, which makes time complexity slightly lower compared to other algorithms. In iVAEs1-iVAEs3, the time complexity of forming the auxiliary data is 𝑂(𝑚1×𝑚2×𝑚×𝑛), where 𝑚1, 𝑚2and 𝑚are the number of segments along each dimension. However, for iVAEs1, where all dimensions are considered jointly, and hence the dimension of auxiliary data can be very large, the space complexity is 𝑂(𝑚1×𝑚2× 𝑚×𝑛). For iVAEs2, the space complexity is 𝑂((𝑚1×𝑚2+𝑚) × 𝑛) and for iVAEs3, it is 𝑂((𝑚1+𝑚2+𝑚) × 𝑛). In terms of computation time, iVAEs3 is the most efficient of segmentation based algorithms as the auxiliary dimension is the lowest. iVAEs2 is also efficient if the number of spatial segments is not very high. In iVAEr, the time and space complexities for forming and storing auxiliary data are 𝑂(𝐾+ 𝐾) × 𝑛, where 𝐾and 𝐾are numbers of spatial and temporal node points, respectively. In all above iVAE variants, the space complexity can be reduced further by constructing auxiliary variables batch-wise during training process. In conclusion, the algorithms are well scalable in terms of sample size 𝑛and relatively well scalable in terms of dimensions of input data, latent data and auxiliary data (linear time complexity). However, if the auxiliary data are not formed batch-wise, memory consumption may grow large if dimension of auxiliary variable is very large. Since the algorithm is essentially composed of three feed forward neural networks, encoder, decoder and auxiliary function, standard parallelization methods such as data parallelism, which distributes the data batch-wise across multiple computation units, can be applied to further reduce the overall computation time. Appendix B. Additional simulation results See Figs. B.8–B.10 References Ansari, S., Alatrany, A. S., Alnajjar, K. A., Khater, T., Mahmoud, S., Al-Jumeily, D., et al. (2023). A survey of artificial intelligence approaches in blind source separation. Neurocomputing,561, Article 126895. http://dx.doi.org/10.1016/j.neucom.2023. 126895. Bachoc, F., Genton, M. G., Nordhausen, K., Ruiz-Gazen, A., & Virta, J. (2020). Spatial blind source separation. Biometrika,107, 627–646. http://dx.doi.org/10. 1093/biomet/asz079. Cappello, C., De Iaco, S., & Palma, M. (2022). Computational advances for spatio-temporal multivariate environmental models. Computational Statistics,37(2), 651–670. Cappello, C., De Iaco, S., & Posa, D. (2020). Covatest: an R package for selecting a class of space-time covariance functions. Journal of Statistical Software,94, 1–42. Chen, W., Genton, M. G., & Sun, Y. (2021). Space-time covariance structures and models. Annual Review of Statistics and Its Application,8, 191–215. Chen, W., Li, Y., Reich, B. J., & Sun, Y. (2024). Deepkriging: Spatially dependent deep neural networks for spatial prediction. Statistica Sinica,34, 291–311. Neural Networks 181 (2025) 106774 16 M. Sipilä et al. Comon, P., & Jutten, C. (2010). Handbook of blind source separation: Independent component analysis and applications. Academic Press, http://dx.doi.org/10.1016/ C2009-0-19334-0. De Iaco, S., Myers, D., Palma, M., & Posa, D. (2013). Using simultaneous diagonalization to identify a space–time linear coregionalization model. Mathematical Geosciences, 45, 69–86. De Iaco, S., Myers, D. E., & Posa, D. (2001). Space–time analysis using a general product–sum model. Statistics & Probability Letters,52(1), 21–28. De Iaco, S., Myers, D. E., & Posa, D. (2002). Nonseparable space–time covariance models: Some parametric families. Mathematical Geology,34, 23–42. De Iaco, S., Myers, D., & Posa, D. (2003). The linear coregionalization model and the product–sum space–time variogram. Mathematical Geology,35, 25–38. De Iaco, S., Palma, M., & Posa, D. (2005). Modeling and prediction of multivariate space–time random fields. Computational Statistics & Data Analysis,48(3), 525–547. De Iaco, S., Palma, M., & Posa, D. (2019). Choosing suitable linear coregionalization models for spatio-temporal data. Stochastic Environmental Research and Risk Assessment,33, 1419–1434. De Iaco, S., & Posa, D. (2012). Predicting spatio-temporal random fields: some computational aspects. Computational Geosciences,41, 12–24. De Iaco, S., & Posa, D. (2013). Positive and negative non-separability for space–time covariance models. Journal of Statistical Planning and Inference,143(2), 378–391. De Iaco, S., Posa, D., Cappello, C., & Maggio, S. (2019). Isotropy, symmetry, separability and strict positive definiteness for covariance functions: a critical review. Spatial Statistics,29, 89–108. Feng, L., Nowak, G., O’Neill, T., & Welsh, A. (2014). CUTOFF: A spatio-temporal imputation method. Journal of Hydrology,519, 3591–3605. Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., et al. (2020). Generative adversarial networks. Communications of the ACM,63(11), 139–144. Hälvä, H., & Hyvärinen, A. (2020). Hidden Markov nonlinear ICA: Unsupervised learning from nonstationary time series. In Conference on Uncertainty in Artificial Intelligence (pp. 939–948). PMLR. Hälvä, H., Le Corff, S., Lehéricy, L., So, J., Zhu, Y., Gassiat, E., et al. (2021). Disentangling identifiable features from noisy data with structured nonlinear ICA. Advances in Neural Information Processing Systems,34, 1624–1633. Hans, W. (2003). Multivariate Geostatistics: An Introduction with Applications. Springer Science & Business Media. Hastie, T., Tibshirani, R., Friedman, J. H., & Friedman, J. H. (2009). vol. 2,The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer. Hyvärinen, A. (1999). Fast and robust fixed-point algorithms for independent component analysis. IEEE Transactions on Neural Networks,10(3), 626–634. Hyvärinen, A., & Morioka, H. (2016). Unsupervised feature extraction by timecontrastive learning and nonlinear ICA. Advances in Neural Information Processing Systems,29. Hyvärinen, A., & Morioka, H. (2017). Nonlinear ICA of temporally dependent stationary sources. In Artificial intelligence and statistics (pp. 460–469). PMLR. Hyvärinen, A., & Pajunen, P. (1999). Nonlinear independent component analysis: Existence and uniqueness results. Neural Networks,12(3), 429–439. Hyvärinen, A., Sasaki, H., & Turner, R. (2019). Nonlinear ICA using auxiliary variables and generalized contrastive learning. In The 22nd international conference on artificial intelligence and statistics (pp. 859–868). PMLR. Khemakhem, I., Kingma, D., Monti, R., & Hyvärinen, A. (2020). Variational autoencoders and nonlinear ICA: A unifying framework. In International Conference on Artificial Intelligence and Statistics (pp. 2207–2217). PMLR. Kingma, D. P., & Welling, M. (2013). Auto-encoding variational Bayes. arXiv preprint arXiv:1312.6114. Kyriakidis, P. C., & Journel, A. G. (1999). Geostatistical space–time models: A review. Mathematical Geology,31, 651–684. Lundberg, S. M., & Lee, S.-I. (2017). A unified approach to interpreting model predictions. Advances in Neural Information Processing Systems,30. Luo, W., & Li, B. (2016). Combining eigenvalues and variation of eigenvectors for order determination. Biometrika,103, 875–887. Luo, W., & Li, B. (2021). On order determination by predictor augmentation. Biometrika, 108, 557–574. Maas, A. L., Hannun, A. Y., Ng, A. Y., et al. (2013). Rectifier nonlinearities improve neural network acoustic models. 30, In Proc. icml (1), (p. 3). Atlanta, GA. Marchini, J. L., Heaton, C., & Ripley, B. D. (2021). FastICA: FastICA algorithms to perform ICA and projection pursuit. URL https://CRAN.R-project.org/package= fastICA. Marcílio, W. E., & Eler, D. M. (2020). From explanations to feature selection: assessing SHAP values as feature selection mechanism. In 2020 33rd SIBGRAPI Conference on Graphics, Patterns and Images (SIBGRAPI) (pp. 340–347). Ieee. Muehlmann, C., Bachoc, F., & Nordhausen, K. (2022). Blind source separation for nonstationary random fields. Spatial Statistics,47, Article 100574. http://dx.doi.org/ 10.1016/j.spasta.2021.100574. Muehlmann, C., Bachoc, F., Nordhausen, K., & Yi, M. (2024). Test of the latent dimension of a spatial blind source separation model. Statistica Sinica,34, 837–865. http://dx.doi.org/10.5705/ss.202021.0326. Muehlmann, C., De Iaco, S., & Nordhausen, K. (2023). Blind recovery of sources for multivariate space-time random fields. Stochastic Environmental Research and Risk Assessment,37, 1593–1613. http://dx.doi.org/10.1007/s00477-022-02348-2. Muehlmann, C., Piccolotto, N., Cappello, C., De Iaco, S., & Nordhausen, K. (2022). SpaceTimeBSS: Blind source separation for multivariate spatio-temporal data. URL https://CRAN.R-project.org/package=SpaceTimeBSS. Nag, P., Sun, Y., & Reich, B. J. (2023). Spatio-temporal DeepKriging for interpolation and probabilistic forecasting. Spatial Statistics,57, Article 100773. http://dx.doi. org/10.1016/j.spasta.2023.100773. Nordhausen, K., Oja, H., Filzmoser, P., & Reimann, C. (2015). Blind source separation for spatial compositional data. Mathematical Geosciences,47(7), 753–770. http: //dx.doi.org/10.1007/s11004-014-9559-5. Nordhausen, K., Taskinen, S., & Virta, J. (2022). Signal dimension estimation in BSS models with serial dependence. In 2022 International Conference on Electrical, Computer, Communications and Mechatronics Engineering (ICECCME) (pp. 1–7). Nychka, D., Bandyopadhyay, S., Hammerling, D., Lindgren, F., & Sain, S. (2015). A multiresolution Gaussian process model for the analysis of large spatial datasets. Journal of Computational and Graphical Statistics,24(2), 579–599. Papalexiou, S. M., & Serinaldi, F. (2020). Random fields simplified: Preserving marginal distributions, correlations, and intermittency, with applications from rainfall to humidity. Water Resources Research,56(2), Article e2019WR026331. Papalexiou, S. M., Serinaldi, F., & Porcu, E. (2021). Advancing space-time simulation of random fields: From storms to cyclones and beyond. Water Resources Research, 57(8), Article e2020WR029466. Porcu, E., Furrer, R., & Nychka, D. (2021). 30 years of space–time covariance functions. WIREs Computational Statistics,13(2), Article e1512. http://dx.doi.org/10.1002/ wics.1512. R Core Team (2023). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing, URL https://www.R-project.org/. Radojičić, U., & Nordhausen, K. (2024). Order determination in second-order source separation models using data augmentation. In J. Ansari, S. Fuchs, W. Trutschnig, M. A. Lubiano, M. A. Gil, P. Grzegorzewski, & O. Hryniewicz (Eds.), Combining, Modelling and Analyzing Imprecision, Randomness and Dependence (pp. 371–379). Cham: Springer. Salvana, M. L. O., & Genton, M. G. (2020). Nonstationary cross-covariance functions for multivariate spatio-temporal random fields. Spatial Statistics,37, Article 100411. Sigrist, F., Künsch, H. R., & Stahel, W. A. (2012). A dynamic nonstationary spatiotemporal model for short term prediction of precipitation. The Annals of Applied Statistics,6(4), 1452–1477. http://dx.doi.org/10.1214/12-AOAS564. Sipilä, M., Muehlmann, C., Nordhausen, K., & Taskinen, S. (2024). Robust second-order stationary spatial blind source separation using generalized sign matrices. Spatial Statistics,59, Article 100803. http://dx.doi.org/10.1016/j.spasta.2023.100803. Sipilä, M., Nordhausen, K., & Taskinen, S. (2024). Nonlinear blind source separation exploiting spatial nonstationarity. Information Sciences,665, Article 120365. http: //dx.doi.org/10.1016/j.ins.2024.120365. Virta, J., & Nordhausen, K. (2017). Blind source separation of tensor-valued time series. Signal Processing,141, 204–216. http://dx.doi.org/10.1016/j.sigpro.2017.06.008. Virta, J., & Nordhausen, K. (2021). Determining the Signal Dimension in Second Order Source Separation. Statistica Sinica,31, 135–156. Wendland, H. (1995). Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Advances in Computational Mathematics,4, 389–396. Xu, H., & Gardoni, P. (2018). Improved latent space approach for modelling non-stationary spatial–temporal random fields. Spatial Statistics,23, 160–181. Xu, R., Vaida, F., & Harrington, D. P. (2009). Using profile likelihood for semiparametric model selection with application to proportional hazards mixed models. Statistica Sinica,19(2), 819. Yan, Y., Huang, H.-C., & Genton, M. G. (2021). Vector autoregressive models with spatially structured coefficients for time series on a spatial grid. Journal of Agricultural, Biological and Environmental Statistics,26(3), 387–408. Yi, M., & Nordhausen, K. (2023). Elasso for estimating the signal dimension in ICA. In 2023 31st European Signal Processing Conference (EUSIPCO) (pp. 2023–2027). http://dx.doi.org/10.23919/EUSIPCO58844.2023.10289956. Neural Networks 181 (2025) 106774 17