Approximate Manifold Regularization: Scalable Algorithm and Generalization Analysis

1 Introduction

Recently, the explosive growth of computing power and applications of the network makes data generation and acquisition more easily. However, most of the collected data are unlabeled, while data annotation is laborious. Further, semisupervised learning (SSL) methods are developed to estimate specific learner from a few labeled samples together with a significant amount of unlabeled data, such as transductive support vector machines [Joachims, 1999] and graphbased methods [Belkin et al., 2006; Camps-Valls et al., 2007]. Graph-based manifold regularization methods draw wide attention of SSL area due to their good performance and relative simplicity of implementation [Belkin et al., 2006]. Despite those advantages of manifold regularization, it remains challenges to process gigantic datasets, for suffering high computational complexity, typically kernel matrix related operations at least O(n2 ) and construction of graph Laplacian at least O(n log n), where n is total sample size.

To tackle those scalability issues, many approaches were proposed [Liu et al., 2012; Jiang et al., 2017; Liu et al., 2019]: (1) Accelerate construction of Laplacian graph. Methods based on the fast spectral decomposition of Laplacian matrix have been well-studied in [Talwalkar et al., 2013], which use a few eigenvalues of graph Laplacian to represent manifold structure. Graph sparsification approaches were devised to approximate Laplacian graph by a line or spanning tree [Cesa-Bianchi et al., 2013] and improved by minimizing tree cut (MTC) in [Zhang et al., 2016]. (2) Accelerate operations associated with kernel matrix. Several distributed approaches have been applied to semi-supervised learning [Chang et al., 2017], decomposing a large scale problem into smaller ones. Anchor Graph regularization (Anchor) constructs an anchor graph with the training samples and a few anchor points to approximate Laplacian graph [Liu et al., 2010]. The work of [McWilliams et al., 2013; Rastogi and Sampath, 2017] applied random projections including Nystr¨om methods and random features into manifold regularization. Gradient methods are introduced to solve manifold regularization on the primal problem, such as preconditioned conjugate gradient [Melacci and Belkin, 2011], stochastic gradient descent [Wang et al., 2012].

∗Corresponding author

In this paper, we focus on the latter scalability issue. With sound theoretical guarantees, we devise a novel graph-based SSL framework, substantially reducing computational time and memory requirements. More precisely, our approach approximates Laplacian regularized least squares (LapRLS) by Nystr¨om methods and then accelerates the solution with preconditioned conjugate gradient methods. It's a non-trivial extension of FALKON [Rudi et al., 2017] to graph SSL with technical challenges in algorithm design and theoretical analysis. Theoretical analysis demonstrates that O(√m) labeled samples and O(logm) iterations (m is the number of labeled samples) can guarantee good statistical properties. Complexity analysis shows our method solve LapRLS with O(n√n) time and O(n) space (n is the number of all samples).

2 RelatedWork

To overcome the computational and memory bottleneck of LapRLS, practical algorithms were developed, including Nystr¨om methods [Williams and Seeger, 2001] of which statistical properties are well studied in [Rastogi and Sampath, 2017], and preconditioned conjugate gradient (PCG) which reduces the number of iterations [Cutajar et al., 2016]. FALKON approach combines Nystr¨om methods and PCG in supervised learning [Rudi et al., 2017]. Further, our work extends the combination to SSL with high computation gains and sound statistical guarantees. The approach improves computational efficiency from O(n3 ) to O(n√n) and reduce memory cost fromO(n2 ) to O(n).

3 Preliminaries

3.1 ProblemDefinition

Assume there is a fixed but unknown distribution ρ on X ×Y, where X = Rd and Y = R. Further, m labeled samples {(x1, y1), · · · , (xm, ym)} ∈ X×Y are drawn i.i.d fromρ and n -m unlabeled samples {xm+1, · · · , xn} ∈ X are drawn i.i.d according to the marginal distribution ρX of ρ.

3.2 Manifold Regularization

Manifold learning methods based on the spectral graph, known as graph-based SSL, is a typical solution to semisupervised learning [Zhu et al., 2003; Belkin et al., 2006], which is to find a smooth low-dimensional manifold embedded in the high-dimensional vector space, based on sample points. Correctly, Laplacian regularization [Belkin et al., 2006] is extensively used in graph-based SSL.

For a Mercer kernel K : X×X → R, there is an associated reproducing kernel Hilbert space (RKHS) H of functions f : X → R with corresponding norm ‖ · ‖H. The following optimization is considered in manifold regularization:

where ℓ is loss function, L is Laplacian matrix by L = D -W, f = [f(x1), · · · , f(xn)]T , λA controls the complexity of the function in the ambient space, and λI controls the complexity of the function in the intrinsic space. Here, W ∈ Rn×n records undirected weight between points and the diagonal matrix D is given by Dii = ∑nj=1 Wij .

The minimizer of the optimization problem (1) admits an expansion in terms of both labeled and unlabeled data

3.3 Laplacian Regularized Least Squares (LapRLS)

With squared loss function, the problem (1) becomes LapRLS

Setting the derivative of the objective function be zero leads a closed form solution

where Kij = K(xi, xj ) is n × n kernel matrix on train data, J = diag(1, · · · , 1, 0, · · · , 0) with the firstmdiagonal entries as 1 and the rest 0, and yn = [y1, y2, · · · , ym, 0, · · · , 0]T with corresponding m labels and the rest filled by 0. Note that when λI = 0, Equation (3) gives zero coefficients over unlabeled data, thus the form reduces to the standard RLS.

4 Algorithm

We devise a fast and scalable graph-based semi-supervised learning framework Nystr¨om-PCG shown as Algorithm 1, which consists of two steps: (1) Nystr¨om with uniform sampling on train data for the LapRLS problem, resulting in a linear system Hα = z. (2) Define a preconditioner P to approximate H, and then solve P-1Hα = P-1z by PCG.

4.1 Nystr¨om Subsampling on LapRLS

We consider Nystr¨om subsampling to reduce memory requirement, which uses a smaller matrix obtained from random column sampling to approximate the empirical kernel matrix. Thus, a smaller hypothesis space Hs is introduced

where s ≤ n and xs = (˜x1, · · · , ˜xs) are Nystr¨om centers selected by uniform subsampling from the training set. The minimizer of LapRLS (2) over the space Hs is in the form:

where H† denotes the Moore-Penorse pseudoinverse of a matrix H, (Kms)ij = K(xi, ˜xj) with i ∈ {1, · · · ,m} and j ∈ {1, · · · , s}, (Kss)kj = K(˜xk, ˜xj) with k, j ∈ {1, · · · , s} and y = [y1, · · · , ym]T ∈ Rm.

4.2 Solving the Linear System by Preconditioning

Nystr¨om subsampling for LapRLS problems resulting solution (4) is also a linear system, so we consider how to accelerate the solution by preconditioning that is

As we all know, the number of iterations for preconditioning methods depends on the condition number cond(P-1H), such that the preconditioner needs to be approximate to H. To obtain a smaller condition number but also avoid inefficient computation, we define the following preconditioners:

• m> √n

In each iteration of any PCG solver, calculation of Hα is needed. To accelerate computation, Hα is decomposed into a series of matrix-vector multiplications

Remark 1. We use LU or QR decomposition to calculate matrix inversion P-1 because they show significant improvement in speed than Cholesky decomposition.

Remark 2. The storage of kernel matrix Kns needs at least O(ns) memory, but it turns to be O(s2 ) when we perform matrix multiplications in s × s blocks.

Algorithm 1 Nystr ¨om LapRLS with PCG (Nystr¨om-PCG)

Input: m labeled samples {(xi, yi)}mi=1 , n - m unlabeled samples {xj}nj=m+1 . Parameters: λA, λI , kernel method K and subsampling size s.

Output: coefficients α

• 1: Construct Laplacian graph matrix L.
• 2: Select s Nystr¨om centers with uniform sampling from training set {˜x1, · · · , ˜xs} ∈ {x1, · · · , xn}.
• 3: Calculate inverse of the preconditioner P-1 (5) or (6).
• 4: Use any PCG solver to solve P-1Hα = P-1z with calculating Hα by Eq. (7) and performing matrix multiplication in blocks.

4.3 Complexity Analysis

Time Complexity

Before the start of iterations, we need to compute the preconditioner P, which needs O(ms2 + s2 + s3 ) time for (5) and O(s3 ) time for (6). The inverse of preconditioner P-1 needs O(s3 ) as well. In each iteration of PCG, calculation of Hα and P-1v are needed, where v represents a vector Hα or z. For P-1v, it just needs O(s2 ) for matrix-vector multiplication. And computation for Hα as in Eq. (7), it needs O(ms + s2 + ns + nk) where L is usually a sparse matrix which is almost at n × k size. Thus, each iteration needs O(ns) time. Combing calculation of P-1 and t iterations of PCG, total time complexity is O(s3 +nst).

Space Complexity

Main memory requirement comes from storing Kns which is at O(ns) and L which is at O(nk), where k ≪ s is a small constant. Indeed, matrix multiplications are performed in a series of s × s blocks so that space complexity can be O(s2 ).

5 Theoretical Result

Define ̂ f s λ,t as the estimator of Nystr¨om-PCG. The goal of nonparametric regression is to minimize the excess risk:

Under some basic conditions in nonparametric learning, we derive the excess risk bound for Nystr¨om-PCG.

Assumption 1 (Basic assumptions). (1) There exists κ ≥ 1 such that K(x, x) ≤ κ2 for any x ∈ X. (2) There exists fH ∈ H, such that E(fH) = inff∈H E(f).

Theorem 1. Let η ∈ (0, 1] and m0 ∈ N. Under Assumption 1 and assume |y| ≤ b, ∀b > 0, if m ≥ m0 and

then the following holds with probability at least 1 - η,

where the constant m0 is independent on λ, s,m, t and the constant c0 is independent on λ, s,m, t, η.

The above result provides the desired bound and proof details are deferred in the appendix. Under basic assumptions, the obtained learning rate is the same as the full KRR estimator and is optimal in a minmax sense [Caponnetto and De Vito, 2007]. Theorem 1 shows that O(√m) Nystr¨om labeled points and O(logm) iterations can obtain O(1/√m) learning rate. In practical, we sampled uniformly O(√m) labeled points and O(√n-m) unlabeled points, thus the size of sampled examples is O(√n). With omitting logarithmic terms, the approach needs O(n1.5 ) time and O(n) space.

6 ComparedMethods

In this part, we introduce compared methods including the standard RLS, LapRLS and Nystr ¨om LapRLS. Results in Table 1 show that the proposed approach Nystr ¨om LapRLS with PCG (Nystr¨om-PCG) can remarkably reduce memory requirements and improve computational efficiency.

Standard RLS (RLS)

Consider the standard RLS only using labeled data

where Kmm is the empirical kernel matrix on labeled data and y is corresponding labels. The computation cost is standard, requiring O(m3 ) time and O(m2 ) memory.

Laplacian Regularized Least Squares (LapRLS)

The linear system of LapRLS has been derived in (3) is

In terms of LapRLS, the space complexity is O(n2 ), but time complexity depends on optimization algorithms, O(n3 ) for matrix inversion and O(n2 t) for gradient methods. Specifically, the number of iterations is t = O(√n logn) for CG and t = O(log n) for PCG [Rudi et al., 2017].

LapRLS with Nystr¨om Subsampling (Nystr¨om)

The linear system of LapRLS with Nystr¨om subsampling to overcome memory bottleneck has been given in (4). According to Theorem 1, we set the number of Nystr¨om centers be s = O(√n). Space complexity of Nystr¨om LapRLS is O(s2 ) by matrix blocks multiplications, that is O(n). For direct methods, time complexity depends on computing H with O(ns2 ), that is O(n2 ). For CG methods, time complexity is O(nst) where t = O(√s log s), that is O(n1.75 ).

7 Empirical Study

In this section, we conduct empirical experiments of some SSL methods on a range of datasets. We compare five approaches, including standard RLS using CG, LapRLS and Nystr ¨om LapRLS in terms of CG and PCG. Nystr ¨om LapRLS using PCG is our proposed method, named Nystr¨om-PCG.

Scalable datasets from thousand to hundreds of thousand are used. We apply 8-NN to construct adjacency matrix with weight Wij = exp-‖xi-xj‖2 /4 as using in many literatures. For each dataset, we use Gaussian kernel K(xi, xj) = exp(-‖xi - xj‖/2σ2 ). We choose kernel parameter σ and regular parameters (λ in standard RLS and λA, λI in LapRLS methods) in 2i , i ∈ {-15,-14, · · · , 14, 15}, by minimizing test error via 10-folds cross-validation.

7.1 Accuracy and Speed

Using the chosen parameters determined by 10-folds crossvalidation, we run all methods 30 times with randomly select 70% for training and 30% for testing on each dataset. Meanwhile, we randomly select 10% samples (m = 0.1n) as labeled data and 10% samples (s = 0.1n) as Nystr¨om centers. The use of multiple training/testing partitions allows an estimate of statistical significance between best one and the remainder referring to 95% level of significance under t-test.

Table 2 reports the average root mean square error (RMSE) and Table 3 reports the number of iterations and running time, can be summarized as follows: (1) RLS is defended by other methods using Laplacian regularization on all datasets, while LapRLS-CG or LapRLS-PCG gives the best results almost on all datasets. (2) There is no significant difference in error rate between LapRLS methods and Nystr¨om LapRLS methods. (3) LapRLS-CG and LapRLS-PCG failed on large datasets because of memory limitation. (4) CG and PCG always result in the same accuracy, but PCG need much smaller iterations thus cost less time. (5) Our method Nystr¨om-PCG achieves similar accuracy as LapRLS methods with tens to hundreds of times speeding up.

7.2 Influence of Label Proportion

To explore influence of different label proportions, we let m vary from n × {1%, 2%, 4%, 8%, 16%, 32%, 64%} and fix sample size as s = 0.1n. For better presentation, we use LapRLS represents primal LapRLS approaches, Nystr ¨om represents approximate LapRLS using Nystr¨om. By repeating the procedure 20 times for different labeled/unlabeled splits on eight datasets, we report results in Figure 1. We can see that: (1) LapRLS methods are always better than the RLS method, but the difference becomes smaller when labeled partition m becomes larger. (2) Average accuracies of LapRLS and Nystr¨om LapRLS are close when the partition of labeled data is big enough. (3) The standard deviation of all methods becomes smaller as the increase of labeled samples.

7.3 Influence of Sample Proportion

To explore influence of different sample proportions, we let s vary from n × {1%, 2%, 4%, 8%, 16%, 32%, 64%} and fixed labeled sample as m = 0.1n. After repeating experiments on three different optimizations 20 times on different sample size but same labeled/unlabeled splits on eight datasets, we report results in Figure 2. We can see that: (1) RLS and LapRLS methods always give the same results because they are run on the same labeled/unlabeled splits. (2) Average accuracies Nystr¨om LapRLS become closer to LapRLS as increase with sample size. (3) Nystr¨om methods can achieve good approximation when the sampled proportion is larger than 10%.

8 Conclusion

In this paper, with sound theoretical guarantee, we use Nystr¨om subsampling on LapRLS and solve the linear systems with PCG, substantially reducing memory and computational costs. More precisely, the theoretical analysis provides good convergence rates for Nystr¨om subsampling and Nystr¨om subsampling with PCG, suggesting that O(√m) sample size and O(logm) iterations can achieve O(1/√m) error bound. Then complexity analysis show our method achieve good statistical accuracy with O(n√n) time and O(n) space. Empirical results show our method achieves similar prediction accuracy to LapRLS with higher computational efficiency and less memory requirement.

9 Proof

Based on integral operator framework and common assumptions, we firstly derive the excess risk bound for LapRLS with Nystr ¨om E ( ˜ f s λ ) - E(fH). And then we prove the excess risk bound E( ̂ f s λ,t ) - E(fH) for Nystr¨om-PCG in Theorem 1.

Firstly, we introduce definition of the integral operator and standard assumptions in statistical learning [Caponnetto and De Vito, 2007; Liu et al., 2014; Li et al., 2018] and approximation theory [Rudi et al., 2015; Rastogi and Sampath, 2017; Liu et al., 2018].

Definition 1 (Integral operator and the effective dimension). Let LK : L2 (X , ρX) → L2 (X , ρX) be integral operator

where L2 (X, ρX) = {f : X → R | ‖f‖2ρ = ∫ |f(x)|2dρX < ∞} and the effective dimension is defined as

The effective dimension N(λ) measures the size of H and there exists Q > 0 such that always holds N(λ) ≤ Q2λ-1 . The following Assumption 2 is satisfied when y is bounded, sub-gaussian or sub-exponential.

Assumption 2 (Moment assumption). Assume there exists M > 0 and σ > 0, such that for all l ≥ 2 with l ∈ N,

Assumption 3 (Regularity assumption). Assume there exists r ∈ [1/2, 1] and g ∈ L2 (X , ρX), such that

Assumption 3 controls the regularity of fH and is common in statistical learning [Caponnetto and De Vito, 2007]. It is always satisfied with r = 1/2.

We consider a more general multi-penalty regularization scheme for Nystr¨om subsampling LapRLS:

where Bj : H → H (1 ≤ j ≤ p) are bounded operators, λ > 0, λj (1 ≤ j ≤ p) are non-negative real numbers.

Theorem 2 (Excess risk bound for Nystr¨om LapRLS). Let η ∈ [0, 1]. Assume there exists some constants M, σ. Under Assumption 1, 2 and 3, for sufficiently large sample according to λ ≥ 8κ2 √m log ( 4 η ) and for subsampling according to s ≥ max {67 log ( 12κ2 λη ) , 2κ2 log ( 12κ2 λη )}, the convergence rate of ˜ f s λ holds with probability at least 1 - η

where the constant c1 and c2 do not depend on λ, s,m, t, η, Bλ = || ∑pj=1 λjB∗ j Bj ||, s is the sample size on labeled data.

Proof. Following proof details of Theorem 3.1 in [Rastogi and Sampath, 2017], there holds

where c1 = 6R+(5+cψ)(3+cφ)R, c2 = (5+cψ)||fH||ρ and Bλ = || ∑pj=1 λjB∗ j Bj||. Assumption 3 is always satisfied with fH(x) = (Lr g)(x) thus φ(λ) = λr-1/2 , ψ(λ) = λ1/2 and ψ(LK) = L1/2 K . There holds

Using the fact ε(f)-ε(fH) = ||L1/2 K (f-fH)||2H for any f ∈ H [Caponnetto and De Vito, 2007], we obtain the result.

Proof of Theorem 1. The proof is mainly following proof techniques of Theorem 8 and Theorem 3 in [Rudi et al., 2017]. Let ˜ f s λ,t be the estimator in Nystr¨om LapRLS optimization after t ∈ N iterations. Theorem 8 of [Rudi et al., 2017] gives the connection between E( ˜ f s λ,t ) - E(fH) and E( ˜ f s λ )-E(fH). Following proof steps in Theorem 8 of [Rudi et al., 2017], there holds

when the Nystr¨om centers are uniformly sampled and

Assumption 3 is satisfied with r = 1/2. From Theorem 3.2 of [Rastogi and Sampath, 2017], we set λj = λr+1 , thus Bλ ≤ cpλr+1 , where cp > 0 only depends on the number of penalty terms. Applying N(λ) ≤ Q2λ-1 , r = 1/2 and λ = cλm-1/2 where cλ = 8κ2 log ( 4 η ), we have

Therefore, the excess risk bound of Nystr¨om-PCG holds

where c0 = (c3c1/2 λ +c2cp+8κM+8σQ)2 . With conditions

satisfying conditions required by Theorem 8 of [Rudi et al., 2017] and Assumption 2 is satisfied by M = σ = 2b, thus

Combining (9), (10) and (11), we complete the proof.

Acknowledgments

This work is supported in part by the National Natural Science Foundation of China (No.61703396, No.61673293), the Youth Innovation Promotion Association CAS, the National Key Research and Development Program of China (No.2018YFC0823104, No.2016YFB1000604), the Science and Technology Project of Beijing (No.Z181100002718004) and the Excellent Talent Introduction of Institute of Information Engineering of CAS (Y7Z0111107).