Is Nonparametric Learning Practical in Very High Dimensional Spaces?

1 Introduction

Nonparametric learning algorithms are designed for learning problems where relatively little information is available a priori about the what type of model struc­ ture is required, or which set of variable inputs are im­ portant. Such problems are common in many domains of computational intelligence. As a simple example, con­ sider the problem of autonomous navigation in a clut­ tered environment . As researchers we must decide which features (inputs) in the environment are important , and

Peter D. Lawrence

Dept. of Elect. and Comp. Eng. University of British Columbia Vancouver, B.C., V6T 1Z4, Canada

peterl@ee.ubc.ca

http :/ /www.ee.ubc.ca/~peterl

what kind of computational models are required. Non­ parametric learning techniques are of potential value to us because they are designed to define their own model structure, and to determine which set of inputs are most relevant to effectively predict the task output. There are a number of effective non-parametric learning techniques in the literature, including C4.5 [Quinlan, 1993], CART [Breiman eta/., 1984], and MARS [Friedman, 1991] (See [Michie et a/., 1994] for a more complete list).

Almost without exception , contemporary implemen­ tations of nonparametric learning algorithms have two common properties. First, some form of variable selec­ tion is done to determine which of all possible inputs best predict the output. And second, they utilize some type of "greedy" search procedure which attempts to pick the best model structure, from some set of possi­ ble structures . Such techniques are feasible when there are relatively few input variables ( < 100). However , neither of these two algorithmic properties are practi­ cal when large numbers of potential inputs are available. Consider the learning example where there are 100 po­ tential inputs, and assume that not all of these inputs are necessary to predict the output. In order to pick the best 50 of these 100 inputs , one would need to construct ( 1 5 ° 0 ° ) =: 1 x 10 29 models , and then evaluate them all to see which set of 50 inputs best predict the output. Such an exhaustive search is currently not feasible , and is unlikely to be in the near future. Similarly, if we use some greedy search procedure to determine which of all possible model structures is the best , we would further increase our search space by the total number of possible model structures .

Certainly, the flexibility of non parametric learning rec­ ommends it as a useful tool for constructing intelligent agents. This has been clearly demonstrated in numer­ ous publications [Michie eta/., 1994] . However, in almost all applications of nonparametric learning, the number of inputs have been relatively low ( < 100), and the number of learning examples have been relatively few ( < 20, 000). Thus, given that many challenges faced by the field of Computational Intelligence in building in­ telligent agents, involve determining mappings between numerous (thousands or more) and varied sensor inputs and complex and flexible action sequences requiring tens of thousands of training examples, we must ask whether nonparametric learning can be usefully applied to such problems? Given our argument in the previous para­ graph, the initial response to this question is no: the computational complexity associated with applying ex­ isting algorithms to such large problems makes them im­ practical. In fact, all currently published nonparamet­ ric learning algorithms, utilize search strategies which are impractical in high dimensional spaces. Even algo­ rithms such as projection pursuit [Friedman and Stuet­ zle, 1981], which do not perform traditional exhaustive searches, build models using all input variables simulta­ neously, making them impractical for very high dimen­ sional problems.

In an effort to develop methods to address the high dimensional mappings required for intelligent agents, we challenge the need for variable selection and greedy func­ tion search in nonparametric learning. Instead, we pro­ pose an algorithm which has two basic characteristics: first, very little computational effort is spent on variable selection; and second, the model selection strategy does not explode in computational complexity, as the dimen­ sion of the problem increases . Our intent is to show that a model having these properties is capable of building practical, effective nonparametric models, for both small low dimensional learning problems, as well as larger, very high dimensional problems .

The basic premise behind the proposed algorithm is that very high dimensional regression can be done us­ ing a finite number of low dimensional structural units , which are added one at a time to the regression func­ tion. These structural units are by necessity low dimen­ sional, because high dimensional structures suffer from the curse of dimensionality (see [Friedman, 1994] for a detailed discussion) . The inputs to new structural units can be the outputs of previously added units as well as input variables. However, both the choice of which struc­ tural unit to add next, and which inputs will act on these structural units, must be done using minimal computa­ tional effort for the algorithm to be practical in high di­ mensional spaces.

In the current paper, we implement and evaluate the simplest algorithmic structure which is consistent with the ideas outlined above. The proposed algorithm has a very simple structure, using only one type of structural unit , and variables are added to the model in a random order. The surprising conclusion of this paper is that this simple, in essence random structure, produces highly effective nonparametric models.

Much of the work in computational learning has been directed towards the classification problem. However, many real problems in AI involve controlling agents which manipulate the world using continuous values con­ trol signals. Hence, in this paper we have chosen to directly address the continuous valued regression prob­ lem, instead of the discrete valued classification problem. However, as discussed in [Friedman, 1994], classification problems can be solved using regression techniques in a number of ways. Therefore the high dimensional re­ gression techniques studied in this paper are, at least in theory, applicable to classification problems.

In the remainder of this paper we give a detailed de­ scription of the proposed algorithm, followed by an ex­ tensive numerical evaluation. Our evaluation is done in 3 parts. First we compare the performance of our algo­ rithm to published results on 10 well known regression problems from the literature . This gives us a good mea­ sure of the performance of the proposed algorithm on some standard regression problems. Next we test the al­ gorithm on 15 large very high dimensional , highly non­ linear regression problems . This allows us to test the efficacy of the algorithm on the type of problems we are most interested in addressing. Our final evaluation of the algorithm is on the 10-bit parity problem (or equiv­ alently the 10 input XOR problem) . This problem was chosen to show that the proposed algorithm is capable of solving problems which have an intrinsic dimension (10 inputs) greater than that of its highest dimensional functional unit (2 dimensional). A common criticism of algorithms which build functions using low dimensional structural units is that they cannot model high dimen­ sional XOR functions . In this paper , we demonstrate that this is not true for the proposed algorithm.

2 The Proposed Learning Algorithm

Let RL(x1, .. . , xN) be anN dimensional regression func­ tion constructed using an L level cascade of 2 dimen­ sional functions: y = RL(x) estimates the dependent variable y given independent variables x = ( x1, ... , x N). Then RL(x) has the following form:

where L defines the number of levels in the cascade or, equivalently, the number of 2 dimensional functions g1( u, v) constructed; O't are real valued scaling factors; and, the subscripts ko , ... , kL E {1 , ... , N} serve to identify the independent variables (x1 , ... , XN)· For the purposes of this paper, the 2 dimensional functions 91 ( u, v) are represented as 3rd order 2 dimensional polynomials:

where the coefficients a1 ... , a10 are defined as each level is constructed. We postulate that this choice for repre­ senting 91 ( u, v) is not unique (initial experiments indi­ cate that higher order polynomials do not significantly affect results) , however , as demonstrated in Section 3, it is an effective one.

We assume that ML learning input/output examples (x;, y;) , fori= 1, ... , ML have been divided into 2 sets: a training set rT containing MT input/output pairs , and a validation set fv containing Mv input/output pairs. In Section 3 we describe how the training and validation sets are obtained. There are 3 preset learning param­ eters, symbolized by depth1 , depth2 , and c, which are used during the construction of the regression function. For all experimental results presented, these learning pa­ rameters are set to the same theoretically motivated and data independent values. The function of each learning parameter is explained in the following algorithmic de­ scription, and is further elaborated on below. A single cascade of gr ( ·) functions is constructed in a series of sections. When the current section i Rj = ~{=i ar · gr(-) of the cascade can no longer reduce mean squared er­ ror , learning outputs y are replaced by residual errors due to i Rj and a new section is begun starting from the last error reducing function gj (-) . Stated in detail, the construction of the regression function, RL(x) ( equa­ tion (1)), proceeds as follows:

STEP 1: Initialize algorithm: Initialize the 2 coun­ ters i = 1 and p = 1, where i and pare used as follows: i indexes the function g;(-) at the start of the current section; p is a subscript indicating that section p of the cascade will be constructed in step 2. Throughout this algorithm, the subscripts ko, ... , kL (used to identify the independent variables) are randomly selected, one at at time, from the set {1, ... , N} , without replacement. When this set is empty, it is re-initialized to {1, ... , N}, and the process of selecting these subscripts continues.

STEP 2: Construct new section: Construct, in or­ der, the functions g;(-), g;+1(-), ... , gj(-), where j- i > depth1, as follows:

1. Obtain a bootstrap sample set r~ (i.e. MT ran­ domly chosen samples with replacement) from the training set fT.
2. Given r~, use Singular Value Decomposition to determine the coefficients a1, ... , a 10 of equa­ tion (2).

The construction stops with function gj (-) when the change in cascade error on the validation set fv, over the past depth1 + 1 level additions, is less than some small, positive real number t. This is defined as follows:

where i Rq = ~r=i ar · gr(-) is the constructed cascade between levels i and q; and MSE( iRq)rv is the mean squared error of this cascade on the validation set fv . With the addition of each level , the scaling factors a1, for I= i, ... ,j, are recomputed as follows :

where MSE(gr(-))rr is the mean squared error of the single function gr(-) on the training set rT . Hence, ar is simply the inverse of the normalized mean squared error of level I, and its purpose is to give more weight to levels which contribute more significantly to error reduction.

STEP 3: Prune section: We now prune the current section back to the level which has the best mean squared error, as follows. Find the minimum level s E { i - 1, i , ... , j} which gives the least mean squared error of the cascade on the validation set fv , in the following sense:

Delete all functions g,+1 (-), g,+2 (-), ... , gi(-)· Set j = s and, using equation (4), recalculate the scaling factors ar. The current state of the cascade is now given by R,(x) = ~:=l ar · gr(-). Note that if s = i- 1 then all the functions just constructed in step 2 are deleted.

STEP 4: Update learning outputs: Update the out­ put values, y;, of both the validation set r v and train­ ing set rT as follows: V(x, y) E fv u rT => y = y- ~:=i ar · gr(-) . The sets fv and rT now contain the residual approximation errors after the cascade has been constructed to pruned level s. Let Pp be the vari­ ance of the output values, y;, of the validation set fv at the completion of section p. Thus pp is equivalent to the mean squared error of the approximation R, ( x) on the validation set. pp is used in the next step to establish a condition under which the construction of the cascade stops.

STEP 5: Check stopping condition: Go back to STEP 2 , constructing at least depth2+ 1 sections, until the following condition is met :

Therefore, the algorithm will terminate when the mean squared error fails to improve by a factor of f., over the last depthz + 1 sections. Before executing STEP 2 , set i = s + 1 and increment the counter p to initialize the next section. Upon termination, set L = s, the total number of levels constructed.

As indicated above, the learning parameters depth 1, depth2, and f., are all used to define conditions under which the construction of the cascade stops. In par­ ticular , f. defines the relative decrease in error, via the addition of further gr (-) functions , required in order for the cascade to continue to grow. Using f, the learning parameter depth1 defines how many levels of the cas­ cade are used to determine the termination condition of STEP 2 . Similarly depth2 and f. are used in STEP 5 to define how many times STEP 2 is executed before the algorithm is terminated. Theoretically, the best re­ gression function (with respect to least squared error) is obtained when f. = 0, and depth1 and depth2 are very large. However, in order for the algorithm to terminate in a reasonable amount of time, for all of the results re­ ported in this paper we fixed learning parameter values to: depth1 = 25, depth2 = 6, and f. = 0.0001. Hence the algorithm described above is completely automated, requiring no a priori parameter tuning.

3 Results and Discussion

3.1 Low Dimensional Data Sets

In this section we evaluate the efficacy of the proposed algorithm by applying it to 10 well known regression problems from the literature (Table 1).

For the data sets found in [Breiman, 1996] and [Ras­ mussen, 1996], the regression functions were constructed using 10 fold cross-validation : the Learning Data set was divided into 10 approximately equally distributed sets, and then 10 regression functions were constructed using, in turn, 9 of these sets as training sets, and the remain­ ing set as a validation set. For each Test Data point (Test Data was not used during learning), the outputs of the 10 regression functions were averaged to produce the final approximation output, for which error results are reported. To test reproducibility, 100 independent approximations (independent with respect to random sequences of input variables and bootstrap samples as defined in Section 2) were generated using the Learning Data: Table 1 reports the best Test Set error , along with the average and standard deviation (s.d.) over the 100 independent runs.

For the data sets found in [Breiman, 1996], we fol­ lowed the reported experimental setup. For each of the 4 data sets, 100 learning and 100 test sets were created randomly according to the guidelines given in [Breiman, 1996] . For each of the 100 learning sets a regression function was constructed (using 10 fold cross-validation as described above) and evaluated on the corresponding test set. Experimental results for the 100 experiments are &iven in Table 1. The previously published error for the lBreiman, 1996] data refers to the average bagged error reported.

For the data sets found in [Rasmussen, 1996], we re­ port results for the largest learning sets only (we used 80 learning examples of auto price data, 104 of cpu data, 256 of housing data, 192 of mpg data, and 88 of servo data). For each of these 5 fixed data sets, the regression functions were constructed using 10 fold cross valida­ tion as described above. We generated 100 independent approximations in order to determine what effect the stochastic aspect of the proposed algorithm has on its performance. As reported in Section 2, the ordering of the independent variables is random, and the construc­ tion of the 2 dimensional functions g1(-) is done using a (random) bootstrap sample of the training data. The best and average relative mean squared test set error, and its standard deviation, are reported in Table 1. From Table 1, it is evident that although there is some varia­ tion in error performance from run to run, the stochastic effect of the algorithm is mostly negligible. The previ­ ously published error given in Table 1 under the [Ras­ mussen, 1996] data sets is the best reported error (indi­ cated by brackets) of the 5 algorithms evaluated, and was obtained from the graphs presented in the paper. For both the [Breiman, 1996] and [Rasmussen, 1996] data sets, the average learning time ranged from about 1 to 10 minutes per approximation (all learning times reported in this paper are for proposed algorithm running on a Pentium Pro 150 using LINUX). The average size of a single cascade was about 90 K-bytes.

Finally, for the forward dynamics data reported in [Jordan and Jacobs, 1994], our evaluation was done us­ ing the experimental setup described by Jordan and Ja­ cobs for the on-line back-propagation algorithm: learn­ ing was done using 15,000 training examples, and learn­ ing stopped when the error could no longer be reduced on 5000 validation examples. The regression functions for this data consisted of a single cascade per output: due to the large data size, 10 fold cross validation was not necessary. In order to estimate the reproducibil­ ity of proposed algorithm on this data, we constructed 10 independent approximations. The best and average relative error on the validation set (in accordance with [Jordan and Jacobs, 1994]) , and its standard deviation over these 10 independent experiments, is reported in Table 1. The previously published error, shown in Ta­ ble 1, is the best relative error (on the validation data) of the 7 algorithms studied in [Jordan and Jacobs, 1994}. The forward dynamics data consists of 12 inputs and 4 outputs, which requires 4 of our regression functions (1 for each output). The algorithm required about 10 hours of computation to build all 4 cascades, and the average size of each cascade was about 1 M-byte.

From Table 1, it is evident that the proposed algo­ rithm demonstrated as good or better error results on all but 1 (the servo data) of the data sets. However, this result should be interpreted with caution. The specific goal of the proposed algorithm is to build a regression function which best represents the learning data in the mean squared error sense. The regression function con­ tinues to grow in parameter size until the mean squared error can no longer be reduced : this is beneficial if one wants the "best" approximation, but detrimental if one wants a representation of fixed size. Most algorithms referred to in Table 1 are parametric and therefore of fixed size. Two exceptions are MARS and CART. How­ ever, even comparing these to the proposed algorithm should be done with caution: unlike our algorithm, both MARS and CART can be used to analyze the signifi­ cance of various inputs and how they interact with one another.

3.2 High Dimensional Data

In this section, our goal is to study the proposed algo­ rithm when applied to very high dimensional regression data, under various noise conditions. Since no large, high dimensional regression data examples were found in the literature, we applied the proposed algorithm to 15 artificial data sets ranging from 100 to 1600 input di­ mensions . For 5 of these data sets no noise was present , while the remaining 10 data sets had noise. For the first 5 of these "noisy" data sets, noise was present only in the output, while the remaining 5 data sets had both input and output noise. Input noise was generated by

~ 600000 tV

simply allowing only half of all inputs to contribute to the output (thus half of the inputs used in constructing the regression function had no effect on output). Out­ put noise was generated using a normal distribution with the standard deviation selected to give a signal to noise ratio of 3 to 1, thus implying that the true underlying function accounts for 90Y. of the variance in learning and test data.

The regression data was generated using the following function:

where N is the input dimension of the data, while r;C) and s;(·) are continuous 1 dimensional functions which are nonlinear and all different from one another . Equation (7) allows us to generate highly nonlinear data of any dimension, while at the same time controlling the complexity of the data via the appropriate selection of functions r;(-) and s;(-). By data complexity, we are referring to the number of random sample points re­ quired to build sufficient non-parametric models of the data, with respect to least squared error. The more complex the generating function, the more data points are required for non-parametric modeling. For the pur­ poses of this paper , we have chosen functions r;(-) and s; (-) such that 40,000 random sample points are suffi­ cient in order to effectively model the data. For spe­ cific details on how r;(-) and s;(-) are constructed see http:/ jwww. ee. ubc. ca/"'gregg for C code and documen­ tation .

The simulation results for the 12 high dimensional data sets are presented in Table 2. For each data set there are 40,000 learning examples and 40,000 testing examples. Each learning set is further divided into 30,000 training examples and 10,000 validation exam­ ples: these are used to construct one cascade which forms the approximation for that data set. The 3 right most columns of Table 2 contain the test set average relative mean squared error and corresponding standard devia­ tions over 10 independent runs. Relative mean squared error is defined, in the usual sense, as the mean squared error on the test data, divided by the variance of the test data. From Table 2 it is evident that the relative error is small when no noise is present. With the addition of noise the relative error approaches the theoretical limit (due to the 3 to 1 signal to noise ratio) of 0.1. Learning time for each data set with no noise was approximately 7 hours and produced cascades which were between 3 and 6 M-bytes in size. In contrast, learning time for each data set containing noise was about 1 hour and produced ap­ proximations of between 200 and 600 K-bytes in size. It is interesting to observe that the dimension of the data did not affect either learning time or the size of the ap­ proximation. This is clearly demonstrated in Figures 1 through 3, which show an increase in approximation size, as a function of the number of input variables. Note that the no noise data generated approximations which where much larger than those generated by data which had noise. This is because the proposed algorithm stops adding levels to the regression function when error can no longer be reduced on the validation set . This con­ dition occurs much earlier when there is noise in the training data, hence generating much smaller regression functions.

3.3 10 Input XOR Problem

Because the proposed algorithm constructs approxima­ tions 2 dimensions at a time, one may be lead to conclude that it is not capable of solving XOR problems of 3 or more inputs, however, this is not the case. The proposed algorithm was able to solve the 10 input XOR problem in 10 minutes of learning time (the size of the resulting re­ gression function was 40 K-bytes). The theoretical basis for this surprising fact is the subject of ongoing thearetical study, however, we postulate that it is the direct result of our use of bootstrap samples during learning.

4 Conclusion

We have demonstrated that nonparametric learning is practical for very high dimensional learning problems. We did this by proposing a simple nonparametric learn­ ing methodology, which neither requires computation­ ally expensive variable selection, nor an expensive search procedure to find the best model representation. An algorithm based on this methodology was successfully applied to several high dimensional learning problems, which no other currently published algorithm could effec­ tively address. The size of the regression functions pro­ duced by the proposed algorithm depended on the learn­ ing data's complexity, and not on its dimension. Our current ongoing efforts are being directed to a more com­ plete theoretical analysis of nonparametric algorithms having these properties, as well as the application of our methodology to real world agents. In addition, we are interested in exploring the efficacy of the proposed methodology as an analytical tool for exploring the rela­ tive importance of input variables, and for defining con­ fidence intervals for model predictions.