Gaussian Process Models of Spatial Aggregation Algorithms

1 Introduction

Many important tasks in data mining, scientific computing, and qualitative modeling involve the successive and system­ atic spatial aggregation and redescription of data into higher­ level objects. · For instance, consider the characterization of WCDMA (wideband code-division multiple access) wireless system configurations for a given indoor environment. In noisy channels, the performance goal is to quantitatively as­ sess the relationship between the signal-to-noise ratio (SNR) and the bit error rate (BER) or bit error probability (BEP) of the realized configuration. To improve performance in of­ fice environments (characterized by doorways, walls, cubi­ cles), a common trick used is to incorporate space-time trans­ mit diversity (STID). Instead of a single transmitter antenna,

Chris Bailey-Kellogg

Department of Computer Sciences Purdue University, IN 47907, USA cbk@ cs.purdue.edu

the base station uses two transmitter antennas separated by a small distance. If the signal from one of the antennas is weak, the signal from another is likely to be high, and the over­ all performance is expected to improve. In this application, it is important to assess how the power-imbalance between the two branches impacts the BEP of the simulated system, across a range of SNRs (see Fig. 1; [Verstak et al., 2002]).

Characterizing the performance of WCDMA systems re­ quires the identification of multi-level spatial aggregates in the high-dimensional configuration spaces of wireless sys­ tems. The lowest level (input) contains individual Monte Carlo simulation runs providing unbiased estimates of BEPs. This space is high-dimensional (e.g. 2:': 10), owing to the mul­ titude of wireless system parameters (e.g. channel models, fading characteristics, coding configurations, and hardware controls). Wireless design engineers prefer to work in at most two or three dimensions (e.g. to study the effect of power im­ balance on system performance) for ease of tunability and de­ ployment. The next level of spatial aggregation thus contains buckets which aggregate data in terms of two dimensions, us­ ing various consistency constraints and design specifications. Finally, the third level aggregates buckets into regions of constrained shape; the shape of the regions illustrates the nature of joint influence of the two selected configuration parame­ ters on performance. Specific region attributes, such as width, provide estimates for the thresholds of sensitivity of configu­ rations to variations in parameter values.

The results of such mining are important for both qualita­ tive and quantitative analysis. For instance, when the average SNRs of the two branches are equal, the BEP is minimal and the width of the mined region in Fig. 1 depicts the largest ac­ ceptable power imbalance (in this case, approximately 12dB). However, the width is not uniform and the region is narrower for smaller values of the SNRs. The qualitative result is that system designs situated in the lower left corner of the con­ figuration space are more sensitive to power imbalance in the two branches.

Each input data point captures the results of a wireless sys­ tem simulation which takes hours or even days (the simula­ tions in Fig. 1 were conducted on a 200-node Beowulf cluster of workstations). Thus it is imperative that we focus data col­ lection in only those regions that are most important to sup­ port our data mining objective, viz. to qualitatively assess the performance in configuration spaces. This requires that we model the functioning of the data mining algorithm, in order to optimize sample selection for utility of anticipated results. Modeling data mining algorithms in this manner is useful for closing the loop, characterizing the effects of the data mining algorithm's parameters, and improving our understanding of how variations in data fields percolate up through the layers. A particularly interesting application is to use such modeled structures to design information-theoretic measures for eval­ uating experimental designs [MacKay, 1992] and for active data selection [Cohn et al., 1996; Denzler and Brown, 2002].

In order to address these goals, this paper develops a novel Gaussian process approach to modeling algorithms that mine spatial aggregates. We first overview the Spatial Aggregation mechanism for spatial data mining and the Gaussian process approach to Bayesian modeling. We then show how to integrate the two approaches in order to achieve our goal of prob­ abilistically modeling spatial data mining algorithms. We il­ lustrate this ability within the context of identifying pockets underlying the gradient in a field - an application that cap­ tures many of the interesting characteristics of more complex studies like the wireless application.

2 Spatial Aggregation

The Spatial Aggregation Language (SAL) [Bailey-Kellogg et al., 1996; Yip and Zhao, 1996], provides a set of operators and data types, parameterized by domain-specific knowledge, for uncovering and manipulating multi-layer geometric and topological structures in spatially distributed data. SAL ap­ plications construct increasingly abstract descriptions of the input data by utilizing knowledge of physical properties such as continuity and locality, expressed with the vocabulary of metrics, adjacency relations, and equivalence predicates. To understand the SAL approach (see Fig. 2), consider a SAL program for analyzing flows in a vector field (e.g. wind ve­ locity or temperature gradient).

In the first level, the goal is to group input vectors (a) into paths so that each sample point has at most one predeces­ sor and at most one successor. SAL breaks the process into two key steps, one capturing locality in the domain space (i.e. sample location), and the other capturing similarity in the feature space (i.e. vector direction). A neighborhood graph aggregates objects with a specified adjacency predicate ex­ pressing the notion of locality appropriate for a given domain. As shown (b), a sample point's neighbors include all other points within some specified radius r. Feature comparison then must consider only neighbors in this graph, thereby ex­ ploiting physical knowledge to gain computational efficiency while maintaining correctness. Here we break feature com­ parison into a sequence of predicates and graph operations. In particular, we first filter the graph (c), applying a predicate that keeps only those edges whose direction is similar enough

(within some angle tolerance B) to the directions of the vec­ tors at the endpoints. The remaining graph has some "junc­ tion" points where vector direction suggests multiple possible neighbors, and the most appropriate path extension from the point must be chosen. A similarity metric sums the distance between the junction and a neighbor, weighted by a constant d, and the difference in vector direction at the junction and the neighbor. The most similar neighbor for the junction is selected ((d) and (e), for successor and predecessor junctions, respectively).

The remaining graph edges are collected and redescribed as more abstract streamline curve objects (f), for the second level of analysis. Again, computation is localized so that only neighboring streamlines are compared. The neighbor­ hood graph here (not shown) uses an adjacency predicate that declares streamlines neighbors if their constituent points were in the first level. It is then straightforward to identify conver­ gent flows (g) with an equivalence predicate that tests when constituent points form a junction in the graph in (c). If de­ sired, these flow bundles can be abstracted and analyzed at an even higher level.

SAL's uniform spatial reasoning mechanism, instantiated with appropriate domain knowledge, has proved success­ ful in applications ranging from decentralized control de­ sign [Bailey-Kellogg and Zhao, 1999; 2001], to weather data analysis [Huang and Zhao, 1999], to analysis of diffusion­ reaction morphogenesis [Ord6fiez and Zhao, 2000]. Recent work has focused on optimizing sample selection for applica­ tions where data collection is expensive, including identifying flows in multi-dimensional gradient fields [Bailey-Kellogg and Ramakrishnan, 2001] and analyzing matrix properties via perturbation sampling [Ramakrishnan and Bailey-Kellogg, 2002]. This paper provides the mathematical foundations necessary for the modeling of such SAL programs to support the meta-level reasoning tasks outlined in the introduction.

3 Gaussian Processes

Gaussian processes have become popular in the last few years, especially as a unifying framework for studying mul­ tivariate regression [Rasmussen, 1996], pattern classifica­ tion [Williams and Barber, 1998], and hierarchical model­ ing [Menzefricke, 2000]. The underlying idea can be traced back to the geostatistics technique called kriging [Journel and Huijbregts, 1992], named after the South African miner Danie Krige. In kriging, the unknown function to be modeled (e.g., ozone concentration) over a (typically) 2D spatial field is expressed as the realization of a stochastic process. A prior is placed over the function space represented by this stochas­ tic process, by suitably selecting a covariance function. Given measured function values at sample locations, kriging then proceeds to estimate the parameters of the covariance func­ tion (and any others pertaining to the random process). Using such values a prediction of the response variable can then be made for a new sample point, typically using MAP or ML in­ ference. This basic approach is still popular in many tasks of spatial data analysis.

Even though parameters are estimated in this approach, it is important to note that kriging is fundamentally a memorybased technique, since the estimated parameters only describe the underlying covariance function of a stochastic process. Thus, predictions of the response variable for new sample points are conditionally dependent on the measured values and their sample points; by unrolling the effect of the pa­ rameters of the random process, we can directly express this dependency.

Kriging is often motivated as a local modeling technique, capable of approximating or interpolating functions with mul­ tiple local extrema, and generalizes well to applications ex­ hibiting anisotropies and trends. The stochastic prior is also viewed as a mathematically elegant mechanism to impart any available domain knowledge to the modeling technique. In 1989, Sacks et al. [Sacks et a!., 1989] showed how krig­ ing can actually be used to model processes with determin­ istic outcomes, especially in the context of computer exper­ iments. The justification for modeling a deterministic code as a stochastic process is often that even though the re­ sponse variable is deterministic, it may 'resemble the sam­ ple path of a suitably chosen stochastic process' [Sacks et al., 1989]. Alternatively, using a stochastic process prior can be viewed as a Bayesian approach to data analysis [Sivia, 1996], and this is the idea emphasized by most recent com­ puter science research in Gaussian processes [Gibbs, 1997; Rasmussen, 1996]. The stochastic process can be suitably formulated to ensure that the model reproduces the same re­ sponse value for repeated invocations of a given sample input (i.e., absence of random error). For instance, the Gaussian prior can be chosen so that the diagonal entries of the covari­ ance matrix are 1, meaning that the model should interpolate the data points.

In the recent past, Gaussian processes have become popu­ lar in the statistical pattern recognition community [MacKay, 1997] and graphical models literature [Jordan (ed.), 1998]. Neal established the connection between Gaussian processes and neural networks with an infinite number of hidden units [Neal, 1996]. Such relationships allow us to take tradi­ tional learning techniques and re-express them as imposing a particular covariance structure on the joint distribution of in­ puts. For instance, we can take a trained neural network and mine the covariance structure implied by the weights (given mild assumptions such as a Gaussian prior over the weight space). Williams motivates the usefulness of such studies and describes common covariance functions [Williams, 1998].

Williams and Barber [Williams and Barber, 1998] describe how the Gaussian process framework can be extended to clas­ sification, where the modeled variable is categorical. Essen­ tially, the idea is to (i) use a logistic function to conduct tradi­ tional Gaussian regression modeling, and (ii) adopt a softmax function to bin the logistic output into a given set of classes. This means that the logistic function uses a "latent variable" as input in its computation, since its values are not provided by the dataset.

4 Gaussian Processes for Spatial Aggregation

SAL programs construct multi-layer spatial aggregates based on specified local adjacency relations, similarity metrics, and consistency checks. We describe here how to capture the qualitative behaviors of such aggregates using Gaussian pro­ cesses. The essence of a Gaussian process is its covariance structure, so we focus on determining covariance structures in a SAL program. For example, in the two-layer SAL pro­ gram of Sec. 2, the parameters (r, (),d) impose a covariance structure by specifying the reach of the neighborhood graph, enforcing the similarity of angles in the vector field, and pe­ nalizing for the distance at decisions involving junctions.

4.1 Covariance Structure

We now describe how to model the covariance structure of a given SAL program. We give the mathematical framework for the case of mining a lD field to determine if there is are­ versal of gradient as we move along the spatial dimension, but essentially the same machinery applies to two and higher di­ mensional spaces. The basic problem is one of classifying 1 D points to determine the qualitative structure of same-direction flows. Fig. 3 (top) depicts the given input field along the x di­ mension. As is shown, the field consists of unit vectors with different orientation. The Gaussian process approach is first to model an underlying regressed variable and then to use a logistic or softmax function to bin the output into classes. In our application, the regressed variable represents the gradient and can be simply summarized as the angle of unit vector ori­ entation yin Fig. 3 (top). In other applications, the regressed variable could be an unobserved 'latent' variable. In either case, it is modeled as a function f of the input x.

First, assume f to be a Gaussian process on x, meaning that the conditional probability distribution of y given a value of x is a Gaussian. For instance, Fig. 3 (bottom) depicts measured values of y superposed with distributions of y at two unseen points. A covariance structure among the y values could, for example, capture the intuition that adjacent values of y should agree more than distant values. The goal of modeling is to determine the extent and stringency of this neighborhood re­ lation - one of the defining parameters of a SAL program. Specifically, we posit a process such as:

The idea then is to estimate a model f' of the same form as f, on the basis of a given set of k observations {(x1,yl), (x2, y2), ... , (xk, yk)} . A typical choice for Z in J' is a random process with zero mean and covariance cr 2 R, where scalar cr 2 is the estimated variance and R is a matrix that captures the correlation between the inputs (i.e., the given locations). Notice that even though the input is one dimen­ sion, the size of R depends on the number of locations for which gradient measurements are available. The above model for f also includes the constant term a:; this can be estimated based on the k observations, or we can substitute more com­ plex terms (e.g. linear), or even omit it altogether.

The functional form of R (including its parameterization) in effect defines the stochastic process and must be carefully chosen to reflect the underlying data's fidelity or any domain­ specific assumptions about local variation. The parameters of the process are then estimated using multidimensional opti­ mization involving a suitable objective function. For instance, given the following form for R:

the problem reduces to estimating p from the given data. Notice that this formulation for R implicitly enforces that the model exactly interpolate the given data points, since R(xi, xi) = 1. A common objective function for estimating pis to minimize the mean squared error (MSE), E{(f'-!) 2 } between f and J'. The p that minimizes MSE is given by the solution to the optimization problem:

where R is the symmetric correlation matrix formed from R. For a new sample point Xk+l> a prediction for the regressed variable is given by:

where r is the correlation vector between the response at Xk+I and all the other k points (derived from R), Ik is the identity vector of dimension k, and 0: is the estimate of a: given by:

The variance in the estimate is given by:

In this case, the optimization is one-dimensional due to the presence of the single parameter p. With a different param­ eterization, we will employ multi-dimensional optimization over the entire set of hyperparameters. When dimensional­ ity is large, the hyperparameters are estimated using MCMC methods. Once such a modeling is complete, as discussed in the previous section, we can relate a categorical class variable to y using softmax functions. For instance, the reversal of the gradient in Fig. 3 can be captured by first using the Gaussian process model to make predictions of the gradient at untested points and then determining if (and where) a zero crossing occurs.

The above equations extend naturally to a 2D case such as that described in Sec. 2. The covariance prior has to be suit­ ably parameterized and we also have the option of taking into account any interactions between the two dimensions (both linear and nonlinear).

4.2 Modeling Many Layers

When SAL programs consist of many layers, we need to de­ velop a sequence of Gaussian process models, each with a suitable covariance function, which can then be superposed to yield a composite covariance function. Recall that while one could simply assess the covariance of the output field for sample values of the parameters and a given input field, the real purpose of a Gaussian process model is to express the co­ variance of the output as a function of the characteristics of the input. This is the key property that allows reasoning about closing the loop and selecting optimal samples. In addition, Gaussian process models help capture the randomness inher­ ent in some of SAL's computations, e.g. non-determinism in labeling, and variations due to how ties are broken for aggre­ gation purposes. Refer again to Sec. 2 for an example of the types of operations that the covariance model must capture.

At the very bottom of the hierarchy is the input data field. For applications characterized by expensive data collection (as in the introduction), it can be advantageous to start with a sparse set of sample data. The Gaussian modeling approach to regression is ideal for creating surrogate representations of data fields from such a sparse dataset. That is, given a sparse set of samples, interpolate a dense field satisfying those val­ ues and incorporating any appropriate domain knowledge, as discussed above regarding kriging. Such surrogate functions can then been used as the starting points for qualitative anal­ ysis [Bailey-Kellogg and Ramakrishnan, 2001].

The operators in a SAL level deal with both locality (which object locations are close to which other ones, as encapsu­ lated in a neighborhood graph) and similarity (which object features are close to which other ones, as encapsulated in met­ rics and predicates). For instance, in the example of Sec. 2, two points are assigned to the same pocket if they are spa­ tially proximate and their flows converge. Here the Gaussian process is classification (or more generally, density estima­ tion). A popular covariance structure for an n-dimensional input field captures locality:

where the expression relates the function values at positions x <k) and x<ll. If TJ E [0, 2], then the covariance function will be positive definite, satisfying the normalization constraints of a posteriori inference.

To see how to capture similarity, consider when two sam­ ple locations are classified into the same trajectory in Sec. 2. In addition to being spatially proximate (as inferred by SAL's neighborhood calculations), the underlying vector fields must also be similar in direction. Expressing the covariance in terms of position alone can cause the resulting estimated hy­ perparameters to be misleading or difficult to interpret, as their effect is confounded with the underlying vector field. One solution is to artificially inflate the dimensionality, so that position and direction together describe the data. Besides increasing the dimensionality, this approach spells trouble for estimation using MCMC methods since significant portions of the sample space will remain unsampled and it would be difficult to assess their effects on the minimized functional. An alternative solution is to use the fact that the vector field is itself a surrogate and add a term to the covariance outside the above structure, capturing the contribution due to similar­ ity in the vector field. We place a Gamma prior on this term with a shape parameter that ensures that its role is secondary to the covariance structure on position (directional similar­ ity alone is not enough for high covariance at the output; the sample locations also must be spatially proximate). This is recognized in the statistics community as a hierarchical prior and described in detail in [Neal, 1997].

5 Experimental Results

In order to test our approach, we studied de Boor's pocket function (see Fig. 4):

where X is then-dimensional point (x1, x2, · · ·, Xn) at which the pocket function p is evaluated, I is the identity n-vector, and II · II is the £2 norm. This function exploits the fact that the volume of a high dimensional cube is concentrated in its comers and pis designed so that it has a "dip" in each comer. It embodies many aspects of datasets like those encountered in the wireless simulation study, including multiple local ex­ trema, non-systematic variation in the location of the pockets, and regional variation. The pocket function is also important as a benchmark for high-dimensional data exploration, where the goal is to identify the most interesting regions of the de­ sign space without necessarily conducting a (costly) global optimization over the entire design space. Data mining pro­ grams are hence required to identify the most promising re­ gions using as few function evaluations as possible.

A SAL program to identify the number of pockets starts with some samples ofthe pocket function. The lowest level of modeling involves a kriging interpolation over a uniform grid (we chose size 13n for testing). Then the approach of Sec. 2 is applied to the gradient vector field of this scalar field: the sec­ ond level bundles points into curves, and the third aggregates these into flows. Each convergent flow represents one pocket. One could mine the covariance structures for each layer sep­ arately; we unfold these mappings here to obtain a single co­ variance structure summarizing all three layers. This is be­ cause the structure (esp. the contributions of each dimension) is easiest to interpret in terms of the original spatial field.

We conducted a parameter sweep over (r, (},d) as:

r E {1, h, 1.5, J3, 2}

0 E {0.7, 0.8, 0.85, 0.9, 0.95}

d E {O.Gl, 0.02, 0.03, 0.04, 0.05}

and used Neal's Bayesian modeling software [Neal, 1997] to construct Gaussian process classifiers for the flow classes. Covariance contributions in the p terms (Eq. 7) from both the dimensions was estimated using hybrid Monte Carlo (ag­ gressive schemes to evolve the system state by adding higher order terms). This procedure uses a leapfrog scheme to sup­ press random walk behavior by selective iteration between Gibbs sampling scans and latent value updates.

Our results indicated a strong positive correlation between the x and y covariance contributions, bringing out the sym­ metry in the underlying SAL computations. The number of pockets mined was constant across the values of d (other pa­ rameters fixed), and one of the goals of our study was to de­ termine if this negligible effect of d is captured in the covari­ ance structures. (The effect of d would actually be more pro­ nounced in other spatial fields but not so much in the pocket function due to the inherent symmetry.) Parameters r and (} produced the most variation in the covariance contributions with(}= 0.95 causing an abrupt jump in the number of mined pockets. This is due to the rather stringent limit imposed on vector similarity arising from the nonlinearity of the cosine metric. Fig. 5 summarizes the results for a 2D pocket func­ tion, where we have averaged the covariance contributions across all ds, for given r and(}.

As the number of pockets increases (Fig. 5(a)), the co­ variance contributions increase (Fig. 5(b,c)) approximately quadratically. In other words, as the underlying latent func­ tion varies rapidly along the given dimensions, we cannot stray "too far" away from a given sample point when mak­ ing predictions at test points. The reciprocal of the covari­ ance scale term is often referred to the characteristic length of a dimension. This gives an estimate of "how far" a given dimension's effect holds. When only four pockets are mined, the characteristic length is about 1, meaning pockets occupy a width of 1 x 1 (exactly one fourth of the total space). As more pockets are mined, the characteristic length drops to about 0.4. It is also interesting to note that the abrupt jump in the number of pockets forB = 0.95 is reflected by a similar in­ crease in the covariance contributions for this value. Essen­ tially, vector and edge directions have to be so similar that few long "runs" can be aggregated as streamlines. This brings out the capability of the Gaussian process approach to capture the essentials of a spatial aggregation algorithm.

6 Discussion

This paper has demonstrated a novel Gaussian process ap­ proach to modeling the qualitative behavior of SAL pro­ grams; in contrast to much of the literature where Gaus­ sian processes are used for pattern classification and regres­ sion [Rasmussen, 1996; Gibbs, 1997], our work takes existing data mining algorithms and recasts them in terms of Gaussian priors. To the best of our knowledge, this is the first study to completely model a qualitative data mining algorithm in terms of a process framework, summarizing the transforma­ tion from data to higher-level aggregates. This is an impor­ tant step in firmly establishing a probabilistic basis for spatial aggregation computations. The modeling undertaken here, while expensive, is justifiable for studies such as the wireless system characterization described in the introduction.

There are several immediate gains from the work presented here; due to space limitations we only mention them briefly. First, the Gaussian process model can characterize experi­ mental design criteria such as entropy as a functional w.r.t. the input space, allowing us to use the mined covariance struc­ ture to focus sampling at the most informative points (e.g., see [Bailey-Kellogg and Ramakrishnan, 2001]). It is impor­ tant to note, however, that the approach taken in [Bailey­ Kellogg and Ramakrishnan, 2001] only addresses the lowest levels of a hierarchy and is unable to reason about higher­ level, more abstract processes of redescription and aggrega­ tion as is done here. Second, Gaussian process models al­ low us to study the effects of different SAL parameters for a given class of datasets, e.g. the inference above of the neg­ ligible role of d in the mining process. Finally, it allows us to take algorithms that function in differing ways (and using different sets of parameters) and places them on a common footing, namely the language of covariance structures. This means that we can reason about the applicability of different algorithms by studying the constraints they impose on spatial locality and field similarity.

Gaussian processes have recently been linked to kernel­ based methods, as used in support vector machines [Cristian­ ini and Shawe-Taylor, 2000]; we intend to explore this con­ nection in future work. Kernel-based methods are attractive in their promise to overcome the curse of dimensionality by the use of nonlinear projections, a facet that is of critical im­ portance for mining data from large parameter sweeps. As the need for data mining in computational science gains promi­ nence, process models will be crucial to achieve effective uti­ lization of data for mining purposes.

Acknowledgements

The authors thank Feng Zhao and Layne Watson for helpful comments. This work is supported by US NSF grants EIA9974956, EIA-9984317, and EIA-0103660.