Reduced Heteroscedasticity Linear Regression for Nyström Approximation

1 Introduction

Kernel method is a powerful and important tool for many areas of machine learning, including support vector machines [Cortes and Vapnik, 1995], kernel principle component analysis [Sch¨olkopf et al., 1998], Guassian process [Williams et al., 2002], manifold learning [Talwalkar et al., 2008], etc. These applications usually involve a large kernel matrix that needs up to O(n3 ) computations. Large matrices also arises in other applications like clustering and matrix completion. For large scale problems, the entries of these matrices can be in the order of billions, leading to serious difficulties in computing and storing them. Various approximation techniques have been developed in machine learning to solve the problem, e.g., incomplete Cholesky decomposition [Bach and Jordan, 2005],

∗To whom correspondence should be addressed. This work was done when J. Wu was with the School of Computer Engineering, Nanyang Technological University, Singapore.

Jianxin Wu∗

National Key Laboratory for Novel Software Technology Nanjing University, Nanjing 210023, China

wujx2001@gmail.com

random Fourier features [Rahimi and Recht, 2007] and the Nystr¨om method [Williams and Seeger, 2001].

The Nystr¨om method is a sampling based low-rank matrix approximation method, utilizing the rapidly decaying spectra of kernel matrices. It is usually considered to be originated from the numerical treatment of integral equations and eigendecomposition of matrices. Sampling techniques play a key role in the Nystr¨om approximation. Therefore, various sampling methods have been developed to improve the approximation quality, to name a few, the diagonal sampling and column-norm sampling that apply a weight to each of sampled column with either diagonal element Kii or the L2 norm of the column [Drineas and Mahoney, 2005; Drineas et al., 2006], adaptive sampling that alternates between selecting a set of columns and updating the distribution over all columns [Deshpande et al., 2006]. Among these methods, K-means sampling [Zhang et al., 2008] shows great results on approximation accuracy but at much higher computational cost [Kumar et al., 2012]. In large scale datasets with high dimensional dense features, the time spent in running the K-means algorithm could be significantly longer than the Nystr¨om approximation itself.

This paper presents a novel point of view for the Nystr¨om approximation, leading to an interpretation of success and deficiency of various sampling methods and a new way to improve the approximation quality. To be specific,

• We demonstrate a novel interpretation of the Nystr¨om method. We show that theoretically the Nystr ¨om method can be regraded as a set of point-wise ordinary least square linear regressions of the kernel matrix sharing the same design matrix;
• By employing the linear regression interpretation, we are able to analyze the approximation quality of the Nystr¨om method from a different angle. We show that empirically the heteroscedasticity problem in the ordinary least square linear regression could be used to explain the success and deficiency of various sampling methods, especially the high approximation quality of the K-means sampling method [Zhang et al., 2008];
• We empirically illustrate that positively skewed explanatory variable distributions coincides and may lead to heteroscedasticity. Inspired by this discovery, we can greatly improve the approximation quality of the

Nystr¨om method on many real world datasets with almost no extra computational cost through a transformation function.

The paper is organized as following: The Nystr¨omMethod and its linear regression interpretation is explained in Section 2. We then demonstrate the heterscedasticity, its relation to skewed regressor distribution, and propose our solutions in Section 3. The empirical results on real world datasets are presented in Section 4. Finally, we draw conclusions and discuss future directions of the research on Section 5.

2 The Nystr¨omMethod Revisited

The Nystr ¨om method [Williams and Seeger, 2001] provides a low-rank approximation of a kernel matrix based on eigendecomposition and numerical treatment of integral.

Let K ∈ Rn×n be a kernel matrix with Ki,j = k(xi, xj) and D = {xi}ni=1 . If C = {ci}mi=1 is a subset of D, which we call anchor points in this paper, the Nystr¨om method generates an (up to rank m) approximation ˜K of K by:

where S is an n × m matrix with Si,j = k(xi, cj); X is an m × m matrix with Xi,j = k(ci, cj) and X+ is the MoorePenrose pseudo-inverse of X.

2.1 A Linear Regression Interpretation

The anchor points can be sampled in various fashions, including random sampling [Williams and Seeger, 2001], diagonal sampling [Drineas and Mahoney, 2005], adaptive sampling [Deshpande et al., 2006], K-means sampling [Zhang et al., 2008], K-means sampling with weight [Zhang and Kwok, 2009], etc. For K-means sampling, the anchor points C are chosen as centroids a K-means clustering of D. These studies [Zhang et al., 2008; Kumar et al., 2012] have shown that K-means based sampling provides much better approximation than random sampling at the price of much higher computational cost. In this section, we propose a different interpretation of the Nystr¨om method, which not only explains the success of K-means sampling, but also leads to a computationally efficient improvement of the Nystr¨om method. Different from the classic eigendecomposition perspective, we interpret the Nystr¨om method as a set of point-wise linear regressions (sharing the same design matrix) of the kernel matrix.

For a kernel function k(q, y), if we treat y as constant, k(q, y) becomes a function of the variable q:

The generalized multiple linear regression [Groß, 2003] will estimate the function fy(·) as:

fy(q) is the dependent variable or response variable, e(q) is is a vector of the independent variables or explanatory variables, and ϵ is the error term. We call e(·) the explanatory function. It is a set of basis functions designed to capture the essence of the kernel function k. Choosing the appropriate explanatory function e(·) is very important as it determines the quality of approximation.

Given a set of anchor points C = {ci}mi=1 for learning the regression models, their true dependent variable values are:

It is well known that the ordinary least square (OLS) solution to Equation (3) is:

where X is called the design matrix and

A natural choice of the explanatory function to capture the characteristic of the kernel function k is to set e(q) = s(q), under which the design matrix becomes

Then, fy(q) is approximated by:

Therefore, if we choose the same set of anchor points for all functions fx1 , . . . , fxn we want to approximate, for arbitrary xi, xj ∈ D:

Clearly, the kernel matrix K is approximated by:

in which S is the matrix of similarity functions s(xi), i = 1, . . . , n, so Si,j = k(cj, xi) = k(xi, cj). Therefore ˜K is exactly the same as the Nystr¨om approximation result.

The efficiency of the approximation is ensured by employing the same design matrix (i.e., same set of anchor points) for all examples in D. For each linear regression, we just change the dependent variable's true values (i.e., the values of similarity functions). The generality of the approximation is also guaranteed, as in Equation (3), the y value can be chosen arbitrarily from the input feature space.

This new interpretation allows us to have an alternative explanation for the success and deficiency of the Nystr¨om and improved Nystr¨om approximation [Zhang et al., 2008]. And, more importantly, we gain the freedom to choose the explanatory function e(·) besides using the same form as the similarity function s(·). We will show in the following sections that the fulfillment for assumptions of OLS, especially homoscedasticity, has obvious impact on the quality of the Nystr¨om approximation, and by employing non-symmetric explanatory functions, i.e., set e(·) ̸= s(·), the approximation quality can be greatly improved without extra computational cost.

3 Heterscedasticity, Skewness, Their Adverse Effects on Nystr¨om, and Cure

As stated, the Nystr¨om approximation is essentially a set of point-wise OLS linear regressions sharing the same design matrix. The error terms of the anchor points ϵi, i = 1, . . . ,m are defined as

OLS is indeed the best linear unbiased estimator (BLUE) by the Gauss-Markov theorem, given that the errors have expectation zero conditional on the independent variables (E(ϵi|c) = 0), are uncorrelated (Cov(ϵi, ϵj) = 0, ∀i ̸= j) and have constant variances (ϵi ∼ N(0|σ2 )) [Groß, 2003]. However, if these assumptions are violated, the robustness of OLS can be doubtful.

The situation where the constant variance assumption (i.e., homoscedastic) is violated is called heteroscedastic. In many real world datasets, we observe that heteroscedasticity indeed exists and hurts the quality of the Nystr¨om approximation if random sampling scheme is applied. Therefore, we mainly deal with the heteroscedasticity problem in this paper.

3.1 Heteroscedasticity and Its Detection

Heteroscedasticity refers to the situation where the error terms have non-constant variances. The presence of heteroscedasticity invalidates many statistical tests of significance. Thus it is a major concern in the application of regression analysis. In the linear regression model, a test for heteroscedasticity is a test of the null hypothesis [Groß, 2003]

Here the error ϵi is the i-th element of the vector of least squares residuals fy -Cˆβ (also cf. Equation (11)). The alternative hypothesis is rather unspecified. Several statistical tests can be applied to detect the existence of heteroscedasticity, including theWhite test and the Breusch-Pagan test. In this paper, we choose the Breusch-Pagan test as the tool, in which the alternative H1 is expressed as [Breusch and Pagan, 1979]:

where zTi , i = 1, . . . ,m, is the explanatory variables and α ∈ Rq is unknown parameters vector, q is the dimension of the explanatory variable. The function h(·) needs not to be specified. The Breusch-Pagan test performs a series of auxiliary regressions to get the test statistics Q, which should be asymptotically χ2(p) distributed under H0. The hypothesis H0 is rejected at level α is Q > χ2(p),1-α. We use the p-value of the χ2 distribution as the indicated value of homoscedasticity; the higher p-value we get from the Breusch-Pagan test, the higher probability Q is asymptotically χ2(p) distributed under H0, the more homoscedastic the regression model is. We will show that, starting with synthetic datasets, severe heteroscedasticity problem exists on many kinds of data distributions if the anchor points C are chosen randomly, and this problem will hurt the linear regression, and consequently the Nystr ¨om method.

We introduce two synthetic datasets to show the existence of heteroscedasticity in the Nystr¨om method with randomly sampled anchor points. We generated 1000 training instances, each instance consists of 100 observations sampled from a univariate Gaussian distribution or a lognormal distribution. The Radial Basis Function k(xi, xj) = exp(-‖xi -xj‖2 /γ) is used as the kernel function, with γ set to be the average square distance between training instances, following [Zhang et al., 2008].

The Breusch-Pagan test is used to test the heteroscedasticity of the Nystr¨om approximation on these two synthetic datasets. The homoscedasticity is measured by p-value, and the approximation error (shown in italic inside the parentheses) is measured by

where ‖·‖F is the Frobenius norm. For both random sampling and K-means sampling, we sampled 100 training instances as anchor points. We repeat the approximation experiment and the test 10 times and use the average values as the final results. Due to space limit, we omit the standard deviations of approximation errors since they are small.

As shown in Table 1, the Nystr¨om method with random sampling has poor homoscedasticity (below 0.5) while the kmeans Nystr¨om has very good homoscedasticity (about 0.9), leading to better approximation results. We also observe the same phenomenon in most real world datasets, which we will present in details in Section 4.

3.2 Heteroscedasticity, Skewness, and Cure

A variety of reasons could account for the heteroscedasticity, including aggregated and grouped data, the presence of outliers [Groß, 2003], the regression model is not correctly specified, and the skewness in the distribution of one or more regressors (explanatory variables) in the model [Fox, 1997].

One common way to overcome heteroscedasticity is to use weighted ordinary least square (WOLS) estimation for linear regression models. The WOLS model is [Groß, 2003]:

Here W is a diagonal weight matrix estimated by the error term ϵi. Unfortunately, if we apply WOLS to solve the heteroscedasticity problem in the Nystr¨om case, we have to estimate a different W for every xi ∈ D (cf. Equation (5)), which will lead to O(mn2 ) more computations. Therefore, it is not feasible to apply WOLS here unless a universally valid W for the whole dataset exists and can be efficiently found.

Alternatively, we conjecture that in the Nystr¨om method, positively skewed explanatory variables will cause severe heteroscedasticity. Skewness is a measure of the asymmetry of the probability distribution of a real-valued random variable.

the mean µ, σ is the standard deviation, S is the explanatory variables (the first term of the right-hand side of Equation (10)), and E is the expectation operator.

As shown in Table 2, positively skewed explanatory variables distributions tend to coexist with heteroscedasticity. Empirical results of many real world datasets also prove that, illustrated in Figure 1. When the skewness of explanatory variables is higher than certain threshold, the homoscedasticity drops remarkably. Empirically, we find out that explanatory variable distributions with skewness higher than 1.5 will have the heteroscedastic problem, and homoscedasticy assumption is satisfied if the p-value is above 0.7 in the Nystr ¨om approximation.

For random sampling, the position of the anchor points is unpredictable and could be in the outskirts of the datasets, which may cause positively skewed explanatory variable distributions on many datasets since kernels are mostly distance based (e.g. RBF kernel). On the other hand, if the anchor points are chosen as the centroids of K-means clusters, they are more evenly distributed, thus the explanatory variable distributions would be less biased. The fact that positions of the anchor points affects the distributions of explanatory variables and positively skewed distributions cause heteroscedasticity explain why K-means sampling has better homoscedasticity.

With the linear regression interpretation, we are able to propose to use non-symmetric explanatory function to solve the heteroscedasticity problem without changing the sampling scheme. In the statistical literature, data transformation is widely used to cure the positively skewed data [Tukey, 1977], especially a family of rank-preserving transformation named the power transformation. We utilize this well-known fact, and design an innovative explanatory function:

Here T(·) is the transformation function and 1 is a bias term for further stabilizing the variance of explanatory variables. Correspondingly, the design matrix is then:

Then the kernel matrix is approximated by:

where ̂E is the matrix of the new explanatory functions ˆ e(xi), i = 1, . . . , n. Notice that all the linear regression functions still share the same design matrix, thus the efficiency of the Nystr¨om method is preserved. Since our explanatory function is not the same as the similarity function, the approximation matrix ˜KT is non-symmetric. Therefore, we use a simple method to symmetries it:

We tested two transformation functions: T1(x) = ln(1+x) and T2(x) = x1/2 . They both are effective on curing the positively skewed data and make it more normal distribution-like. As a result, the heteroscedasticity problem is greatly relieved and the approximation quality is improved. The effects of our method can be demonstrated in the two synthetic datasets as show in Tables 1 and 2. Clearly, both T1 and T2 have good effects on correcting the positively skewed data, and since T2 are more powerful in the correction, it has better approximation quality (almost as good as k-means) than T2. Details of our method can be found in Algorithm 1. We will show the experimental results on real world datasets in Section 4.

4 Experimental Results

In this section we presents empirical evaluations of our transformation methods. We mainly examine the performance of the low-rank approximation methods by measuring their approximation errors. The methods we tested are the Nystr¨om method, K-means sampling Nystr¨om, and T1 and T2 transformed Nystr¨om. The computational complexities of different method are summarized in Table 3, where d is average non-zero dimensions of the dataset. Since high dimensional data naturally arises from many applications, such as computer vision, d is likely to be much larger than m, leading to extremely expensive K-means cost. For our method, the only extra computation cost is the transformation function T, which is just O(nm), much cheaper than K-means sampling and WOLS. We also shown the real CPU time for the Algorithm 1 Reduced Heteroscedasticity Linear Regression Nystr ¨om Approximation

• 1: Given a dataset D = {xi}ni=1, randomly sample m anchor points as D = {ci}mi=1
• 2: Calculate explanatory matrix E (which is the same as S) and X.
• 3: Calculate the skewness of E, denoted as SKE.
• 4: if SKE > thres then
• 5: Apply transformation function T1 or T2 to E and X to get ̂E and ̂X.
• 6: The approximation is then ˜KT = ̂ÊX+ST .
• 7: Symmetrize ˜KT with Equation (19).
• 8: else
• 9: Use the original Nystr¨om approximation ˜K= EX+ST .
• 10: end if
• 11: Output: The approximated matrix ˜KT or ˜K .

scene15 dataset, a moderate size vision dataset, with 100 anchor points. Our method has almost the same computational time as random sampling, while K-means is about 10 times slower.

We choose a number of benchmark datasets from the LIBSVM archive. 1 The vision dataset scene15 is from [Lazebnik et al., 2006]. We generate the dataset using Bag-ofWords model and Spatial Pyramid with libHIK [Wu, 2010]. The datasets summary is shown in Table 5. Gaussian kernel k(xi, xj) = exp(-‖xi - xj‖2 /γ) is used here with γ set to be the average square distance between training instances. To produce statistically reliable results, we repeat all the methods 10 times to get the mean and standard deviations.

1 http://www.csie.ntu.edu.tw/∼cjlin/libsvmtools/datasets/

The approximation errors are shown in Figure 2. The xaxis shows the ratio of m/n, i.e. the sample rate, which we choose from 0.01 to 0.1, following [Zhang et al., 2008]. The y-axis shows the approximation errors as defined in Equation (14). The experiments demonstrate that K-means sampling still provides the best approximation results. However, as we have discussed before, K-means clustering itself is very expensive and can be much slower than the Nystr¨om method. For large scale problem, it is not even feasible to performK-means clustering [Kumar et al., 2012]. Our methods are shown as LogNytr¨om and SqrtNytr¨om, corresponding to the T1 and T2 transformation functions. These two transformations both show good improvements on the random sampling with almost no extra cost, especially the SqrtNytr¨om. For several datasets, the SqrtNytr¨om exhibits comparable results to K-means sampling. We have also tested a simple implementation of WOLS. The weight matrix is evaluated from error terms of each linear regression. Preliminary results show that WOLS is more costly than K-means and inferior in accuracy, making it unsuitable to be used in practical situations.

We again examine the relationship between homoscedasticy, skewness and the approximation quality. We sampled 100 anchor points from each dataset to do the BreuschPagan test and repeat 10 times for each method. The results shown in Table 4 again supports our analysis: Positively skewed explanatory variable distributions from random sampling scheme leads to heteroscedasticity, and thus affects the approximation quality. Apart from sampling good anchor points using K-means clustering to avoid the skewness problem, applying transformation functions can serve well for the purpose, too, with almost no cost.

4.1 Failure Mode Analysis

Our transformation works well for positively skewed explanatory variable distributions, which widely exists in real world datasets if random sampling scheme is applied. However, for datasets with normal-like explanatory variable distributions or if anchor points are K-means sampled, our method may not have positive effects, because they do not suffer from the heteroscedasticy problem.

We have tested over 16 datasets, among which 6 are without the heteroscedasticity problem. An example of such datasets is svmguide1, whose approximation results are illustrated in Figure 3. We do not show others due to space limit. As shown in Table 4, the distribution of explanatory variables is not positively skewed and there is no heteroscedasticy issue. Therefore, the transformations do not improve the approximation accuracy. K-means sampling, on the other hand, still have positive effect, but the improvement is not as significant as on the heteroscedastic datasets. This is also an evidence for our reasoning for the success of K-means sampling.

Discussions

Learning from the failure cases, we suggest to calculate the skewness of explanatory variable distribution before applying data transformations. Since calculating skewness only involve the kernel values, which are just scalars, it is also very cheap compared to K-means clustering. As we have discussed in previous section, the empirical threshold we suggest is 1.5. By setting up the threshold, we can successfully avoid the failure cases and make good use of the transformation to get better approximation results, thus Algorithm 1 can be safely applied on almost all datasets.

5 Conclusions

In this paper, we present a linear regression interpretation for the Nystr¨om approximation and propose to use nonsymmetric explanatory functions to improve the quality of the Nystr¨om approximation with almost no extra computational cost. To be specific, we show that theoretically the the Nystr¨om method can be regraded as a set of point-wise ordinary least square linear regressions of the kernel matrix, sharing the same design matrix. We analyze the approximation quality based on the fulfillment of the homoscedasticity assumption and explain the success and deficiency of random and K-means sampling. Empirically, we demonstrate that positively skewed explanatory variable distributions can lead to heteroscedasticity, thus we propose to use the transformation functions to normalize the skewed distribution and overcome heteroscedasticity. In large scale problems, when K-means clustering is not feasible, our method can still be applied efficiently. The experimental results on various real world datasets prove the effectiveness of our method.

In the future, we want to explore the possibility to apply transformations to the ensemble Nystr¨om method [Kumar et al., 2012] to further improve our method. The ensemble Nystr ¨om method combines multiple Nystr ¨om approximations, which is called expert, with several weight functions. If we apply data transformation to each expert and also combine them with appropriate weight function, we expect to get better results on heteroscedastic datasets. We also want to found a theoretical analysis on the relationship between positively skewed regressors and heteroscedasticity so that our method is more theoretically sound.