extfieldsofstudy:
Computer Science
Mathematics

text:
Deep Markov Spatio-Temporal Factorization

We introduce deep Markov spatio-temporal factorization (DMSTF), a deep generative model for spatio-temporal data. Like other factor analysis methods, DMSTF approximates high-dimensional data by a product between time-dependent weights and spatially dependent factors. These weights and factors are in turn represented in terms of lower-dimensional latent variables that we infer using stochastic variational inference. The innovation in DMSTF is that we parameterize weights in terms of a deep Markovian prior, which is able to characterize nonlinear temporal dynamics. We parameterize the corresponding variational distribution using a bidirectional recurrent network. This results in a flexible family of hierarchical deep generative factor analysis models that can be extended to perform time series clustering, or perform factor analysis in the presence of a control signal. Our experiments, which consider simulated data, fMRI data, and traffic data, demonstrate that DMSTF outperforms related methods in terms of reconstruction accuracy and can perform forecasting in a variety domains with nonlinear temporal transitions.

Introduction
Analysis of large-scale spatio-temporal data is relevant to a wide range of applications in biology, marketing, traffic control, climatology, and neuroscience. Due to the high dimensionality of the data, methods for spatiotemporal analysis either exploit smoothness and correlations in the signal [4,26], or map high-dimensional data onto a lower-dimensional representation [37,44,39,2,6,45,8].
In this paper, we combine these two principles by building on recent approaches for hierarchical probabilistic factor analysis [25,24,35]. These models represent data Y ∈ R T ×D with T temporal and D spatial dimensions as a product Y W F between temporal weights W ∈ R T ×K and spatial factors arXiv:2003.09779v1 [cs.
LG] 22 Mar 2020 F ∈ R K×D . Hierarchical approaches to probabilistic factor analysis explicitly model the correlations between multiple instances {Y 1 , . . . , Y N } by representing W n and F n in terms of some set of latent variables Z n that are drawn from a shared prior. The result is a modeling framework that performs two levels of dimensionality reduction. At the instance level, the model represents highdimensional spatio-temporal data Y n in terms of lower-dimensional matrices W n and F n . At the level of the corpus, we further reduce dimensionality by learning a mapping from latent variables Z n to weights and factors.
The main contribution of this paper is that we develop methods for hierarchical probabilistic factor analysis that uncover common patterns of temporal variation across instances and are amenable to temporal forecasting. To this end, we introduce deep Markov spatio-temporal factorization (DMSTF), a model that learns a deep generative Markovian prior in order to reason about (potentially nonlinear) temporal dynamics forecasting.
DMSTF is related to recent studies that employ temporal smoothness assumptions [9] or more explicitly model dynamics [37,44,39,33,17,2,6] for matrix factorization of multidimensional times series. While most of these methods are not essentially probabilistic (providing only point estimates for imputation/prediction tasks), some Bayesian probabilistic matrix factorization methods have been proposed [34,45,8], and linear temporal dymanics on factor latents, W n , have been adapted to these methods in [43,7,38]. However, these methods do not explicitly adopt a priori assumptions about functional form of spatial factors, F n , when available. This makes them intractable for extremely high dimensional spatial data such as neuroimaging data. Moreover, the linear dynamical assumptions in these methods fall short of modeling potentially nonlinear transitions.
DMSTF builds on topographic factor analysis (TFA) [25] and neural topographic factor analysis (NTFA) [35], which were designed with neuroimaging applications in mind. However, the DMSTF generalizes from both these models, specifically by modeling temporal dynamics, hence, can be applied to spatiotemporal analysis more generally. We demonstrate this in our experiments, which evaluate model performance on simulated data, two large-scale fMRI datasets, and four traffic datasets. Our experiments demonstrate that DMSTF achieves better predictive performance for unseen data relative to existing baselines, both when evaluating the log-likelihood for new instances and when performing next time-point forecasting.
We summarize the contributions of this paper as follows: -DMSTF learns a deep Markovian prior that represents high-dimensional spatio-temporal data in terms of low-dimensional latent variables that can capture nonlinear dynamics of the data. -DMSTF learns a mapping from latent variables to spatial factors and temporal weights, which provide an intermediate representation for downstream regression or classification tasks. -DMSTF is able to perform clustering in the low-dimensional temporal latent space, which can provide visual insights about data.
-DMSTF introduces a control input that can be used for semi-or fullysupervised factorization.
The remainder of this paper is structured as follows. Section 2 presents an overview of related work. Section 3 described the DMSTF model architecture. Section 4 provides details on model training. Section 5 describes experiments on simulated data, neuroimaging datasets and traffic datasets.

Related Works
Factor analysis has been extensively used for reducing dimensionality in spatiotemporal data. Principal component analysis (PCA) [29] and independent component analysis (ICA) [11] are among the most well-known classical factor analysis methods. To accommodate tensor data, and mitigate scalability issues, multilinear versions have been proposed [15,32,10,41]. These methods do not naturally handle missing data. They are also permutation invariant along the batch dimension, and therefore cannot capture temporal dynamics [44]. Spatial factors obtained by these methods are also unstructured and difficult to interpret in many applications. In an early attempt to get temporally smooth structures, Chen and Cichocki [9] developed a non-negative matrix factorization model and applied temporal smoothness and spatial decorrelation regularizers to achieve physiologically meaningful components.
Since then, several matrix/tensor factorization approaches have been proposed for modeling temporal dynamics in multivariate/multidimensional time series. Sun et al. [37] presented a dynamic matrix factorization suited for collaborative filtering setting in recommendation systems using a linear-Gaussian dynamical state space model. Yu et al. [44] proposed to use an autoregressive temporal regularizer for matrix factorization to describe temporal dependencies in multivariate time series. [39] added an additional graph Laplacian regularizer to learn spatial autocorrelations and perform prediction on unknown locations. [33] applied multilinear dynamical systems to the latent core tensor obtained from Tucker decomposition of tensor time series data, similarly to [17]. Bahador et al. [2] enforced a low rank assumption on coefficient tensor of vector autoregressive models, and used a spatial Laplacian regularization for prediction in spatio-temporal data. [6] developed a probabilistic temporal tensor decomposition that models temporal dynamics in latent factor using a multilinear Gaussian distribution with a multilinear transition tensor (using multilinear dynamical system) with additional contextual constraints.
In contrast to the the methods above, which provide point estimates for imputation/prediction tasks, Bayesian probabilistic matrix/tensor factorization methods have been proposed (see [34,45,8]). In these approaches, latent factors have Gaussian priors, and Markov chain Monte Carlo (MCMC) methods are used for approximate inference in training and imputation. However, these models essentially focus on global matrix/tensor factorization without explicitly modeling the local temporal and spatial dependencies between factors. Hence, linear temporal dymanics on factor latents have been adapted to Bayesian Gaussian tensor factorization in [43,7,38]. While some of these methods have been effective in dynamical modeling of multidimensional time series, they do not explicitly adopt any functional form for the spatial domain, which makes them less effective for high dimensional spatial data such as fMRI. Moreover, the linear dynamical assumptions in these methods lack the capacity to characterize complex nonlinear dependencies.
Motivated by recent advances in deep learning, several papers have studied incorporation of neural networks into Gaussian state space models for nonlinear state space modelling [21,42,22,13,18,3]. A common practice is to learn a low-dimensional temporal generative model with an encoding scheme for performing amortized inference. However, these models are not applicable to high dimensional spatio-temporal data without factorization assumptions, and are not tractable in the case of missing data.
A number of fMRI-specific generative models have been proposed in which spatial factors are parameterized by Gaussian kernels to enhance interpretability [24,25,35]. The most closely related of these is NTFA [35], which represents the neuroimaging data for each trial, Y n , in terms of an embedding for the participant, Z p(n) , and an embedding for the stimulus, Z s(n) , which allows the model to capture commonalities in trials associated with a particular participant and stimulus. However, NTFA and other TFA variants assume a prior in which weights are conditionally independent as a function of time, which means that these models do not explicitly encode temporal dynamics. In this sense, DMSTF, which does reason about temporal dynamics, but does not explicitly model the interaction between participants and stimuli, is complementary to NTFA.

Deep Markov Spatio-Temporal Factorization
Model Structure and Variational Inference: DMSTF defines a hierarchical deep generative model for a corpus of instances {Y 1 , . . . Y N }, Here p θ (Z) is a deep generative prior over a low-dimensional set of local (instancelevel) variables Z n . The temporal weights W n are sampled from a Gaussian distribution that is parameterized by neural networks µ w θ and σ w θ . The spatial factors are defined as a deterministic transformation Φ(H n ) of a set of factor parameters H n , which are sampled from a distribution that is paremeterized by neural networks µ f θ and σ f θ . All networks have parameters, which we collectively denote by θ.
We train this model using stochastic variational methods [14,30,20,31]. These methods approximate the posterior p θ (W, H, Z | Y ) using a variational distribu- By maximizing this bound with respect to the parameters θ, we learn a deep generative model that defines a distribution over datasets p θ (Y ), which captures correlations between multiple instances Y n of the training data. By maximizing the bound with respect to the parameters φ, we perform Bayesian inference by approximating the distribution q φ (W, H, Z) p θ (W, H, Z | Y ) over latent variables for each instance.
Model Variants: We will develop three variants of DMSTF, which differ in the set of latent variables Z that they employ. The graphical models for these model variants are shown in Fig. 1. All three models incorporate a deep generative Markovian prior to represent temporal variation in W . To do so, they introduce a set of time-dependent weight embeddings z w t to define the distribution on each row w t of the weight matrix W . Moreover, the K rows of the factor matrix are sampled from a shared prior, which is conditioned on a factor embedding z f . The models differ in the following ways: - Fig. 1(a): This model simplifies the structure described above by assuming that the factor parameters H n and embeddings z f n do not vary across instances, but are instead shared at the corpus-level. Markovian Temporal Latent: We assume a Gaussian distribution for the latent transition probability. Given the latent z w t−1 we parameterize the mean and covariance of the diagonal Gaussian distribution p θ (z w t |z w t−1 ) using a neural network. In the model in Fig. 1(c) this distribution p θ (z w t |z w t−1 , u t−1 ) is additionally conditioned on a control variable u t−1 . Concretely, we pass z w t−1 (and u t−1 when applicable) to a multilayer perceptron (MLP) for estimating the Gaussian parameters. We combine a linear transformation of z w t−1 with the estimated mean from the neural network to support both linear and nonlinear dynamics [22]: where L θ (·) is a linear mapping, F θ (·) is the nonlinear mapping of neural network, and g ∈ [0, 1] is a scaling vector which itself is estimated from z w t−1 using a neural network. We show the neural network architectures we employ in Table 1.
Clustering Latent: In the models in Fig. 1(b) and Fig. 1(c), we assume that each sequence Y n belongs to a specific state out of S possible states, and is determined by the categorical variable c n in our temporal generative model. This is sampled from a categorical distribution c n ∼ Cat(π), where π = [π 1 , · · · , π S ] specifies cluster assignment probabilities. To this end, we assume that the the first temporal latent z w n,0 is distributed according to a Gaussian mixture: where the cluster assignment c n enforces µ s and diagonal covariance Σ s .

Temporal & Spatial Factors:
As with the transition model, we assume Gaussian distributions for temporal weights, and spatial factors. We parameterize the diagonal Gaussian distributions for temporal weights p θ (w t |z w t ) and factor parameters p θ (H | z f ) with neural networks. z f itself is sampled from a normal distribution: z f ∼ N (0, I). Introducing z f as a low dimensional spatial embedding in the model encourages estimation of a multi-modal distribution among spatial factors.
The form of the spatial factor parameters H depends on the task at hand. In the case of fMRI data, we employ the construction used in TFA [25,24,35], which represents each factor as a radial basis function with parameters H k = {ρ k , γ k }: This parameterization represents each factor as a Gaussian"blob". The vector r d ∈ R 3 denotes the position of voxel with index d. The parameter ρ k ∈ R 3 denotes the center of the Gaussian kernel, whereas γ k ∈ R controls its width.
In the case of traffic forecasting experiments, we learned spatial factors without any functional form constraints, hence, H directly parameterized {F kd } K,D k=1,d=1 by mean and covariance of a Gaussian distribution (i.e., Φ(·) is basically an identity mapping here).
Variational Distribution: We introduce trainable variational parameters (mean and sigma of a diagonal Gaussian) for each data point in our dataset to define a fully-factorized (i.e. mean-field) variational distributions on the variables We use a structured variational distribution for the variables q φ (z w 1:T | w 1:T ) in the form of a one-layer bidirectional recurrent neural network (BRNN) with ReLU activation, which is then combined with z w t−1 through another neural network to form distribution parameters of z w t : where f φ is a nonlinear mapping parameterized by an MLP, and Σ is diagonal. Although we can define variational parameters for categorical distributions q(c n ), we approximate it with posterior p(c n |z w n,0 ), as suggested in [16], to compensate information loss induced by the mean-field approximation: .
Kullback Leibler Divergence Terms: A significant advantage of our formulation is that we can analytically calculate the KL terms of evidence lower bound for both multivariate Gaussian, and categorical distributions, which leads to lower variance gradient estimates, and faster training as compared to e.g., noisy Monte Carlo estimates often used in literature.

Training Details
We describe the network architectures for the neural networks used in DMSTF in Table 1, where D z is the dimension of z w and z f , and D t and D e are the dimensions of hidden layers for Markovian and temporal latents respectively. For all the experiments in this paper, we assume D z = 2 (except for traffic dataset). We did all the programming in PyTorch v1.3 [28], and used the Adam optimizer [19] with learning rate of 1 × 10 −2 . We initialized all the parameters randomly except for spatial locations of fMRI data which we set the initial values to the local extrema in their averaged fMRI data. We clipped spatial locations and scales to the confines of the brain if needed. We used KL annealing [5] to suppress KL divergence terms in early stages of training, since these terms could be quite strong in the beginning, and we do not want them to dominate the log likelihood term (which controls reconstruction) in early stages. We used a linear annealing schedule to increase from 0.01 to 1 over the course of 100 epochs. The number of learnable parameters in our model is dominated by w t parameters, and therefore will be O(KNT). We run and learn all of the models on an Intel Core i7 CPU @3.7 GHz with 8 Gigabytes of RAM, which proves tractability of the learning process. Per-epoch training time varied from 0.1 in small datasets to 6.0 minutes in larger experiments.
Test/Prediction: We evaluate the performance of DMSTF for short-time prediction tasks (one-point-ahead forecasting) by adopting a rolling prediction scheme as in [8,44]. We predict observations for the next time point using the generative model and spatial factors:Ŷ t+1 =ŵ t+1 F , whereŵ t+1 ∼ p(ŵ t+1 |ẑ t+1 ), andẑ t+1 ∼ p(ẑ t+1 |z t ). Then, we fix the parameters of the generative model and variational parameters for spatial factors, and update variational parameters for z t+1 , w t+1 using Y t+1 , the actual observation at t + 1 (running model for a few epochs). We repeat these steps to make predictions in a rolling manner across a test dataset.

Toy Example
We generated N = 100 synthetic spatio-temporal data with T = 15, K = 2 using a nonlinear dynamical model (motivated by [22]): where ρ = 0.2, α = 0.5, β = −0.1. For spatial factors, F , we picked two Gaussian blobs centered at ±(7.5, 7.5, 7.5) with scales of 3, 4.5 respectively in a box of 30 × 30 × 30 at origin. And finally generated {Y n = W n F } N n=1 with additive noise. We train DMSTF with this synthetic dataset, and estimate the parameters of the model given the functional forms of the generative model. As depicted in Fig. 2 (a) DMSTF is able to recover the actual values of the parameters.

Synthetic Data
We generate synthetic fMRI data using a MATLAB package provided by [25], which is known to be useful for analysing fMRI models. The synthesized brain image for each trial (time point) is a weighted summation of a number of radial basis functions (spatial factors) randomly located in the brain. The synthesized fMRI data is then convolved with a Hemodynamic Response Function (HRF), and finally we add zero-mean Gaussian noise with a medium-level signal-tonoise ratio. Here, we considered 30 activation sources (spatial factors) randomly located in a standard MNI-152-3mm brain template with roughly 270, 000 voxels, and 150 trials. We randomly split these 30 activation sources into 3 groups, each having 10 of the Gaussian blobs. These three groups of sources are periodically activated in turn (according to some random weights) for 5 trials. We generated non-overlapping sequences of T = 5 time points from this synthetic fMRI data. This resulted in 10 data points for each activation group (N = 30). In order to train DMSTF, we set T = 5, K = 30, D t = 2, D e = 8, S = 3, and σ 0 = 1 × 10 −2 . As depicted in Fig. 2 (b), our model was able to successfully recover the 3 clusters of activation that were present in the dataset in a totally unsupervised manner. Predictions of learned generative model are visualized in Fig. 2 (d), which proves that we were able to capture nonlinear hemodynamic response function by using multi-layer perceptrons in our temporal generative model.

Neuroimaging Datasets
Autism Dataset: We used the publicly available preprocessed resting state fMRI (rs-fMRI) data from the Autism Brain Imaging Data Exchange (ABIDE) collected at 16 international imaging sites [12]. This dataset includes rs-fMRI imaging from 408 individuals suffering from Autism Spectrum Disorder (ASD), and 476 typical controls. Each scan has T = 145 ∼ 315 time points at TR = 2, and D = 271, 633 voxels. We split the signals into sequences of 75 time points. We take two approaches to evaluate the performance of our model in separating ASD from control: (1) Cluster data directly in the low dimensional latent, z w , using the clustering feature of DMSTF (called DMSTF), (2) Extract functional connectivity matrices from learned weights, W , followed by a 10fold SVM for classification (called DMSTF+SVM). As a baseline, we performed a 10-fold SVM classification on extracted connectivity matrices from averaged signals of 116 regions of interest (ROIs) in AAL atlas [40] (called ROI+SVM). Several studies have been done on this dataset to differentiate ASD group from control, all of them using supervised methods, and could achieve accuracies up to 69% (by carefully splitting data to be as homogeneous as possible, and reducing site-related variability) using the signals extracted from anatomically labeled regions in the brain [1,27,36]. We set T = 75, K = 100, D e = 15, D t = 5, S = 2, σ 0 = 1 × 10 −2 , and trained DMSTF for 200 epochs on the entire dataset (Full), and also datasets from 9 sites (with more balanced datasets) separately: Caltec, Leuven, MaxMun, NYU, SBL, Stanford, UM, USM, Yale. As shown in Fig. 3 (c) DMSTF and DMSTF+SVM outperform ROI+SVM in Caltec, MaxMun, SBL, Stanford, and Yale, while ROI+SVM only performs better in NYU dataset. DMSTF+SVM performs slightly better than ROI+SVM on the entire dataset (Please note that DMSTF is a clustering approach, hence, no error bars are provided in Fig. 3 (c)). Clustering results for Caltec, Maxmun, SBL, and Stanford are shown in Fig. 3 (a) in which ASD, and control seems to be partially separable (see Fig. S1 in Appendix B for more visualization results).
Depression Dataset: In this dataset [23], 19 individuals with major depressive disorder (MDD) and 20 never-depressed (ND) control participants listened to standardized positive and negative emotional musical and nonmusical stimuli during fMRI scanning. Each participant underwent 3 musical, and 2 nonmusical runs each for 105 time points at TR=3 with D = 353, 600 voxels. During each run, each stimulus type (positive, and negative) was presented for 33 seconds (∼ 11 time points) interleaved with instances of neutral tone of the same length. We discard instances of neutral tone, and split each run into non-overlapping sequences of T = 6 time points in agreement with stimuli design (each stimuli block is split into two sequences). In other words, each run has 4 sequences associated with "positive stimuli", and 4 with "negative stimuli" resulting in a total of 8 data points for each run. In the first experiment, we trained DMSTF on the entire musical runs (N = 39 × 3 × 8 = 936) by setting T = 6, K = 100, D e = 15, D t = 5, σ 0 = 1×10 −3 for 200 epochs. The results are shown in Fig. 4 (a,  Left). We observed that DMSTF fully separated data points associated with each subject into distinct clusters across the low-dimensional temporal latent space. In other words, DMSTF was able to re-unite pieces of signals associated with each subject without any kind of supervision. More importantly, DMSTF was able to partially separate data points associated with MDD group from control. As seen in Fig. 4 (a, Left), MDD group data points are fairly populated in the center of temporal latent while control group are dispersed across latent space. However, DMSTF was not able to meaningfully separate "negative" and "positive" music pieces in latent embedding from a subject-level perspective, since the variation between runs of a subject dominates stimulus-level variation. For this reason, in a second experiment, we focused on 5 subjects, and their first musical run from both MDD and control group and trained DMSTF respectively. Again, as expected, data points from each subject were distinctly clustered in latent space (see middle and right columns in Fig. 4 (a)). Additionally, DMSTF was able to fit two partially separating Gaussians to "positive", and "negative" stimuli per subject. However, since the number of data points for each subject and run is Table 2. Held-out Log-Likelihood. DMSTF results in models with higher held-out likelihood, and therefore better fit comparing to HTFA.

Dataset
HTFA DMSTF Autism (Caltech) −2.82 × 10 6 −2.33 × 10 6 Depression −6.64 × 10 5 −5.71 × 10 5 limited it is not clear how significant these clusters are. A dataset with longer runs could possibly answer that. In a third experiment, we incorporate control inputs u t , and evaluate prediction performance of DMSTF. We train DMSTF on 2 musical and 1 nonmusical runs for a subject with depression (8 × 3 = 24 sequences) using u t = 1, and u t = 0 respectively (i.e., u t is encoding musical, nonmusical stimuli). We predict the remaining musical and nonmusical runs (8 × 2 = 16 sequences) once with u t = 0, and another time with u t = 1.
As reported in Fig. 4 (d), nonmusical sequences are better predictable with u t = 0 than u t = 1 with p-value= 0.011 (vice versa for musical sequences with p-value=0.036). Sample predicted fMRI time series from both musical, and nonmusical runs are shown in Fig. 4 (c) (see Fig. S2 in Appendix B for more visualization results).
Comparison with HTFA [24]: We further evaluated DMSTF against HTFA, an established probabilistic model for multi-subject fMRI analysis, which uses unimodal Gaussian priors for both temporal weights, and spatial factor parameters, in terms of held-out log-likelihood (see Table 2). For autism, we used Caltec site dataset, and split each subject's fMRI into two half (each T = 70). We trained DMSTF on the first half, and test on the second half. For depression, we consider 4 sequences from each subject's run for training, and test on the remaining 4 sequences. To this end, after training DMSTF on each training set, we fixed all parameters of the generative model and only tuned variational parameters of the test set for temporal latents z w t , w t . And finally computed an importance sampling-based estimate of the log-likelihood [22]. The results are shown in Table 2, which proves that DMSTF results in models with higher likelihood on test set, hence a better fit as compared to HTFA.  . We picked the data from January 1 to January 28 (28 days) as in [8], and organized it into a tensor of 28 × 288 × 323. We compared DMSTF (variant (a)) with three state-of-the-art baselines on our short-time prediction task: TRMF [44], BayesTRMF developed in [8], and BTMF [8]. We picked the last seven days from the Birmingham dataset, and the last five days from the Guangzhou, Hangzhou, and Seattle datasets for prediction, then trained the models on the rest for each dataset with K=10, 30, 10, 30 respectively (consistent setup with [8]). For DMSTF, we additionally set {D z , D t , D e } = 5, σ 0 = 0 for all datasets, and learned spatial factors without any functional form constraints. We trained DMSTF for 500 epochs. We report root mean squared error (RMSE), and mean absolute percentage error (MAPE) for all models on the testsets in Table 3. DMSTF outperforms in short-time prediction on all datasets, doing significantly better particularly on the Birmingham dataset. Testset predictions for two sample locations from each dataset are shown in Fig. 5 (see Fig. S3 in Appendix B for a number of long-term prediction visualizations).

Conclusion and Future Work
We presented deep Markov spatio-temporal factorization, a new probabilistic model for robust factor analysis of high-dimensional spatio-temporal data. We employed a chain of low-dimensional Markovian latent variables connected by deep neural networks as a state-space embedding for temporal factors in order to model nonlinear dynamics in data, account better for noise and uncertainty, and enable generative prediction. We also employed a low-dimensional spatial embedding to generate a multi-modal distribution of spatial factors. We then demonstrated the tractability of DMSTF on fMRI data (with very high spatial dimensionality) by incorporating functional form assumptions, and on traffic data with high temporal dimensionality. DMSTF enables clustering in the lowdimensional temporal latent space to reveal structure in data (e.g., cognitive states in fMRI), providing informative visualizations about the data. We plan to extend DMSTF to accommodate higher order dynamics for long-term prediction tasks. We can readily achieve this in our setting by conditioning temporal latents on a time-lag set, such as conditioning z w t on z w t−1 , z w t−2 in a 2-Markov chain. A Derivation of Evidence Lower BOund (ELBO) for DMSTF We derive the lower bound on the log-likelihood of observations by incorporating conditional independences inferred from the graphical model of DMSTF in Fig. 1 (c) as follows: Y n ⊥ ⊥ Y ¬n |w n,1:T , H n , z w n,t ⊥ ⊥ z w n,¬(t,t−1) |z w n,t−1 , z w n,t ⊥ ⊥ z w ¬n |z w n,t−1 z w n,0 ⊥ ⊥ z w ¬n,0 |c n , w n,t ⊥ ⊥ w n,¬t |z w n,t , w n,t ⊥ ⊥ w ¬n |z w n,t , H n ⊥ ⊥ H ¬n |z f n Considering these conditional independencies, the joint distribution of observations and latent variables will be (denoting Z = {z w , z f }): We assume a fully factorized variational distribution on {c, z f , z w 0 , w}, and a structured variational distribution on q φ (z w 1:T |w 1:T ), hence: We then derive the ELBO by writing down the log-likelihood of observations, and plugging in p(·) and q(·) from Equation (1) and Equation (2) p(Y n |w n,1:T , H n )p(H n |z f n )p(z f n )p(z w n,0 |c n )p(c n ) q(H n )q(z f n )q(z w n,0 )q(c n ) · · · T t=1 p(w n,t |z w n,t )p(z w n,t |z w n,t−1 , u n,t−1 ) q(w n,t )q(z w n,t |z w n,t−1 , w n,1:T ) · · · z ={w n,1:T ,Hn} c =∅ q(w n,1:T , H n ) log p(Y n |w n,1: q(w n,1:T )q(z w n,t−1 |w n,1:T ) · · · q(z w n,t |z w n,t−1 , w n,1:T ) log p(z w n,t |z w n,t−1 , u n,t−1 ) q(z w n,t |z w n,t−1 , w n,1:T ) dz + T t=1 z ={z w n,t ,wn,t} c =∅ q(z w n,t )q(w n,t ) log p(w n,t |z w n,t ) q(w n,t ) dz where the integration over distributions of latent variables outside of the summation over n is c\c q(c\c ) Z\z q(Z\z )dZ\z = 1 for {c , z } of each term inside the summation (we are abusing notation here, for the sake of brevity). Considering that KL(q, p) = q log q p , and E q(x) [f (x)] = x q(x)f (x), we can rewrite each term of the summation in Equation (3) to summarize the ELBO: where, L rec n = E q φ (w n,1:T ,Hn) [log p θ (Y n |w n,1:T , H n )] L C n = −KL(q φ (c n )||p θ (c)) − cn q φ (c n )KL(q φ (z w n,0 )||p θ (z w n,0 |c n )) L z w n,t = −E q φ (w n,1:T ) E q φ (z w n,t−1 |w n,1:T ) KL q φ (z w n,t |z w n,t−1 , wn,1:T ) p θ z w n,t |z w n,t−1 , un,t−1) L W n,t = −E q φ (z w n,t ) [KL(q φ (w n,t )||p θ (w n,t |z w n,t ))] Kulback Leibler (KL) Divergence Terms: We can analytically calculate the KL terms of ELBO in Equation (4). For two multivariate (d-dimensional) Gaussian distributions q(·), and p(·), the KL divergence is: which would be further simplified here, in the case of diagonal covariances. The KL divergence between two S-dimensional categorical distributions q(·) and p(·) can be readily computed: where the subscript s denotes the probability of state s out of S possible states.

B More Visualizations for Autism, Depression, and Traffic Datasets
We have visualized real and reconstructed brain images from the nine subsets of autism dataset (Caltex, Leuven, MaxMun, NYU, SBL, Stanford, UM, USM, and Yale sites) along with z w 0 after training DMSTF on the full autism dataset in Fig. S1. DMSTF clustered together temporal latent variables associated with each acquisition site without supervision in z w 0 . As depicted, the variation among different acquisition sites dominates the cognitive differences between ASD group and control, hence, a downstream connectivity matrix classification (using learned temporal weights, W ) helps better in differentiating ASD group from control in multi-site analysis. In Fig. S2, we have visualized example predicted fMRI time-series from both musical and non-musical trials in the test-set of depression dataset using control variable u t = 1 for musical and u t = 0 for non-musical trials. In Fig. S3, we have visualized next-day (long-term) prediction results for Birmingham and Huangzhou subsets of traffic data for four sample locations. These predictions are purely obtained from the trained temporal generative model. Please note that the actual values for the predicted days are not available. showing the smoothing given by sparse factorization. Visualizing z w 0 after training DMSTF on the full autism dataset. DMSTF clustered together temporal latent variables associated with each acquisition site without supervision. As depicted, the variation among different acquisition sites dominates the variation in cognitive state of the brain (ASD group vs. control), hence, a downstream connectivity matrix classification helps better in differentiating ASD group from control in multi-site analysis.   S3. Visualizing next-day (long-term) prediction results for Birmingham and Huangzhou subsets of traffic data for four sample locations. Although DMSTF is wellsuited for short-time prediction, next-day forecasts (purely from the trained temporal generative model) show its capability in long-term predictions. Please note that the actual values for the predicted days are not available.


text:
Iterative Method for Tuning Complex Simulation Code

Tuning a complex simulation code refers to the process of improving the agreement of a code calculation with respect to a set of experimental data by adjusting parameters implemented in the code. This process belongs to the class of inverse problems or model calibration. For this problem, the approximated nonlinear least squares (ANLS) method based on a Gaussian process (GP) metamodel has been employed by some researchers. A potential drawback of the ANLS method is that the metamodel is built only once and not updated thereafter. To address this difficulty, we propose an iterative algorithm in this study. In the proposed algorithm, the parameters of the simulation code and GP metamodel are alternatively re-estimated and updated by maximum likelihood estimation and the ANLS method. This algorithm uses both computer and experimental data repeatedly until convergence. A study using toy-models including inexact computer code with bias terms reveals that the proposed algorithm performs better than the ANLS method and the conditional-likelihood-based approach. Finally, an application to a nuclear fusion simulation code is illustrated.

Introduction
Modern computer simulation codes contain various unknown parameters. Assuming the validity of the simulation code, we can adjust or estimate such parameters using the nonlinear least squares estimation (NLSE) method, which minimizes the sum of the squared differences between computer responses and real observations. This procedure is called calibration (e.g., Kennedy and O'Hagan 2001; Higdon et al. 2008; Tuo and Wu 2018) or code tuning (e.g., Cox, Park, and Singer 2001;Kumar 2015). It is formally defined as the process of improving the agreement of a code calculation or set of code calculations with respect to a chosen and fixed set of experimental data via adjustment of the parameters implemented in the code (Trucano et al. 2006). Han, Santner, and Rawlinson (2009) differentiated between tuning parameter and calibration parameter. In this study, however, the two parameters are treated as the same, henceforth it is referred to as the tuning parameter.
If a simulation program is complex, with one execution requiring several hours, the NLSE method may not be computationally feasible. In this case, a statistical metamodel can be built to approximate the unknown functional relationship between a set of controllable input variables and a simulated response variable. This metamodel is then employed in place of the original simulation code in the NLSE method, making the problem solvable and computationally feasible. This method was described by Cox, Park, and Singer (2001), where a Gaussian process (GP) model was employed as the metamodel of a complex simulation code. It is called the approximated NLS (ANLS) method. As alternatives to the ANLS, Cox, Park, and Singer (2001) proposed likelihood-based methods. Kennedy  A potential drawback of the ANLS method is that the metamodel is built only once and not updated thereafter; thus, the computer data are no longer used live. To address this, in the present report, an iterative algorithm is proposed, in which the likelihood function is maximized, and the squared distance is minimized iteratively until convergence is achieved. That is, the tuning parameters of the simulation code and model parameters of the GP are repeatedly reestimated and updated. The parameters of the GP model are estimated from updated combined data using the maximum likelihood approach. Then, the tuning parameters are re-estimated using the ANLS method. These two optimizations are iteratively executed until convergence is achieved. We call this method the "Max-min" algorithm.
The remainder of this report proceeds as follows. Section 2 presents a GP model to approximate a complex simulation code and real experimental data. Section 3 describes the calibration methods, including our proposed method, and the calculation of standard errors. Section 4 presents a toy-model simulation study. Section 5 applies the proposed methods to a nuclear fusion simulator, and Section 6 summarizes the results and provides suggestions for future research. Some details of this study are provided in the Supplemental Material.

Gaussian process model as a surrogate
For the metamodel of the simulation code, we use a GP or a spatial regression model that treats response y(x) as a realization of a random function superimposed on a regression model: where f s are known functions, and βs are unknown regression coefficients. Here, the random process Z(.), which represents a departure from the assumed linear model, is assumed to be a GP with mean zero and covariance cov(t, u) = V (t, u) = σ 2 R(t, u) between Z(t) and Z(u) for t = (t 1 , ..., t d ), u = (u 1 , ..., u d ), where σ 2 represents the process variance (a scale factor), and R(t, u) is the correlation function. When the response of a computer code is stochastic, the random component term is added to the model (1). However, we do not include in this study because a computer code is assumed to be deterministic. Some possible choices for the correlation function are obtained from the Gaussian correlations denoted by where θ ≥ 0. These are special cases of a power exponential family with a power of 2. The non-negative parameter θ determines the covariance structure of Z: a small θ reflects high correlations between nearby observations, whereas a large θ reflects low nearby correlations. One may consider a different version of (3) by taking several θ values as follows: where θ i ≥ 0, for i = 1, 2, · · · , d. This is a legitimate separable correlation function because it is a product of valid correlation functions. We refer the readers to Santner, Williams, and Notz (2018) for more information on this GP model and its application to the design and analysis of computer experiments. Once the data have been collected at the design sites or at "training" inputs X, the parameters are estimated via the maximum likelihood estimation (MLE) method and are then plugged in to predict y(x 0 ) as in (5), where x 0 is an untried site or a "test" input. This prediction is called Kriging. The empirical best linear unbiased prediction  with the MLEs of the parameters is denoted bŷ where f 0 is the known linear regression function vector, F is a design matrix, r t 0 is the correlation vector between Y (x 0 ) and model outputs Y (X), y is the vector of observations collected at the design sites, andβ is the generalized least squares estimator of β (see the Supplemental Material for the details).
The combinations of βs and θs determine the model, but the following simple GP model is considered first: with the correlation function (3) or (4). When a common θ of (3) is used, we call (6) "Model 1" in this study. When the several θ i s of (4) are used, we call (6) "Model 2".
3 Methods for code tuning

Data structure
For notational convenience, experimental data is denoted by the subscript "E," computer simulation data by the subscript "C," and "both" computer and experimental data by "B." Let τ be an adjustable parameter vector to be estimated. Let T be the input variables of the computer code corresponding to τ . The original experimental input variables are denoted by X. Let q and p be the dimensions of τ and X. Further, let n C andn E be the number of observations; then, n B = n C + n E . The details of the data structure for calibration are described in the Supplemental Material.

Approximate nonlinear least squares
In this subsection, the GP model approach to be used in the ANLS procedure provided in Cox, Park, and Singer (2001) is described. If the parameters σ 2 , β, and θ are known, then, for a given value of τ , a "prediction" of y E (τ, x iE ) can be calculated for the given computer data. This is obtained using equation (5) and X E , F E , y E , and V EE . Here, the computer data alone are used to calculate β andθ, and the data are not used thereafter. Note that X E , F E , R(X E , X C ), R(X E , X E ), and β are functions of τ . A design site selected for a computer experiment is denoted by (T , x C ). Then, the computer response y (or y C ) at (T , x C ) is where Y represents the expected value of the output from computer code. Y can be viewed as the value produced from a theory. Since we assume that the computer code is close to the real experimental data with some variation if the tuning parameters in the computer code are optimal, the response of the real experiment y E at (τ , x E ) is modeled by Here Y (τ , x E ) also represents the expected value of the response in the real experiment. This common Y in the two abovementioned equations (7) and (8) connects the computer code and the real experiment. The stochastic term E is assumed to be independent and identically distributed with mean zero and variance σ 2 E . On the other hand, a computer code is sometimes considered inexact because of computer model bias, which is the systematic difference between the model and the truth (Kennedy and O'Hagan 2001; Gramacy 2016; Plumlee 2017). In that case, the response of the real experiment can be modeled by where b(τ , x E ) represents the model bias or discrepancy. A drawback of the computer experiment is that one run of a complex simulator sometimes requires several minutes, say, m minutes. Then, τ is usually estimated by minimizing the residual sum of the squares: where y E i is an observed response from the real experiments, and Y (τ , x E i ) is the expected value of the output from the computer code at the experimental point (τ , x E i ). One evaluation of RSS(τ ) requires approximately n E × m minutes. It is thus computationally infeasible to run the code as many times as needed for an iterative nonlinear optimizer to find τ . The ANLS method first fits the GP Model 1 or Model 2 as defined in equation (6), with MLE using computer data alone. Then, upon treating the fitted prediction model as if it were the true model, the method attempts to determine the τ that minimizes the residual sum of squares with predictors: whereŶ (τ , x E i ) is the empirical best linear unbiased prediction of Y (τ , x E i ), as in (5). Since it is difficult to have a closed-form minimizer of (11), a numerical optimization routine is necessary to determine τ . Note thatŶ is a computationally cheaper emulator (surrogate or metamodel) of the expensive simulation code. This makes the problem computationally feasible. It is applied for each GP model (Model 1 and Model 2) stated in (6). The advantages of this method are that it is reasonably economical and easy to implement, and the computer and experimental data are uncoupled. The ANLS method does not require the functional relation between the inputs and output to be the same for the computer code and real experiments. In contrast, the likelihood-based approaches described in the next subsection do require this, but they employ marginal likelihoods to estimate model parameters. The prediction residuals, r i = y E i −Ŷ ( τ , x E i ), can be used to check the validity of prediction model, τ , and of the ANLS method. A potential drawback of the ANLS method is that it does not account for uncertainty in the approximation of Y C byŶ . Another difficulty is that the metamodelŶ is built only once and is not updated thereafter.

Likelihood-based tuning methods
Given a computer code and our GP approach for y, a unified statistical approach is available. We have the likelihood for all the parameters, including the tuning parameters τ ; the error term parameter σ 2 E , and the random function parameters β, θ, and σ 2 . Thus, all parameters can be estimated using the MLE method. We refer to this method as the full MLE. It can be done for each GP model defined in (6). The −2 times concentrated log likelihood function (except for constants) of all parameters for the combined data withβ B andσ 2 B plugged in is where Some other approaches based on the likelihood function, including the full MLE, were also proposed by Cox, Park, and Singer (2001). One of them is the "Separated MLE" (SMLE) method, which maximizes the conditional likelihood function of experimental data when the computer data is given. Here, the parameters β, θ, and σ 2 are estimated by maximizing the marginal likelihood for the computer data only. These are plugged into the conditional likelihood of the experimental data given the computer data. This is then maximized with respect to γ E and τ to obtain estimates of those parameters.
One advantage of the likelihood-based method is that it can simultaneously use both computer and experimental data to estimate τ , whereas ANLS method uses computer data only. These likelihood-based approaches enrich the tuning methods. Cox, Park, and Singer (2001), found the SMLE method to be better than the full MLE method. Thus, SMLE is compared with the proposed method in this study. The details of the SMLE are provided in the Supplemental Material.

Proposed method 4.1 Max-min algorithm
The following are the steps for the proposed tuning method. We call it a Max-min algorithm because it uses maximization and minimization iteratively.

Algorithm Max-min: iterate Step 3 and Step 4 until convergence is achieved
Step 1 (model building): build a surrogate (6) using the MLE for the given computer data only.
Step 3 (maximization): build a new surrogate (6), using the MLE for the combined data with the fixed τ obtained in the previous step.
Step 4 (minimization): set iteration i = i + 1, and find τ by minimizing RSS p (τ ) in (11) using the surrogate built in Step 3. If τ satisfies the stopping rule, then stop; otherwise, go to Step 3. Note that in each iteration of Steps 3 and 4, τ is updated; thus, the estimates of the parameters of θ, β, σ 2 E , and σ 2 are updated. We expect this to positively influence the finding of τ of Step 4. Steps 2 and 4 are the same in terms of minimizing RSS p (τ ), but Step 2 uses only computer data, while Step 4 uses either the combined or computer data. Steps 1 and 3 are the same in terms of obtaining the MLE by maximizing the likelihood function, but Step 1 uses only computer data, while Step 3 uses the combined data. The likelihood function in Step 3 is L(η B ; τ , y B , X B ), which is used for the full MLE method. But here, the tuning parameters are fixed as the τ that was obtained in the previous step. For the optimizations in Steps 2, 3, and 4, quasi-Newton numerical algorithms were used.
On using the combined data in Step 3, we assume that the functional relation between the inputs and output is the same for the computer code and physical process. The computer data and the experimental data are linked by the use of this common response function, with τ being the value for the tuning parameters in data from the physical process.
The Max-min algorithm stops when one of the following rules is satisfied: for i = 1, 2, · · · , 1 (maximum iteration): Number of iterations reaches the pre-assigned maximum number, where 'ftol' is a pre-assigned small value for tolerance. Here RSS p ( τ 1 ) is the minimum value of RSS p obtained in Step 2, and RSS p ( τ i ) is the minimum value obtained in Step 4 in the i-th iteration.
When the RSS p in Step 4 is greater than that of Step 2 or that of the last iteration, a small random fluctuation on τ is given before Step 3. Without this fluctuation setting, the algorithm stopped within four iterations. Based on our experience, the random fluctuation caused a reduction in RSS p . However, it did not make the algorithm execute more than 20 iterations.
When x 0 is a given site representing ( τ , x E i ) in the experimental data (X E ), the prediction formula needed for computing RSS p in Step 2 iŝ where f 0 is the known linear regression function vector; r 0C is the n C × 1 correlation vector between Y (x 0 ) and Y (X C ), and The prediction formula needed for computing RSS p in Step 4 iŝ where r 0B is the (n E + n C ) × 1 correlation vector between Y (x 0 ) and Y (X B ). Here,θ B , β B , andγ E are the MLEs from the combined data in Step 3. Note thatθ B andγ E are needed for constructing V B . Even though x 0 and X E are included in the construction of r 0B and V B , the predictionŶ (x 0 ) in (16) is not an exact interpolation because of the positiveγ E . Our approach was motivated by the iteratively re-weighted least squares method in regression analysis. The alternating estimation of parameters recursively is similar to the EM (expectationmaximization) algorithm. One can view this method as similar to a frequentist version of the Bayesian modularization approach (Liu, Bayarri, and Berger 2009). Some practical suggestions for the calibration of large-scale simulations and optimization are provided in Gramacy (2016, Section 4). However, our study was conducted independently from their study.
The advantages of the proposed Max-min algorithm are as follows. The uncertainty in the approximation of Y usingŶ in Step 4 is smaller than that in the ordinary ANLS, because the MLE method in Step 3 and the prediction model in Step 4 use combined data (a larger sample size), while the ordinary ANLS method uses only the computer data in Step 1 for MLE and in Step 2 for prediction. In addition, the Max-min algorithm accounts for the covariance between the computer data and experimental data in building the metamodel. Moreover, our test function experience revealed that the solutions from the Max-min algorithm are less influenced by the initial value of τ than by the ANLS method because the estimated τ in the Max-min is iteratively updated several times.

Approximated confidence region of estimates
Once τ has been estimated, some indications of the accuracy of the estimates are generally necessary. Thus, we rely on the asymptotic theory for the nonlinear least squares estimator that appears in the regression analysis (Draper and Smith 1981). The approximate 100(1 − α)% confidence region of τ is obtained by where F α (q, n E −q) is the upper (1−α) percentile of the F distribution with q and n E −q degrees of freedom (q is the number of parameters in τ ). Wong, Storlie and Lee (2017) considered a Bootstrap approach to compute the uncertainty of the estimates. 5 Toy-model study 5

.1 Exact computer models
In this section, we apply the methodology outlined in the previous sections to the seven test functions Y (τ , x) that are easy to compute. In the first five test functions, the experimental data with sample size n E are generated by where the random variable e follows a normal distribution with mean zero and variance σ 2 e . These five examples assume that the functional relation between the output and inputs is the same for both the computer code and the physical process. By contrast, in the last two test functions (6 and 7), the model bias term b(x) is added to the function Y (τ * , x) in (18). Here, τ * stands for the true tuning parameters and is specified for each toy model. No random error is given to the computer responses. Sample sizes for both data are set to 30 (n C = n E = 30) for the first five test functions. For generating the computer and experimental data, the same uniform distributions were used for the random numbers of x variables. The details on the eight test functions are described in the Supplemental Material.
To address the uncertainty resulting from the designs for the input variables, we repeated the τ estimation using 30 random Latin-hypercube designs. Thus, 30 sets of estimates of τ were obtained, and the averages and standard deviations are reported. Two different GP models of (6) were used with the correlation functions of (3) and (4), which are called "Model 1" and "Model 2," respectively. Figure 1 shows box plots of the distance to the true value in five test functions in which Model 2 is employed as a surrogate. The Max-min algorithm generally works better than the ANLS and SMLE methods. Tables S1 to S8 in the Supplemental Material present the results for each toy model, each method, and Model 1 and Model 2. Some columns show the average of Euclidian distances (Dist) between the true τ and estimates, with standard deviations in parentheses. The last column shows the mean squared error (MSE) of the estimates obtained using the following formula: where std(τ i ) is the standard deviation of each estimate computed from 30 repetitions. Figures S1 to S5 show the box plots for each test function and method. In terms of the average distance to the true value and MSE in Tables S1 to S5, the Max-min algorithm works better than ANLS for all toy models. Max-min also works better than SMLE for test functions 1, 2, and 4. Finally, SMLE works slightly better than ANLS, except in test function 2.    Figure S6 shows the typical convergence of the Max-min algorithm for 10 trials of test function 1. In most cases, the algorithm stops at the third or fourth iteration, except for a few cases wherein it stops at the sixth or seventh iteration. This means that the improvement of τ in the second iteration is significant, while it may not be significant after the second iteration. The small random fluctuation on τ , as mentioned in subsection 4.1, was not applied in this computation.

Another prediction
A potential drawback of the Max-min algorithm is that it partially fails to decouple between the estimation of τ and building of the metamodel because it uses the combined data in both Steps 3 and 4. To address this difficulty, instead of using equation (16) in Step 4, one can employ prediction (15) but withθ B and β B (or β C ): This leads to another version of the Max-min algorithm, which is applied to the next toy models. Note that this version does not useγ E . We denote this prediction asŶ C|B (x 0 ). We can employ β C instead of β B in (16), but we have not yet attempted to do so.

Inexact computer models
A computer model is often considered inexact. This means that the computer model does not perfectly match the real system even if some parameters included with it are optimal (Plumlee 2017). This is the computer model bias; it is the difference between the model and the truth. The model-bias term b(x) is included in the following two test functions. These are employed to check the performance of the proposed method for cases where the functions generating the output for the computer code and for the physical experiment differ, but we incorrectly assume that (7) and (8) hold.
Test function 6: n C = n E = 20 Test function 7: n C = n E = 20 Computer data : Test function 6 is modified from Plumlee (2017). Test function 7 was used in Bastos and O'Hagan (2009), in Goh et al.(2013), and in Gramacy (2016). Twenty random Latin hypercube designs were used repeatedly to address the uncertainty.
In the previous exact toy-model study, we used the difference between τ and true τ for the physical process as a measure of performance. This requires that there be a true τ , but there would be no such true τ in an inexact computer model. In this situation, the minimum value of RSS p would be a more appropriate measure.
In calculating RSS p in the Max-min algorithm for the above two test funtions, we set ftol = 1.e-4, and maxagain = 7. For the small fluctuation of τ , random numbers from N (0, σ 2 τ ) were used where σ τ = max( τ × 0.1, 0.3). Figure 2 shows parallel coordinated box plots of RSS p values for test function 6 computed from 20 Latin hypercube designs. RSS p values are calculated using the ANLS method and the Max-min algorithm with GP Model 1 and Model 2 via two different predictions. The predictions byŶ B in (16) andŶ C|B in (20) in calculating RSS p in the Max-min algorithm were attempted for the sake of comparison. The Max-min algorithm worked better than ANLS in most cases. The left panel of Figure 3 is an ANOVA-type main-effect plot based on mean values, which shows improvements from ANLS to Max-min, from Model 1 to Model 2, and fromŶ C|B toŶ B . The mean values from the last two predictions (IB0 and IB1) are calculated within the Max-min algorithm.
Relative improvement (RI) from ANLS to Max-min is provided in Table 1. RI is defined by     For tuning test function 7, we considered a simple bias correction (BC) technique. Two constants for additive and multiplicative corrections were set for the predictor (Fernandez-Godino et al. 2016) asŶ This naive BC method may not be better than sophisticated methods such as that using ρ(x) and δ(x), or the approach by Kennedy and O'Hagan (2001), but we hope for this to work better   Figure 4: This is the same as Figure 2, but for test function 7, calculated without (non-BC) and with bias correction (BC).
than the non-BC case. The objective function RSS p in (11) is now changed to Then, τ , ρ, and δ are estimated simultaneously by minimizing RSS bc p in both the ANLS and the Max-min methods. Figure 4 shows parallel coordinated box plots of the RSS p values for test function 7, computed from 20 Latin hypercube designs. The left panel is for non-BC, while the right one is for BC. The acronyms in this figure are the same as in Figure 2. This figure shows that the Max-min algorithm worked better than ANLS in most cases. The right panel of Figure 3 is a main-effect plot for test function 7. It shows improvements from ANLS to Max-min, Model 1 to Model 2, non-BC to BC, andŶ C|B toŶ B . Note that the patterns of the main effects for test functions 6 and 7 are very similar. The RI from ANLS to Max-min is provided in Table 1. The bias correction was most effective when used with Model 2 and with the predictionŶ B .
6 Application to nuclear fusion model 6

.1 Nuclear fusion data
A simple measure of energy efficiency in a nuclear fusion device (called a tokamak, from the Russian language) is the global energy confinement time τ E . The theoretically based confinement model can be written as follows (Kay and Goldston 1985): P, I, N, B) , (24) where f is a known function calculated using a complex simulation code called Baldur, P is the total input power, I is the plasma current, N is the electron density, B is the magnetic field, and τ = ( τ 1 , τ 2 , τ 3 , τ 4 ) are the following adjustable parameters that determine energy transfer, that is, drift waves, rippling, resistive ballooning, and the critical value of η i (which provokes increased ion energy losses for the drift waves), respectively. The experimental data comprises only P, I, N, B, and the real observation y E , whereas the computer data comprises eight independent variables (T , P, I, N, B) and computer response y C obtained using Baldur. The experimental data were drawn from the database of S. Kaye: 42 observations from the PDX (Poloidal Divertor Experiment) tokamak in Princeton and 64 from the Baldur simulator (Singer et al. 1988). Because the Baldur simulator requires five CPU minutes on a Cray supercomputer for one execution, a careful selection of input points is required, which is a statistical design problem for a complex simulation code. For this purpose, we used a data-adaptive sequential optimal design strategy, which is described briefly in the Supplemental Material.  Table S9 in the Supplemental Material presents the MLE of the parameters of GP Model 1 obtained from the combined data. Table 2 provides the results of τ estimation obtained using ANLS, SMLE, and the Max-min methods on the basis of Model 1. A quasi-Newton optimization routine was employed in searching for τ . Several starting values were tried to avoid the local minima. The last column in Table 2 shows the value of RSS p at the convergence of the algorithms. The RSS p for SMLE was obtained by calculating RSS p for the estimated τ . The RSS p for the Max-min with Model 1 is the smallest among the three methods. Figures 5 and 6 show the residual plot (residual vs. predicted values) and confidence regions of the tuning parameters (τ 1 and τ 4 ) for the tokamak data, which were obtained using the Maxmin algorithm with Model 1. In Figure 5, the predicted values for the computer data (circles) had a wider horizontal range than those for the experimental predicted values (crosses). In addition, the residuals for the computer data had a narrower vertical range than those for the experimental data. These results indicate that the computer data are fitted relatively well, and that good coverage of the range of experimental observations is ensured.

Summary and Discussion
Using a GP model, we considered an iteratively re-estimated ANLS method to tune complex computer code to data, namely a Max-min algorithm. This method is an extension of the ANLS method. A simulation study using toy functions suggests that the proposed Max-min algorithm works better than the ANLS and SMLE methods. We applied this technique to a computational nuclear fusion model. The proposed method can be useful for other applications in disciplines in which unknown theoretical parameters must be estimated using complex computer codes and real experimental data. Moreover, the Max-min algorithm is also applicable to other metamodels, such as spline (Wong, Storlie and Lee 2017), support vector machines, and neural networks. The Max-min algorithm requires more computing time than the ANLS method. This is because Max-min uses the combined data, resulting in a larger correlation matrix that must be inverted. It also requires more iterative minimizations of RSS p compared to the ANLS method, which requires just one iteration. Two versions of predictions were tried in the last two toy models wherein bias functions were added. The prediction based on both data sources (Ŷ B ) reduced RSS p (residual sum of squares of prediction) more than that based on computer data given the parameter estimates obtained from both data sources (Ŷ C|B ). Even thoughŶ B reduced RSS p more thanŶ C|B , the improvement of the estimation of true τ compared toŶ C|B is not known. A simple bias correction was considered in test function 7. This bias correction was the most effective when used with GP Model 2 and with the prediction based on both data sources.
We tried two different GP models in the toy-model study, and found that the estimates of tuning parameters were changed according to the selected GP model. Model 1 in (6) is a superimposition of a GP with correlation function (3) on the first-order regression model. We believe that the first-order regression part sometimes works well to fit the dominant relationship, even though the isotropic correlation (3) may be unrealistic for capturing the variability of responses from the dominant relationship. It is computationally easier to obtain the MLE when the number of input variables is large. Model 2, with (4) in (6) assumes a separable Gaussian correlation, which is a popular family of correlation models as per the literature (Santner, Williams and Notz 2018), but it is computationally expensive to obtain the MLE if many inputs exist.
A model selection procedure among many GP models with various combinations of nonzero β j s and θ j s in (6), as in Marrel et al. (2008), Chen et al. (2016), and Lee and Park (2017), may lead to realization of a better surrogate model. In this study, we did not consider the identifiability of the tuning parameter, which is an important focus of recent developments (e.g., Plumlee 2017; Tuo and Wu 2018).
There are, however, basic limitations with tuning computer code to real-world data regarding the experimental design. We found that the performance of tuning methods is significantly dependent on the designs for both the computer experiments and the physical experiments. Some authors, including Cailliez, Bourasseau, and Pernot (2014), and Beck and Guillas (2016), explored this topic. We agree that a sequential tuning approach is practically useful, as in Pratola et al. (2013), Kumar (2015), and Damblin et al. (2018). Further research on relevant designs under the sequential tuning approach will be helpful.

A.1 Maximum likelihood estimation in Gaussian process model
Once the data have been collecte at the observation sites {x 1 , ..., x n }, we use the maximum likelihood estimation (MLE) method to estimate the parameters in linear model part and covariance function. Since we assume that y(x) is a Gaussian process with mean F β and covariance matrix σ 2 R, the likelihood function of y is where F is a so-called design matrix. When the covariance parameters θ and γ C are specified, the MLEs of σ 2 and β are denoted by Here,β is the generalized least squares estimator of β. Since the likelihood equations do not lead to a closed-form solution, a numerical optimization procedure is required. The Cholesky decomposition V = U t U is used as a major computation in calculating the likelihood function, where U is an upper triangular Cholesky factor. The computational details of calculating and minimizing negative log-likelihood function are provided in Park and Baek (2001). One can use a R program DiceKrig (Roustant 2012).

A.2.1 Computer and experimental data
For notational convenience, experimental data is denoted by the subscript "E" and computer simulation data by subscript "C." Let τ be an adjustable parameter vector to be estimated. Let T be the input variables of the computer code corresponding to τ . Here, τ is a vector of the deterministic tuning parameters and T is a vector of the random variables. The original experimental input variables is denoted by X. Let q and p be the dimensions of τ and X. Further, let n C , n E be the number of observations. Then, we have the data matrix of the independent variables: X C and X E for computer and experimental data; Here, t ij in X C represents the j-th value of the i-th T variable (T i ) and x Eij and x Cij denote the j-th value of the i-th X variable of experimental (X Ei ) and computer (X Ci ) input. X E is a n E × (q + p) matrix and X C is a n C × (q + p) matrix. Note that the first part of X E is composed of the unknown parameters τ 1 , · · · , τ q , while the corresponding part of X C comprises input values (t ij ).

A.2.2 Combined data
The following notations for combined computer and experimental data are introduced: for the data matrix of the independent variables; the so-called "design matrix, defined as the functions of the values of input variables; and the computer responses and real observations, respectively. Here, the subscript "B" indicates the combined "both" computer and experimental data. Note that variables T are incorporated in the simulation code as design sites. X C and F C contain T , while X E and F E are the functions of the unknown parameters τ . The Gaussian process model is now simultaneously applied to computer and experimental data. Let η = (τ , θ, γ C , γ E , σ 2 , β), where γ C = σ 2 C /σ 2 and γ E = σ 2 E /σ 2 , which are the variance ratios for the computer and experimental data. Here, σ 2 C and σ 2 E are the variances of error term ( ) in the Gaussian process model for the computer and experimental data, respectively. When necessary, β C and β E are used to denote the regression coefficients for the computer and real experimental data. Then, given the independence and normality assumptions, we have where where R(X C , X E ) represents a n C × n E matrix composed of the correlations computed between X C and X E . Note that V B is a n B × n B positive definite covariance matrix for the combined data, where n B = n C + n E . We set γ C = 0 because only a deterministic computer model is considered in this study.

A.3 Separated MLE for code tuning
For the details of SMLE, we make use of the conditional distribution of the experimental data given the computer data, which is normally distributed with mean and covariance where covariance matrices V 's are given as in (32). In these formulae, we suppressed the parameter dependencies in µ E = F E β E , µ C = F C β C , V CE , V CC , and V EE . Now the −2 times concentrated log conditional likelihood function (except for constants) withβ andσ 2 E|C plugged in is −2 log L(τ , γ E ; X E , y E | y C ,γ C , β C ,θ, X C ) = n E logσ 2 E|C + log |V E|C |, A. 4 Test functions for toy model study Test function 1: nC = nE = 30 for all five test functions.
In each test function, n C values of (T , x) were selected as inputs for the "computer code", that is, the function Y (T , x) is evaluated at n C values. The inputs were chosen to be well spread around a reasonable space known to potentially contain the true parameter value. n C computer data points and n E experimental data points were generated by using random Latin-hypercube designs, except the τ values were used instead of T . In real application situation, one would consider more sophisticated designs such as data-adaptive sequential optimal experiments as reported in the next section.

A.5 Sequential designs in tuning a nuclear fusion simulator
For given a Gaussian process model, the A-optimal design is obtained by minimizing the integrated mean squared error of prediction (MSEP) with respect to a design S , where Q is the design region, and µ is a "weight function" which may be the empirical measure of uniformly distributed random points. Note that neither M SEP nor IM SE depend on the unknown parameters β and σ 2 , but depend on θ, γ C and design S. This makes it possible to design an experiment (for specified values of θ and γ C ) before taking the data. Because θ is generally not available for the initial design stage, in our example, we used an rough estimate of θ based on a previous similar work given by a Baldur specialist. Initially we found 10 optimal design points for eight variables (T , P, I, N, B) which minimize IM SE over the design region Q, with θ = .5 and γ C = .001.
The following is the data-adaptive sequential optimal design procedure that we have used in this study.
Step 1. Collect computer observations based on the given optimal design (n 1 points, say).
Step 2. Find an appropriate model and estimates of parameters (θ, β, σ 2 , γ C ) from the computer data.
Step 3. Check the MMSE, and stop constructing the next stage design if MMSE is smaller than a preassigned target value. Otherwise, go to the next step.
Step 4. Use the estimates and model found in Step 2 to choose the next stage optimal design (n 2 points, say) under the condition that the previous design is given (i.e, update n 2 more points to the previous design to make n 1 + n 2 points).
Step 5. Collect (n 2 ) more observations, and go to Step 2.
The MMSE (maximum mean squared error of prediction) is defined as where x i , i = 1, 2, · · · , K, are d-dimensional random vectors. Note that the MMSE is used here as a measure of accuracy of a given prediction model.        Regression coefficient for T 1 -0.027 β 2 Regression coefficient for T 2 -0.010 β 3 Regression coefficient for T 3 0.001 β 4 Regression coefficient for T 4 -0.015 β 5 Regression coefficient for P -0.031 β 6 Regression coefficient for I 0.009 β 7 Regression coefficient for N -0.012 β 8 Regression coefficient for B 0.017 σ 2 Variance of Y 1.506E-04 γ E Variance of σ 2 E /σ 2 0.454  Figure S1 but the test function 2.  Figure 9: Same as Figure S1 but the test function 3.  Figure 11: Same as Figure S1 but the test function 5.