{'fraction_non_alphanumeric': 0.03977314040572966, 'fraction_numerical': 0.010906711262997163, 'mean_word_length': 4.458821194681485, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Invariant Gaussian Process Latent Variable Models and Application in Causal Discovery

In nonlinear latent variable models or dynamic models, if we consider the latent variables as confounders (common causes), the noise dependencies imply further relations between the observed variables. Such models are then closely related to causal discovery in the presence of nonlinear confounders, which is a challenging problem. However, generally in such models the observation noise is assumed to be independent across data dimensions, and consequently the noise dependencies are ignored. In this paper we focus on the Gaussian process latent variable model (GPLVM), from which we develop an extended model called invariant GPLVM (IGPLVM), which can adapt to arbitrary noise covariances. With the Gaussian process prior put on a particular transformation of the latent nonlinear functions, instead of the original ones, the algorithm for IGPLVM involves almost the same computational loads as that for the original GPLVM. Besides its potential application in causal discovery, IGPLVM has the advantage that its estimated latent nonlinear manifold is invariant to any nonsingular linear transformation of the data. Experimental results on both synthetic and realworld data show its encouraging performance in nonlinear manifold learning and causal discovery.


INTRODUCTION
Latent variable models, including dynamic models, aim to find a low-dimensional underlying manifold that helps understand the structure of the data. For instance, factor analysis and probabilistic principal component analysis (PPCA, Tipping & Bishop, 1999), as linear models, find the subspace of the common factors such that the error terms are uncorrelated from each other. Nonlinear latent variable models, such as the Gaussian process latent variable model (GPLVM, Lawrence, 2005), look for an intrinsically simple nonlinear manifold. In dynamic models, such as the Kalman filtering (Kalman, 1960) and Gaussian process dynamic model (GPDM, Wang et al., 2008), the latent processes are expected to be modelled well by a given temporal model.
In nonlinear latent variable models or dynamic models, usually the observation noise is assumed to be independent across data dimensions. This implicitly assumes that the interesting manifold coincides with the common subspace of the data, which is not necessarily true in practice. Consequently, the modelling performance and recovered latent variables are sensitive to the dependencies of the noise terms; in particular, they may be quite different for the data after linear transformations. This can be avoided by a nonlinear latent variable model that is able to adapt to arbitrary noise correlations, which will be developed in this paper.
On the other hand, the dependence between the noise terms is useful for data modelling, prediction, and causal discovery in the presence of latent common causes (confounders), 1 which is of particular interest to us. For example, in human motion tracking, different parts of humans may be closely related, or even causally related, almost instantaneously. By exploiting these relations, one could model or predict the joint distribution of the human parts better. From the causal discovery point of view (Spirtes et al., 2001;Pearl, 2000), making use of the error dependencies enables us to identify the linear causal relations among the observed variables, with the latent variables considered as confounders.
In this paper the latent variable model of interest is GPLVM, a fully probabilistic one. Due to the Gaussian process (GP) prior (Rasmussen & Williams, 2006) on the nonlinear mapping, it is able to control the complexity of this mapping automatically, and can easily give probabilistic reconstructions of the data. It has been reported to work well with small data sets (Lawrence, 2005). However, to avoid the computational difficulties associated with a large kernel matrix which is caused by a dependent noise model, GPLVM assumes that the noise is independent across the dimensions. We circumvent these difficulties by putting the GP prior on a certain linear transformation of the latent functions. Consequently, we extend GPLVM (as well as GPDM, which is a combination of GPLVM with the GP prior on the temporal structure of the latent variables) to allow arbitrary noise covariances, while the algorithm still has similar computational loads as the original GPLVM.
The extension of GPLVM proposed by us has some appealing properties. First, its modelling performance and estimated latent variables are invariant to any nonsingular square linear transformation of the original data. Second, it admits a hierarchical way to explain the data generating process: the latent variables are a common cause to all observed variables, and the estimated noise dependencies imply further linear relations between the observed variables. If one believes such relations to be causal, the proposed method provides a suitable tool for causal discovery in the presence of nonlinear confounders.

NONLINEAR LATENT VARIABLE MODELS AND CAUSAL DISCOVERY
Let y t = (y 1t , y 2t , ..., y Dt ) T , t = 1, 2, ..., N , be a sequence of vector-valued observed variables. The nonlinear latent variable model considered here is where x t = (x 1t , x 2t , ..., x dt ) T are latent/hidden variables, and the observation noise e t is independent from the latent variables. In dynamic models, there is another model to specify the dynamic structure of the hidden variables. The kth order dynamic can be described as x t = f (x t−1 , ..., x t−k ; η) + ε t , where ε t denotes the process noise.

RELATION TO CAUSAL DISCOVERY WITH NONLINEAR CONFOUNDERS
In some scenarios, we are more interested in the linear (causal) relations between the observed variables y it than the common effect of the latent variables. To this end, we describe the generating process of y t as where B contains the coefficients of the linear instantaneous influences, and the entries on its diagonal are constrained to be zero. If B can be transformed to strict low-triangularity by row and column permutations, its implied linear causal relations are acyclic. The model above is an extension of the linear causal structural equation model (Pearl, 2000) in that it has incorporated the effect of latent variables. Here it is assumed that the disturbancesẽ it are independent from the latent variables (or confounders), as well as from each other.
Equation (2) can be written as where g (I − B) −1 ·g and e t (I − B) −1 ·ẽ t . This is exactly the latent variable model given by (1), but the noise terms e it are generally not independent.
One can then see that the linear causal relations between the observed variables are implied in the structure of the noise. Suppose that the noise in the latent variable model (1) has been estimated. If the disturbancesẽ it are non-Gaussian, according to the independent component analysis (ICA, Hyvärinen et al., 2001) theory, (I−B) can be found up to certain permutation and scaling indeterminacies. Moreover, if the causal relations are acyclic, B is fully identifiable (Shimizu et al., 2006). In fact, this condition can be relaxed; under some assumptions, B may be identifiable even if there exist feedbacks between y it (Lacerda et al., 2008).
Ifẽ it are Gaussian, the observation noise e t is Gaussian with the covariance matrix cov(e t ) = (I − B) −1 Σẽ(I − B) −T (where Σẽ is the covariance matrix ofẽ t , which is diagonal). In this case, (I − B), as well as B, is not identifiable, i.e., the whole causal model can not be fully determined. One can then use conditional independence-based approaches to find the Markov equivalence class of the causal model (Spirtes et al., 2001;Pearl, 2000), or the linear-Gaussian Markov network (Lauritzen, 1996) of the observed variables, with the effect of the latent variables excluded.

A GENERAL PROCEDURE FOR LINEAR CAUSAL DISCOVERY WITH NONLINEAR CONFOUNDERS
To estimate the linear (causal) relations between observed variables, with the latent variables considered as confounders, we first need to estimate the latent variable model (1) which is able to adapt to arbitrary noise dependencies. The estimated noise termsê it are then obtained.
Next, depending on the problem and target, one can find the graphical model or causal model implied by the errorsê it . One can always obtain the Gaussian Markov graph of y it , which is implied by the precision matrix (the inverse covariance matrix) ofê it . If we would like to find an acyclic causal model, we need to test ifê it are Gaussian. If they are, traditional conditional independence-based approaches can be used to find the Markov equivalence class of the causal model. If they are not, one can apply the linear, non-Gaussian, and acyclic model (LiNGAM, Shimizu et al., 2006) based analysis onê t to find the causal influence matrix B. The case where the influences are not necessarily acyclic is more difficult to handle, although there exist certain methods, such as that proposed in Lacerda et al. (2008), to identify B.

RELATION TO PREVIOUS WORK
Causal discovery in the presence of confounders is a very challenging problem, since it is difficult to find the confounders. In the linear confounder case, Hoyer et al. (2008) exploited overcomplete ICA to solve this problem. In the nonlinear case, Janzing et al. (2009) proposed a heuristic but useful method to examine if there exists a nonlinear confounder for two given observed variables. The model they exploited is essentially a nonlinear latent variable model with independent noise terms, and it inspired our work.

INVARIANT GAUSSIAN PROCESS LATENT VARIABLE MODELS
We aim to develop a reliable nonlinear latent variable model which can adapt to arbitrary noise correlations. Now the problem is to find a suitable form for the nonlinear mapping g in (1). Recently, the GP prior has been adopted to model it, resulting in GPLVM (Lawrence, 2004;Lawrence, 2005). This model is typically different from traditional nonlinear dimension reduction models in that it marginalizes the mapping instead of the latent variable. The latent variables are estimated together with the parameters appearing in the GP prior by maximizing the marginal likelihood.
In the original GPLVM (Lawrence, 2005), the GP prior is directly put on the nonlinear function g in (1). If the noise terms e it are dependent, the involved kernel matrix is of the dimensionality DN × DN , which is very large and causes computational difficulties. Hence, the noise e t in GPLVM is assumed to have uncorrelated Gassian components (or more frequently, to have isotropic Gaussian components). As explained in Section 2, this prevents GPLVM from being a suitable tool for causal discovery with nonlinear confounders.
To tackle this problem, we extend GPLVM to allow arbitrary noise correlations. From the nonlinear dimension reduction point of view, the latent nonlinear manifold estimated by the extended model is invariant to any nonsingular square linear transformation of the data. Hence we name the extended model invariant GPLVM (IGPLVM). Two approaches for IGPLVM are given below.
3.1 Approach I: AN EFFICIENT SCHEME

Model Formulation
Here we still assume that the noise e t is Gaussian; treatment of the non-Gaussian noise will be discussed later. For a parsimonious parametrization, we use the Cholesky factorization to decompose the noise covariance matrix as cov(e t ) = LL T , where L is a lowertriangular matrix with strictly positive diagonal entries. That is, the noise is written as where e * t has independent components with an identity covariance matrix.
Combining (1) and (3) gives Let Y i be the column vector of the values of y it , i.e., If one assumes that the components of L −1 g(x t ; θ), as nonlinear functions of x t , can always be modelled with independent GP priors, i.e., we have p( , y I t can then be modelled by the original GPLVM (Lawrence, 2005), with independent noise terms and the noise variance fixed to one. Since y t = Ly I t , we have Y = (L ⊗ I N )Y I , where ⊗ denotes the Kronecker product and I N denotes the N -dimensional identity matrix. Consequently, log p(Y|X) = log p(Y I |X) − N log |L|. Finally, the data log likelihood of (4) is where r and γ are parameters involved in the kernel matrix K y I , which is given below. Bear in mind that the noise in y I it has a unit variance; we can then use a "RBF kernel" of the following form for Y I i : 2 is the (t, t ′ )th entry of K y I . Parameters X, L, r, and γ can then be estimated by maximizing the likelihood.
Thanks to the linear transformation L in (4), one can treat the components of y I t independently, and consequently the likelihood (5) involves the N × N kernel matrix K y I . If one directly applies independent GP priors on the components of g(x t ; θ) (as in the original GPLVM), the kernel matrix appearing in the likelihood would be of the dimensionality DN × DN .

Learning and Discussions
Define L L −1 . It is also lower-triangular with positive diagonal entries. We alternate between maximizing the likelihood (5) w.r.t. r, γ, and X using the scaled conjugate gradient (SCG), which is similar to the parameter estimation procedure in the original GPLVM, and that w.r.t. L, which is obtained in closed-form in each iteration. Let y I (Y I 1 , Y I 2 , ..., Y I D ) T , which is of the size D × N , and similarly for y. The second term of l 1 given in (5) can then be written as − 1 2 Tr( LyK −1 y I y T L T ). In each update iteration, at the stationary point w.r.t. L, we have Multiplying the above equation by L T from the right gives LyK −1 y I y T L T = N I N . That is, given the estimate of the remaining parameters, the optimal value of L is the inverse of the (lower-triangular) Cholesky factorization matrix of 1 N yK −1 y I y T . In our implementation, PCA is used for initialing the latent variables x t , and both r and γ are initialized to one.
2 Strictly speaking, y I it , i = 1, 2, ..., D, may have different signal-to-noise ratios (SNR's), and consequently the parameter r in the kernel may vary for different dimensions of y I it . Here for simplicity we assume that all y I it have the same SNR, like in Lawrence (2005); Wang et al. (2008).
Compared to the scaled GPLVM (Grochow et al., 2004), Approach I for IGPLVM involves D(D + 1) − D = D(D−1)/2 more parameters, due to L. But since these parameters are obtained in closed-form in each iteration, the computational load remains almost the same.
3.2 APPROACH II: AN EXACT BUT LESS EFFICIENT TREATMENT

Model Formulation
In the original GPLVM (Lawrence, 2005), as well as GPDM (Wang et al., 2008), the components of g are treated independently; in Approach I for IGPLVM (Subsection 3.1), we put independent GP priors on the components of L −1 g(x t ; θ) in (4). In practice these nonlinear functions are not necessarily independent.
To make the modelling performance of g completely insensitive to linear nonsingular transformations of the data, we further simply model g = (g 1 , ..., g D ) T as a linear transformation R of some functions with independent GP priors, denoted by g * = (g * 1 , ..., g * D ), 3 i.e., This idea has been used in the semiparametric latent factor models for GP regression with multiple response variables (Teh & Seeger, 2005). However, the linear transformation R adopted here is different from Teh and Seeger (2005) in two aspects. First, for computational reasons, we constrain R to be square. Second, we give R a parsimonious parametrization. Without loss of generality, we parameterize R as a lower-triangular matrix with positive diagonal entries. This scheme has been frequently used to parameterize the loading matrix in factor analysis; for instance, see Geweke and Zhou (1996).
Combining (4) and (7) gives the observation equation: We construct G * i and G * from g * it in the way that Y and X were constructed in Subsection 3.1.1. Let K g * be the kernel matrix shared by all g * i , i = 1, ..., D. We have p(G * i ) ∼ N (0, K g * ). As g * i are assumed to have independent priors, one can see p(G * ) ∼ N (0, K G * ), where K G * = I D ⊗ K g * = diag(K g * , ..., K g * ). We use a "RBF kernel" for the kernel shared by all g * i : As in Lawrence (2005), there is a redundancy between γ and the overall scale of x t ; however, we still keep γ, because it is expected to compensate the overall scale of x t quickly.
Let y * t g * (x t , θ) + R −1 Le * t . Equation (8) can then be written as y t = Ry * t . For simplicity of the notation, we define L R R −1 L and R R −1 . Let Σ = L R L T R . One can see 4 Furthermore, since y t = Ry * t , we have Y = (R ⊗ I N )Y * , such that log p(Y|X) = log p(Y * |X) − N log |R|. Therefore, the data log likelihood is

Learning and discussions
All involved parameters, including X, L R , R, and γ, are learned by maximizing the likelihood (10). We divide all parameters into two subsets. One contains free parameters in R and L R (both of which are constrained to be lower-triangular), and the other contains X and γ, which is involved in the kernel defined in (9). We alternatively update them by maximizing (10) using SCG for 10 iterations.
This approach is not suitable for relatively large datasets due to the high computational load in calculating the inverse and determinant of K Y * , whose dimensionality is DN × DN . Generally speaking, the complexity of finding K −1 Y * or det(K Y * ) would be O(D 3 N 3 ). One can exploit the informative vector machine (IVM, Lawrence et al., 2003) framework or other more recent techniques to select an active set of the training data and reduce the computational load, by extending the methods given by Lawrence (2005); Teh and Seeger (2005); Lawrence (2007). Furthermore, Approach II has D(D +1)/2−1 more parameters than Approach I given in Subsection 3.1. Consequently, it might be more prone to local optima and overfitting. However, this approach is still presented here, as an exact formulation for IGPLVM. In practice we found that Approach I plays a nice trade-off of the model capacity and simplicity.

NON-GAUSSIAN CASE
As claimed in Section 2, in the case that the observation noise is non-Gaussian, it might be possible to fully identify the linear causal relations among the observed variables, with the latent variables treated as confounders. Here we give the basic idea of how to estimate IGPLVM with non-Gaussian noise.
When the noise is Gaussian, the marginal distribution p(Y|X, θ) is in closed form; see (5) or (10). However, generally speaking, with a non-Gaussian distribution for e it , this marginal distribution could not be expressed in closed-form any more, causing difficulties in parameter estimation and inference. A straightforward way is to resort to the Laplace approximation, which approximates the posterior of the latent functions with a Gaussian density. Although theoretically rather simple, this approach would lead to complicated algorithms for IGPLVM.

Two-Step Method
For practical reasons we propose a simple heuristic two-step method to solve this problem. In the first step, we decompose the data into the nonlinear manifold and the noise effect, i.e., estimate the latent processes x t and the noise e t , with the Gaussianity-based IGPLVM. In the second step we consider the non-Gaussianity of the noise. We test ifê it are Gaussian. If they are not, apply ICA (Hyvärinen et al., 2001), which makes use of the non-Gaussianity, to further find their linear structure, i.e., to estimate the matrix L in (3). 5 (In this case, L is a full matrix, instead of a lower-triangular one.)

SIMULATIONS
Invariant GPDM (IGPDM) can be obtained straightforwardly by combining the proposed IGPLVM with the autoregressive GP regression model for the latent processes. The temporal structure of those processes helps estimate them more reliably. We use simulations to compare the performance of IGPDM with the original GPDM (Wang et al., 2008) and demonstrate the behavior of the two-step method for IGPDM with non-Gaussian noise proposed in Subsection 3.3.
The simulated data were constructed in the following way. The two-dimensional underlying latent pro-cess x t smoothly moved on a circle; see Figure 1(a) for its trajectory. The sample size was N = 400. The eight-dimensional observed process y t was generated according to the nonlinear latent variable model (1); components of g were random mixtures of various smooth functions, including linear, quadratic, thirdorder polynomial, and tanh functions. The error e t has a special structure; it was constructed as a linear transformation of some independent variables, i.e., e t = As t . s it were obtained by passing independent and i.i.d. Gaussian samples through power nonlinearities with the exponent between 1.5 and 2. Thus the noise is non-Gaussian (supper-Gaussian, in fact). The variances of s it were randomly chosen between 0.2 and 1. The mixing matrix A was lower-triangular with random entries (between -1 and 1) in the strictly lower-triangular part, and the entries on its diagonal were set to one. Consequently, if x t is considered as a confounder for all y it , the causal relations between y it are linear, acyclic, and fully connected. Furthermore, to make some components of e t significantly dependent on each other, we made three rows of A have similar but different entries.
Our goals were 1) to see if IGPDM could successfully estimate the nonlinear manifold and noise, 2) to examine if the two-step method for IGPDM with non-Gaussian noise could recover the noise structure, and 3) furthermore, to see if the causal relations between y it introduced by the noise could be identified by the causal discovery procedure given in Subsection 2.2. We used two measures for performance evaluation. One is the mean square error (MSE) between the true nonlinear functions g i and the estimated ones; it measures how accurately the nonlinear manifold is recovered. The other is the Amari performance index (Hyvärinen et al., 2001), denoted by P err , which is frequently used for source separation performance evaluation. It is used here to measure how accurately the matrix A, which generated the noise, can be estimated by applying ICA onê it (which is step 2 of the two-step method given in Subsection 3.3). The smaller this index, the better the performance.
We repeated the simulations for 5 random replications (we did not run more replications because Approach II for IGPDM given in Subsection 3.2 is timeconsuming). All methods were initialized by PCA. Table 1 reports the performance. One can see that due to the impact of the noise dependencies, it is not easy for the original GPDM to recover the smooth latent manifold. In contrast, incorporation of the noise structure enables the proposed methods, especially Approach I (Subsection 3.1), to estimate the latent variables as well as the noise very well. Approach II seems not as stable as Approach I, meaning that it may be prone to local optima or over-fitting. Approach I plays a nice trade-off of the simplicity and capacity of the model. The LiNGAM (Shimizu et al., 2006) software 6 was then used to analyze the estimated noise to find the causal relations between y it . For the noise estimated by Approach I for IGPDM, in all five replications, the software always reported thatê t nicely follows the LiNGAM model. In average 71% of the causal relations were recovered, and the chance of producing spurious relations was 7%. This is consistent with the good performance in the recovery of A, as shown by the small values of P err . For the noise estimated by the original GPDM, in all replications, LiNGAM reported thatê t does not follow LiNGAM or that the causal relations are only somewhat acyclic. Figure 1 (b,c,d) show the the learned latent processes in the first replication. Note that the comparison is fair: different methods employed the same prior on the latent process. Both methods for IGPDM successfully recovered the subspace of the latent process. Approach I was the best. We then applied a random rotation to the observed processes y t and learned the models for another time. The leaned latent processes are shown by crosses. Clearly the results of IGPDM remained the same (note that for Approach II, the results slightly changed, which was caused by different parameter initializations.) However, the results learned by the original GPDM changed a lot. This verifies that the performance of the proposed IGPDM (as a special case of IGPLVM) is rotation-invariant.

ON MOTION TRACKING
We then applied IGPLVM on a human motion capture data set extracted from CMU mocap database. We used the data of two walk cycles (subject 7, trial 1) for analysis. The same data have been used by Wang et al. (2008). Our goal was to see if the proposed IGPDM could find appealing relations between different joints. To make such relations more obvious, we further simplified the skeleton used by Wang et al. (2008), only keeping the limbs, as shown in Figure 3. Finally 22 variables were used.
To save space, the learned latent process is not given. As it is difficult to believe that the joints have acyclic causal relations, we are interested in the learned noise covariances and the implied Gaussian Markov network. Figure 2 shows the precision matrices (inverse covariance matrices) of the observed variables, of the error estimated by Approach I for IGPDM, and of the error estimated by the original GPDM. Zero entries in the precision matrix indicate missing edges in the Markov network. One can see that the network implied by

DISCUSSIONS
In this paper we proposed an extension of the Gaussian process latent variable model (GPLVM), namely, invariant GPLVM (IGPLVM). It has the appealing property that its estimated nonlinear manifold is invariant to any nonsingular linear transformation of the data. Since this model allows arbitrary noise correlations, we show that it also provides a suitable tool to solve an inherently difficult causal discovery problem -causal discovery in the presence of nonlinear confounders. The proposed Approach I for IGPLVM is recommended; it involves almost the same complexity as the original GPLVM, and we found that in practice this approach, in combination with suitable priors on the temporal structure of the latent variables, is able to estimate the latent functions and noise structure (including the implied causal relations).
In the original GPLVM as well as the proposed extensions, the model parameters in the nonlinear mapping and the latent variables are estimated simultaneously. As the latent variables consist of many values, the learning process may be prone to multiple local optima and overfitting, if ones does not have suitable temporal constraints on the latent variables. As reported by Lawrence (2005), the result of GPLVM may be sensitive to initializations. This is inherently related to the difficulties in non-Gaussianity-based causal discovery with linear confounders . However, suitable priors on the temporal structure of the latent processes help alleviate this problem significantly, as shown in our experiments, in which au- toregressive Gaussian process models were used for the latent processes. To improve the applicability of IG-PLVM for causal discovery in various practical problems, one line of our future research is to investigate the possibility of developing more reliable algorithms for IGPLVM without temporal constraints on the latent variables.