Data Sparseness in Linear SVM

1 Introduction

Large sparse datasets are common in many real-world applications. They contain millions of data instances and attributes, but most of the attribute values are missing or unobserved. A good example is the user-item data of various popular on-line services. Netflix has published part of its movie rating dataset [Bennett and Lanning, 2007] for algorithm competition. This data consists of 100,480,507 ratings given by 480,189 users to 17,770 movies, which amounts to a sparseness of 98.822%. Data sparseness becomes even higher in other domains. For example, the Flickr dataset collected by [Cha et al., 2009] contains the 'favorite' marks given by 497,470 users on 11,195,144 photos, its sparseness reaches 99.9994%. Such high sparseness is understandable: compared to the massive amount of available items, each user could only have consumed or interacted with a tiny portion of them. In many situations, we need to classify these datasets to make predictions. For example, in personalized recommender systems and targeted advertising, we would like to predict various useful features of a user (gender, income level, location and so on) based on items she/he has visited.

1 Charles X. Ling (cling@csd.uwo.ca) is the corresponding author.

Linear SVM such as Pegasos [Shalev-Shwartz et al., 2011] and LibLinear [Fan et al., 2008] are popular choices to classify large sparse datasets efficiently, because they scale linearly with the number of non-missing values. Nowadays, it is possible to train linear SVM on very large datasets. Theoretically, it is known that larger training data will lead to lower generalization error, and asymptotically it will converge to the lowest error that can be achieved [Vapnik, 1999]. However, it is still hard to answer the following important questions about linear SVM and data sparseness:

• (1). If we put in effort to reduce data sparseness, would it decrease the asymptotic generalization error of linear SVM?
• (2). Would data sparseness affect the amount of training data needed to approach the asymptotic generalization error of linear SVM?

These questions essentially concern the convergence behavior of learning, which has been addressed in previous works. In Statistical Learning Theory [Vapnik, 1999], PAC bound gives high-probability guarantee on the convergence between training and generalization error. Once the VCdimension of the problem is known, PAC bound could predict the amount of training instances needed to approach the asymptotic error. Other works [Bartlett and ShaweTaylor, 1999] have shown that the VC-dimension of linear SVM is closely related to the hyperplane margin over the data space. As an initial step of understanding the impact of data sparseness on the margin of hyperplanes, a recent work [Long and Servedio, 2013] has given bounds for integer weights of separating hyperplanes over the k-sparse Hamming Ball space x ∈ {0, 1}m≤k (at most k of the m attributes are non-zero). There also exist several SVM variants that could deal with missing data [Pelckmans et al., 2005; Chechik et al., 2008]. However, still no previous work could explicitly predict the convergence behavior of linear SVM when data is highly sparse.

In this paper, we will answer this question by systematic experiments and then verify our findings through theoretic study. We propose a novel approach to generate large and sparse synthetic data from real-world datasets. First, we infer the statistical distribution of real-world movie rating datasets using a recent Probabilistic Matrix Factorization (PMF) inference algorithm [Lobato et al., 2014]. From the inferred distribution, we then sample a large number of data instances following the PAC framework, so we can study the generalization error with various training sizes.

In order to study the effect of data sparseness, we consider a simple data missing model, which allows us to systematically vary the degree of data sparseness and compare the learning curves of linear SVM. To follow the PAC framework, we study the curves of training and testing error rates as we keep increasing the training size. We have made several important observations about how data sparseness would affect the asymptotic generalization error and the rate of convergence when using linear SVM; see Section 3.

We then analyze our observations with detailed theoretic study. For asymptotic generalization error, we study the change of Bayes risk as we vary data sparseness. With proper assumptions, we have proved that higher sparseness will increase the Bayes risk of the data; see Section 4.1. For the asymptotic rate of convergence, we observe that different sparseness would not change the amount of training data needed to approach convergence. We study its theoretic reason using the PAC bound; see Section 4.2.

Our results are very useful for using linear SVM in realworld applications. Firstly, they indicate that sparser data would generally increase the asymptotic error, which encourages practitioners to put effort in reducing data sparseness. Secondly, our results also imply that sparser data will not lead to slower learning curve convergence when using linear SVM, although the asymptotic generalization error rate would increase.

The rest of the paper is organized as follows. Section 2 gives details of our data generation approach. Section 3 describes our experiment and observations. Section 4 proves our observations. The final section concludes the paper.

2 ANovel Approach to Generate Sparse Data

To systematically study the impact of data sparseness on linear SVM, in this section, we first propose a novel approach to generate sparse data.

2.1 Basic Settings

To simplify the problem setting, we only consider binary attribute values x ∈ {0, 1}m and binary label y ∈ {+,-}. Our data generation approach follows the standard setting of the PAC framework [Vapnik, 1999], which assumes that data (x, y) are independently sampled from a fixed distribution P(x, y) = P(y)P(x|y). In addition, our approach uses distribution inferred from real-world datasets, so the generated data are realistic in a statistical sense. Specifically, we use datasets of recommendation systems such as movie ratings to infer P (x, y).

Recently, datasets of recommendation systems are widely studied. They usually contain ordinal ratings (e.g., 0 to 5) given by each user to each item. We consider each user as a data instance and its rating on items corresponds to each attribute. To make the data binary, we only consider the presence/absence of each rating. We choose the gender of each user as its label y ∈ {+,-}.

To infer the distribution P(x, y), we use Probabilistic Matrix Factorization [Mnih and Salakhutdinov, 2007] which is a widely-used probabilistic modeling framework for user-item ratings data. In the next subsection, we will briefly review the PMF model for binary datasets proposed by [Lobato et al., 2014], which will be used in our data generation process.

2.2 Review of PMF for Sparse Binary Data

Consider L users and M items, the L × M binary matrix X indicates whether each rating exists. In this case, X could be modeled by the PMF given in [Lobato et al., 2014] which has good predictive performance:

where a Bernoulli likelihood is used along with the Matrix Factorization assumption X ≈ U · V to model each binary entry of X:

where δ(·) is the logistic function δ(x) = 1 1+exp(-x) to squash a real number into the range of (0, 1), Ber(·) is the pdf of a Bernoulli distribution and z is a global bias parameter in order to handle data sparsity.

This PMF model further assumes that all latent variables are independent by using fully factorized Gaussian priors:

where N(·|m, v) denotes the pdf of a Gaussian distribution with mean m and variance v. Given a binary dataset X, the model posterior p(U, V, z|X) could be inferred using the SIBM algorithm described in [Lobato et al., 2014]. The posterior predictive distributions are used for predicting each binary value in X:

2.3 The Distribution for Data Sampling

Based on the above PMF model, we next describe the distribution for sampling synthetic data. Since the distribution is inferred from a given real-world dataset X, we denote it as

For p(y), we simply use a balanced class prior, p(+) = p(-) = 0.5. This ensures equal chances of getting positive and negative instances when sampling.

Now we describe how we get p(˜x|y,X). We first divide the real-world dataset X into the positive labeled set X+ and the set of negative instances, X-. We infer the model posteriors p(U,V, z|X-) and p(U,V, z|X+) separately using the SIBM algorithm.

Notice that these inferred models can not be directly used to generate infinite samples: suppose the number of users in X+ (and X-) is L+ (and L-), the posterior predictive (6) only gives probabilities for each of these L+ (L-) existing users. In other words, the predictive distribution could only be used to reconstruct or predict the original X, as did in the original paper [Lobato et al., 2014].

In order to build a distribution for sampling an infinite number of synthetic users, we employ the following stochastic process: whenever we need to sample a new synthetic instance ˜x, we first randomly choose a user i from the L+ (or L-) existing users, and we then sample the synthetic instance using the posterior predictive of this user. Using this process, the pdf of p(˜x|y,X) is actually

where y ∈ {+,-}, p(i, y) = 1 Ly and

is the likelihood for each instance (existing user) of Xy, which is equivalent to (1). The process of sampling data from p(˜x, y|X) can be implemented by the following pseudo code.

Algorithm 1 Data Sampling

infer p(U, V, z|X-) and p(U, V, z|X+) using SIBM; sample V+, z+ from p(U, V, z|X+); sample V-, z- from p(U, V, z|X-); for each new instance (˜x, y) do randomly sample y from {+, -}; randomly sample i from {1, ..., Ly}; sample ui from p(U,V, z|Xy); for j in 1...M do sample ˜ xj fromBer(˜ xj|δ(ui · vj,y + zy)); end for end for

This algorithm allows us to sample infinite instances from the fixed distribution p(˜x, y|X), to be used for generating data of various sizes. Next, we will discuss how to systematically vary sparseness of the sampled data.

2.4 DataMissingModel

We employ a simple data missing model to add and vary data sparseness. Our data missing model assumes that each attribute has the same probability to become 0, which follows the Missing Completely At Random assumption [Heitjan and Basu, 1996].

Given the probability of missing s, the data missing model will transform an instance ˜x = (˜ x1, ˜ x2, ..., ˜ xm) to x = (x1, x2, ..., xm) following:

To ease our illustration, we hereafter call this process as dilution, and the resultant data as the diluted data. Now if instances are originally sampled from p(˜x|y), the distribution of the diluted data can be computed by p(x|y, s) = ∫ p(x|˜x, s)p(˜x|y)d˜x. When we apply this dilution process to data generated from Algorithm 1, we denote the resultant distribution as p(x|y, X, s) = ∫ p(x|˜x, s)p(˜x|y, X)d˜x.

Using the above data missing model, additional sparseness will be introduced uniformly to each attribute and higher s will lead to sparser data. This enables us to vary s and study the impact of data sparseness systematically, as we will describe in the next section.

3 Experiment

In this section, we describe how we use the proposed data generation and dilution process to study our research question with experiments.

Data. We use two movie rating datasets: the movielens 1M dataset 1 (ml-1m) and the Yahoo Movies dataset 2 (ymovie). As mentioned earlier, we only consider absence/presence of each rating. The preprocessed ml-1m dataset is 3, 418 (instances) ×3, 647 (attributes) with balanced classes and an original sparseness of 0.9550. The preprocessed ymovie dataset is 9, 489 (instances) ×4, 368 (attributes) with balanced classes and an original sparseness of 0.9977.

Training. To study the impact of data sparseness, we use various values of s; see the legends of Figure 1 and 2. For each s, we generate training samples of various sizes l from p(x, y|X, s) using the proposed data generation and dilution process. We use Lib-linear [Fan et al., 2008] with default parameters to train the classifier and get the corresponding training error ϵtrain(l, s).

Testing. It is empirically impossible to test a classifier on the whole distribution and get the true generalization error. For this reason, we generate a large test dataset from p(x, y|X, s) for each setting of s. For the ymovie experiment, the size of test sets is 7.59 million instances (800 times of the real data size). For the ml-1m experiment, the test size is 0.68 million (200 times of the real data size). We test the classifiers on the corresponding test set and get the test error ϵtest(s).

For each setting of s and training size l, we repeat the data generation (including dilution), training and testing process for 10 times, and record the average training and testing error rates. The most time consuming step in our experiment is data generation: the largest training dataset (l = 220 in the ymovie experiment) costs us one day on a modern Xeon CPU to generate each time; Meanwhile, training on this particular dataset only costs several minutes, thanks for the efficiency of linear SVM. To speed up, we straight-forwardly parallelize the data generation jobs on a cluster.

As we systematically vary s, the resultant learning curves of averaged training and testing error rates are given in Figure 1. Figure 2 shows the difference between the two errors in order to show the rate of convergence. We have two important observations which are obvious in both experiments:

1 http://grouplens.org/datasets/movielens/

2 http://webscope.sandbox.yahoo.com/catalog.php?datatype=r

Observation 1. asymptotic generalization error: Higher sparseness leads to larger asymptotic generalization error rate;

For example, in the ymovie experiment, the asymptotic (training size l > 106 ) generalization error rates are 16.6%, 22.4%, 33.7% and 39.2% for s = 0.5, 0.7, 0.9 and 0.95, respectively.

Observation 2. asymptotic rate of convergence: the asymptotic rate of convergence is almost the same for different sparseness.

In other words, different sparseness would not change the amount of training data needed to approach convergence. For example, in the ymovie experiment, the minimum training size l for the two error rates to be within 1% is almost the same (l = 370K), regardless of the value of s.

In the next section, we will study the theoretic reasons behind these observations.

4 Theoretic Analysis

In section 4.1, we study the theoretic reason for observation 1, and in section 4.2 we show that the theoretic reason of observation 2 can be easily explained using PAC bound.

4.1 Asymptotic Generalization Error

Suppose: (1) the original data distribution is p(x|y); 3 (2) after applying the data missing model described in Section 2.4, the diluted data distribution is p(x|y, s); (3) class prior is balanced, i.e., p(+) = p(-) = 0.5. From the Bayesian Decision Theory [Duda et al., 2012], we know that the asymptotic generalization error rate of linear SVM can be lower bounded by the Bayes risk, which is the best classification performance that any classifier can achieve over a given data space:

where R(y(x)|x) is the loss for predicting the label of instance x as y, which in our case is the 0-1 loss. We denote the Bayes risk for p(x|y, s) as R(s), thus the asymptotic generalization error rate is lower bounded:

Notice x lives in the discrete data space {0, 1}m, which eases us to write the Bayes risk as the following form:

3 Our proof does not require p(x|y) to be the specific distribution p(˜x|y,X) introduced earlier.

We will next prove that higher s leads to larger R(s) using three steps. We first consider the case if only one of the attributes is diluted. In this case, we could prove that the Bayes risk will not decrease by any chance, and with high probability it will increase (Lemma 4.1). We next consider the case if we still dilute only one attribute but with different s, we prove that higher sparseness will lead to larger Bayes risk (Lemma 4.2). Based on these results, we finally prove that higher s leads to larger Bayes risk R(s) (Theorem 4.3).

When we only dilute one of the m attributes, xj , we denote the rest of the attributes as x-j = (x1, ...xj-1, xj+1, ..., xm). Since the order of attributes does not matter, we denote x as (x-j, xj). We denote the corresponding distribution as p(j) (x|y, s) and Bayes risk as R(j) (s), we now prove:

Lemma 4.1 R(j) (s) ≥ R(0) always holds. Specifically, with high probability, R(j) (s) > R(0).

Proof Since we have only diluted xj, we first expand the computation of Bayes risk (13) along xj :

where Z(x-j) denotes the sum of probability mass of p(x|y) at (x-j , 0) and (x-j , 1), each minimized among classes:

Z (j) (x-j, s) is also defined accordingly for p(j) (x|y, s). ∑x-j denotes the summation over all x-j in the {0, 1}m-1 space. Now define ∆Z(s) := Z (j) (x-j , s)-Z(x-j ), we will next prove that ∆Z(s) ≥ 0 holds for all x-j ∈ {0, 1}m-1 :

Since xj is diluted with s, according to (10), this means ∀x-j ∈ {0, 1}m-1 , we have

We denote yh ∈ {-, +} as the label that has a higher probability mass at (x-j , 1), and we denote the other label as yl: 4

and (17)(18) lead to:

In order to write out miny∈{+,-} p(x-j, 0|y) and miny∈{+,-} p(j) (x-j , 0|y, s), we define

4 Notice that yh and yl will change for different x-j. Strictly, we should use the notation yh(x-j, 1) and yl(x-j, 1) if not for the purpose of brevity.

For each x-j there are in total three different situations to consider:

We could straight-forwardly compute ∆Z(s) for each case. Case 1:

Case 2:

Case 3:

For space limit, detailed derivation of (25)(26)(27) is omitted.

Now that ∆Z(s) ≥ 0 always holds, R(j) (s) ≥ R(0) is true because of (14). Moreover, R(j) (s) = R(0) is true only if Case 1 happens for all x-j ∈ {0, 1}m-1 , which has low probability.

In the next lemma, we consider the case when we vary the sparseness s on the one diluted attribute.

Lemma 4.2 ∀s1, s2 s.t. 1 ≥ s2 > s1 ≥ 0, we have R(j) (s2) ≥ R(j) (s1). Specifically:

(1). only when both s1 and s2 are close enough to 1.0, R(j) (s2) = R(j) (s1) will be true with high probability.

Its proof could be derived by studying the value of ∆Z(s) in different situations. It is straight-forward after we have derived equations (25)(26)(27) in the proof of Lemma 4.1.

Proof We first prove that ∆Z(s2) ≥ ∆Z(s1) always holds. For any given x-j ∈ {0, 1}m-1 :

(1). If Case 1 is true, we will always get ∆Z(s2) = ∆Z(s1) = 0 regardless of the value of s1 and s2.

(2). If Case 1 is false and 1 ≥ s2 > s1 ≥ g(x-j), then Case 2 holds for both s1 and s2, and we get

(3). If Case 1 is false and s2 ≥ g(x-j) > s1, then Case 2 holds for s2 and Case 3 holds for s1. Thus

(4). Finally, if Case 1 is false and g(x-j) > s2 > s1 ≥ 0, then Case 3 holds for both s1 and s2. From (27) we know that ∆Z(s2) > ∆Z(s1).

Now we have proved that ∆Z(s2) ≥ ∆Z(s1) always hold, substitute into (14) leads to

From the above analysis, we further notice that R(j) (s2) = R(j) (s1) only happens if for ∀x-j ∈ {0, 1}m-1 either Case 1 happens or Case 2 holds for both s1 and s2. Notice for given p(x, y), it is very unlikely that all x-j could satisfy Case 1. For those x-j that Case 1 is not satisfied, it thus requires Case 2 to hold for both s1 and s2, which is only likely when both s1 and s2 are close enough to 1.0.

The next theorem will generalize to the case when sparseness is introduced to all attributes, which is the goal of our proof:

Theorem 4.3 ∀s1, s2 s.t. 1 ≥ s2 > s1 ≥ 0, we have R(s2) ≥ R(s1). Specifically:

(1). only when both s1 and s2 are close enough to 1.0, R(s2) = R(s1) will be true with high probability.

(2). other wise, R(s2) > R(s1).

The basic idea of our proof is to find a hypothetic process that varies the dilution of data from s2 to s1 one attribute at a time, so we could leverage the result of lemma 4.2.

Proof We use F to denote the total attribute set {x1, ..., xm}. At a certain state, we use F(s2) ⊂ F to denote the set of attributes that have been diluted by s2; and F(s1) ⊂ F the set of attributes that have been diluted by s1. We denote the corresponding distribution as pF (s1),F (s2) (x|y, s1, s2) and its Bayes Risk as RF (s1),F (s2) (s1, s2).

Now we consider the following process that iteratively changes the elements of F(s1) and F(s2): we start from F(s1) = ∅ ∧ F(s2) = F, and move each attribute in F(s2) to F(s1) one after another. After m such steps we come to F (s1) = F ∧ F (s2) = ∅.

After step i, we denote the current F(s1) as Fi(s1), F(s2) as Fi(s2). Using this notation, the initial state is F0(s1) = F ∧F0(s2) = ∅, and final state is Fm(s1) = ∅∧Fm(s1) = F .

We next prove that for ∀i ∈ {0, ...,m- 1}:

Notice that in each step i, there is only one attribute (e.g., xj) changes fromFi(s2) to Fi(s1), i.e., F-[Fi(s1)∪Fi+1(s2)] = {xj }.

Now consider the distribution pFi(s1),Fi+1(s2) (x|y, s1, s2), which corresponds to an intermediate state between step i and i + 1 with xj being not diluted. Now if we dilute xj by s2, we get pFi+1(s1),Fi+1(s2) (x|y, s1, s2); Instead, if we dilute xj by s1, we get pFi(s1),Fi(s2) (x|y, s1, s2). Using Lemma 4.2, it now becomes obvious that (29) is true, which further leads to

We further notice that R(s2) = R(s1) only holds when all m equal signs hold simultaneously. As illustrated in Lemma 4.2, each equal sign holds with high probability only when both s2 and s1 are close enough to 1.0.

This theorem tells us that higher sparseness leads to larger Bayes risk, which explains our observation on the asymptotic generalization error.

This result is important for understanding how data sparseness affects the asymptotic generalization error of linear SVM.

Since the Bayes risk of the data space is irrelevant to the classifier in use, our theoretic result is also applicable for other classifiers.

4.2 Asymptotic Rate of Convergence

Our observation on the asymptotic rate of convergence could be explained by the PAC bound [Bartlett and Shawe-Taylor, 1999]. From the PAC bound, we know that the convergence of generalization error rate ϵ and the training error rate ϵtr is bounded in probability by training size l and the VCdimension |d|:

For fixed δ, the rate of convergence is approximately:

Though different data sparseness could change |d|, for linear SVM, |d| < m+ 1 holds true regardless of data sparseness. 5 From (32) we could see that asymptotically (when 2l3 ≫ ln m+1 δ ), the rate of convergence is mostly influenced by the increase of training size l rather than by the change of |d|, since |d| is always upper bounded by m+ 1. In other words, varying sparseness will have little impact on the asymptotic rate of convergence, which verifies our observation.

When using linear SVM in real-world applications, it is important to know whether sparser data would lead to slower convergence rate. If so, practitioners will have to collect more training instances in order for linear SVM to converge on highly sparse datasets. Our observation and analysis shows that the rate of convergence is almost not affected by different data sparseness for linear SVM.

5 Conclusion

Linear SVM is efficient for classifying large sparse datasets. In order to understand how data sparseness affects the convergence behavior of linear SVM, we propose a novel approach to generate large and sparse data from real-world datasets using statistical inference and the data sampling process of the PAC framework. From our systematic experiments, we have observed: 1. Higher sparseness will lead to larger asymptotic generalization error rate; 2. Convergence rate of learning is almost unchanged for different sparseness. We have also proved these findings theoretically. Our experiment and theoretic results are valuable for learning large sparse datasets with linear SVM.

5 Notice that VC-dimension can be very large for some non-linear versions of SVM, how their convergence rate will be affected by data sparseness could be an interesting future work.

Acknowledgments

This work is supported by National Natural Science Foundation of China (grants 61432020, 61472430, 61202137) and Natural Sciences and Engineering Research Council of Canada (NSERC).