{'fraction_non_alphanumeric': 0.022755717373990215, 'fraction_numerical': 0.0036978040732734098, 'mean_word_length': 4.946887686062246, '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': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Learning non-linear structures with Gaussian Markov random fields

Nowadays, one of the most changeling points in statistics is the analysis of high dimensional data. In such cases, it is commonly assumed that the dimensionality of the data is only artificially high: although each data point is described by thousands of features, it is assumed that it can be modeled as a function of only a few underlying parameters. Formally, it is assumed that the data points are samples from a low-dimensional manifold embedded in a high-dimensional space. In this paper, we discuss a recently proposed method, known as Maximum Entropy Unfolding (MEU), for learning non-linear structures that characterize high dimensional data. This method represents a new perspective on spectral dimensionality reduction and, joined with the theory of Gaussian Markov random fields, provides a unifying probabilistic approach to spectral dimensionality reduction techniques. Parameter estimation as well as approaches to learning the structure of the GMRF are discussed © 2015 The Authors. Published by Elsevier B.V. Peer-review under responsibility of Spatial Statistics 2015: Emerging Patterns committee.


Introduction
The aim of statistical methods for dimensionality reduction (DR) is to detect and discover low dimensional structures in high dimensional data. Many high-dimensional data in real-world applications can be modeled as data points lying close to a low-dimensional nonlinear manifold. The key observation is that, even if the dimension of the embedding spaces is very high, the intrinsic dimensionality of the data points could be rather limited. Traditional DR techniques, such as principal component analysis and multidimensional scaling, usually work well when the data points lie close to a linear (affine) subspace in the input space. They cannot, in general, discover nonlinear structures embedded in the set of data points.
Different methods have been proposed in literature for nonlinear dimensionality reduction 1 .
In this paper we discuss a spectral nonlinear approach in the framework of Maximum Entropy Unfolding (MEU) 2 . It is shown that a low-dimensional representation of the data can be achieved through a spectral decomposition of the precision matrix of an Intrinsic Gaussian Markov random field (GMRF), which represents an important class of spatial models 3 .
In general, the correlation structure of a GMRF model is hard to determine (except numerically), but the inverse correlations can be directly specified. On a lattice, this gives most, but not all, elements of the inverse dispersion matrix which is required for Gaussian maximum likelihood estimation of the proposed model. The Markov property of a GMRF makes it possible to utilize numerical methods for sparse matrices to construct fast algorithms for sampling and evaluation of the likelihood. Based on the theory of graphical models, we discuss both parameter estimation and possible approaches to retrieve the underlying structure (graph) of a GMRF.
The paper is organized as follows. Section 2 introduces the dimensionality reduction problem, while Section 3 provides details on MEU. In Section 4 we discuss parameter estimation based on Maximum Pseudolikelihood which, resulting in a least square estimator, favours the use of algorithms, such as Elastic Net and Lasso, to learn the graph structure of the Markov random field. Finally, Section 5 concludes the paper by discussing a Bayesian approach to parameter estimation.

The Manifold learning problem
Dimensionality reduction is the transformation of high-dimensional data into a meaningful representation of reduced dimensionality. Ideally, the reduced representation should have a dimensionality that corresponds to the intrinsic dimensionality of the data, defined as the minimum number of parameters needed to account for the observed properties of the data 4 .
Specifically, let be a data matrix of input coordinates consisting in data-vectors, , with . We further assume that the data matrix has intrinsic dimensionality , with . In this context, the aim of dimensionality reduction methods is that of finding a matrix of the coordinates in the reduced space, or reconstructed embedding coordinate matrix, with , while retaining the geometry of the data as much as possible. The assumption that the observed data have intrinsic dimensionality , implies that the data points lie on (or near) a manifold , with dimensionality , that is embedded in the -dimensional input space . To retrieve a faithful low dimensional representation of the data, which accounts for their intrinsic geometry, it is, then, important to determine the manifold metric.
Formally, the dimensionality reduction problem can be stated as follows: find a mapping , where is an affine space, called feature space of dimension , such that the reconstructed embedding coordinates represent the metric structure in well. In other words, given points that lie on a -dimensional manifold that can be described by a single coordinate chart , find , where 5 . Solving this problem is referred to as manifold learning, since we try to learn the manifold structure starting from a sample of data points. Hence, assuming the map to be non-linear allows capturing the nonlinearity of the underlying structure.

Maximum entropy unfolding
In recent years, several spectral nonlinear methods have been developed to address the manifold learning problem, including Isomap 6 , Locally Linear Embedding 7 , Laplacian eigenmaps 8 , Kernel Principal Component Analysis 9 and Maximum Variance Unfolding 10 .
Specifically, the spectral approach to dimensionality reduction involves the eigendecomposition of a similarity matrix, also know as the Gram matrix, whose principal eigenvectors are extracted in order to retrieve a lower dimensional representation of the high dimensional data.
These nonlinear dimensionality reduction techniques are closely related and can be unified in a classical multidimensional scaling (CMDS) 11 perspective, where the Gram matrix is obtained from the pairwise squared Euclidean distances between the data points.
Starting form the perspective of CMDS, Lawrence 2 shows that the main difference between the nonlinear spectral DR techniques is in the distance matrices they use. Moreover, by following a different approach for constructing the distance matrix, Lawrence 2 proposes a new method, called the Maximum Entropy Unfolding, whose main feature is that of providing a probabilistic model. Exploiting the maximum entropy formalism 12 , MEU specifies the probability density by a free form maximization of the entropy subject to constraints imposed on the expectations of the squared distances between two neighboring data points sampled from the model. For any two samples, it thus holds where represents the set of neighboring points of , and is the squared Euclidean distance between and In order to maximize the entropy in a continuous system, Lawrence 2 proposes to maximize the negative Kullback-Leibler divergence (KLD), or relative entropy, between a base density, and the density of interest, ) where is defined as a very broad spherical Gaussian density with covariance , with typically assumed to be close to zero. By writing the set of constraints as where are the Lagrange multipliers, it can be shown that the probability distribution corresponds to a zero-mean Gaussian Markov random field (GMRF) ( 2 ) with precision matrix Note that, in order to capture the underlying structure of the manifold, is defined as the Laplacian matrix, with off diagonal elements given by if is neighbor of and 0 otherwise, and its diagonal elements defined as . Hence, is symmetric and constrained to have a null space in the constant vector, , where is the -dimensional vector of ones. Equation (2) enlightens that MEU implies independence across data features, which is due to the imposed constraints, and that the parameters, , can be estimated through maximum likelihood. In fact, the gradient of each Lagrange multiplier is given by 2 which can be evaluated by computing the expectation of the squared distance as (3) where is the th entry of the kernel (i.e. covariance matrix) . Equation (3) represents a scaled version of the standard transformation between distances and similarities. This relationship arises naturally in the probabilistic model since, in general, every GRF has an associated interpoint distance matrix and it is this matrix that is used in CMDS 2 . The parameter ensures that the precision matrix is positive definite. This implies that, knowing that the Laplacian has a null space in the constant vector, then . This reflects an insensitivity of the covariance matrix to the data mean, and this in turn arises because that information is lost when we specify the expectation constraints only through interpoint distances. In practice, is always centred before its eigenvectors are extracted, , where the centering matrix is , with a identity matrix. Once the maximum likelihood solution is recovered, the data can be visualized by looking at the eigenvectors of the Gram matrix .

Parameter estimation
Optimal parameter estimation techniques, such as maximum likelihood (ML) estimation, can be computationally expensive, especially if noisy data are available and the noise is non-Gaussian. Also, there are many practical applications when fast solutions are required, even at the expense of a non-optimal estimate.
Hence, in light of equation (2), which tells us that our model is a GMRF, and according to graphical models theory, we propose a different approach for parameters estimation based on the maximization of the pseudolikelihood (PL) subject to constraints. Maximum pseudolikelihood (MPL) estimation provides an alternative, quick, and often reasonably efficient method of parameter estimation. Specifically, we discuss constrained MPL estimation to help deal with the problem of estimating Gaussian intrinsic autoregressions for which the model parameters are constrained to be on the boundary of the valid region.
In accordance with earlier discussion, we assume henceforth that each data point has a limited number of neighbors and that the conditional distribution of is fully specified in terms of a vector consisting of a few unknown parameters. Given the data, we shall use to denote the conditional probability of observing , given all other values. The primary objective is to obtain a reasonable estimate of . The most naive approach to the estimation of the unknown parameters in the terms would be to take the vector which maximizes the quantity ( 4 ) with respect to . From a graphical model perspective, this implies the specification of a family of probability distributions defined in terms of directed or undirected graphs. The nodes in the graph are identified with random variables, and joint probability distributions are defined by taking products over functions defined on connected subsets of nodes 13 .
Although maximum pseudolikelihood estimation is intended to have fairly widespread applicability, it is of special interest in the Gaussian case in which we assume that (5) where is the conditional variance. According to Rue and Held 14 , in the case of GRFs, the connected subsets of nodes are specified through the elements of the precision matrix, . The pairwise conditional independence properties of is contained solely in the covariance matrix and for detecting conditional independences, one must investigate : each vertex represents a data point, and an undirected edge connects two data points, , if the corresponding element of the precision matrix, , is non-zero. In order to analyse the nonzero pattern of the precision matrix , and hence, the connectivity of the underlying graph we follow an approach based on the Cholesky decomposition of : where is a lower triangular matrix and represent a weighted adjacency matrix from a directed acyclic graph -DAG. This implies that when constructing the neighborhood, the triangular form for this matrix can be achieved by first imposing an ordering on the data points. Then, when seeking the nearest neighbors for , we only consider a candidate data point . Applying recursively this procedure one vertex at time, we obtain that the -th row of provides an alternative parametrization of the conditional distribution of . These conditional statements are described by a linear recursive system where the zero pattern for the regression coefficients is the same as in the concentration matrix. Hence, in light of these results and from equations (4) and (5), it can be shown that the estimated 's are the result of the solution of a system of recursive regressions ( 6 ) where is a indicator matrix which selects the neighbors of . The 's are strictly related to the coefficients 's of the precision matrix through MPL thus reduces to the ordinary method of least square which, for a design matrix , gives ( 7 ) and ( 8 ) is the estimate of the conditional variance.

Learning the graph structure
In the MEU framework, an important part of the model specification is the choice of the structure of the GMRF or, equivalently, of the precision matrix , which has the form of the Laplacian matrix.
As seen in the previous section, we parameterize the GMRF according to the Cholesky decomposition of , and, in order to constrain the Laplacian matrix to be positive semidefinite we may have to guarantee that assume that , which satisfy . If we force and set the diagonal elements , we will have a Laplacian matrix which is positive semidefinite without need of any further constraint on . One possibility to force a sparse structure in is to work with a nearest neighbor rule, which assigns the same number of (nearest) neighbors to each data point.
The constraint can be imposed by ensuring that sum of the off diagonal elements from each column is equal to one, i.e.
. Under general linear constraints we thus have that the negative pseudolikelihood can be reformulated as ( 9 ) where is a unit matrix which defines the linear system of the parameters which must be solved accordingly to the constraints in the -vector , and is a -vector of Lagrangian multipliers. It is easy to show that the minimum of equation (9) is given by ( 1 0 ) where . Equation (10) provides a constrained MPL estimator for the regression weights, which guarantees invariance to all the similarity transformations of the data.

The Elastic Net algorithm
An alternative procedure to learn the graph structure in is that of using the Elastic Net (EN) algorithm. This is a regularization and variable selection method that adds the penalty on the L2-norm to the L1 penalty of the Lasso. The Elastic Net is particularly useful when the number of predictors ( ) is much bigger than the number of observations ( ), which is a case found in many examples. For any fixed non-negative and the elastic net criterion is defined as (12) where is the data matrix without the th input coordinate, and . We refer to Hastie et al. 15 for know results on EN and its practical implementation.

Conclusion
In this paper, we considered MEU approach 2 that represents a novel approach to dimensionality reduction based on Gaussian random fields. The main advantage of MEU is that it provides a probabilistic model and a unifying perspective to spectral DR methods.
Given that we are required to estimate the parameters of a GMRF, we showed that, instead of applying the ML, one may use the PL approximation for obtaining these estimates. This is due to the conditional independence properties of GRMFs. We showed that the optimization of the PL turns out to be equivalent to solve a system of recursive equations, and hence we derived the maximum PL estimator through the least square method. Hence, PL allows for explicit estimation of the parameters through the regression formalism.
In order to guarantee the constraint concerning the null space of the Laplacian matrix, we derived also a constrained maximum PL estimator.
Moreover, in order to avoid user-defined parameters for the graph building, typical of rules like -nn algorithm, we saw that combining PL with Elastic Net algorithm allows us to perform DR while learning the neighborhood relations directly from the data. Generally, by a good tuning of and , EN can increase both sparsity and accuracy of the final solution.
Developing a Bayesian approach to achieve sparsity will be of interest for future works. One promising approach would be that of considering Bayesian spike and slab priors, in which the sparsity is induced by placing a mixture prior on the regression coefficient. Specifically, we will exploit the results obtained by George and McCulloch 16 , who proposed a stochastic search variable selection (SSVS) algorithm for normal linear regression using Gibbs sampling to search for high posterior probability models. Their approach relies on a mixture of a low and high variance normal prior centered at zero for each of the regression coefficients, with the low variance component corresponding to a predictor being effectively excluded due to the coefficient being close to zero. For each independent regression in Equation (6), we suspect only a subset of the elements of are non-zero. We assume that arises from one of two normal mixture components, depending on a latent variable : where is positive but small s.t. is close to zero when ; is large enough to allow reasonable deviations from zero when . In addition, the prior probability that is: To obtain the normal mixture prior for , George