An N logN Parallel Fast Direct Solver for Kernel Matrices Chenhan D. Yu∗, William B. March† and George Biros§ ∗Department of Computer Science ∗†‡ Institute for Computational Engineering and Science The University of Texas at Austin, Austin, Texas, USA ∗chenhan@cs.utexas.edu, †march@ices.utexas.edu, †gbiros@acm.org Abstract—Kernel matrices appear in machine learning and and may require 1000s of iterations. This complexity bar- non-parametricstatistics.GivenN pointsinddimensionsand rier has limited the use of kernel methods for large-scale 7 a kernel function that requires O(d) work to evaluate, we 1 problems [6], [20]. present an O(dNlogN)-work algorithm for the approximate 0 Contributions. We exploit hierarchically low-rank ap- factorization of a regularized kernel matrix, a common com- 2 putational bottleneck in the training phase of a learning task. proximations in which we assume that K can be approxi- n With this factorization, solving a linear system with a kernel matedwellbyD+UV,whereDisblock-diagonalandtheU a matrixcanbedonewithO(NlogN)work.Ouralgorithmonly andV matriceshavelowrank.Usingsuchadecomposition, J requireskernelevaluationsanddoesnotrequirethatthekernel weimprovethefactorizationalgorithmpresentedin[36].In matrix admits an efficient global low rank approximation. 9 thatpapertheblock-diagonalplussparsedecompositionwas Instead our factorization only assumes low-rank properties for the off-diagonal blocks under an appropriate row and done using ASKIT1, a method introduced in [21], [22] (see ] C columnordering.Wealsopresentahybridmethodthat,when §II-A). ASKIT approximates K in O(dNlogN) time and the factorization is prohibitively expensive, combines a partial thealgorithmin[36]factorizestheASKITapproximationin D factorizationwithiterativemethods.Asahighlight,weareable O(Nlog2N) time (see §II-B). Roughly speaking, ASKIT . to approximately factorize a dense 11M×11M kernel matrix s in2minuteson3,072x86“Haswell”coresanda4.5M×4.5M is based on the approximation of K as the sum of a c matrix in 1 minute using 4,352 “Knights Landing” cores. block-diagonal matrix and a low-rank matrix followed by [ recursion for each diagonal block. We refer this process 1 as the construction of the hierarchical representation of K. v I. INTRODUCTION Oncewehavethisrepresentation,wecanfactorizeK byap- 4 plyingrecursivelytheSherman-Morrison-Woodbury(SMW) 2 LetX beasetofN points∈Rd andletK(x ,x ):Rd× ¯i ¯j formula.Thefactorizationhastobedonefordifferentvalues 3 Rd →Rbeagivenkernelfunction.Thekernelmatrixisthe of λ during cross-validation studies. Therefore optimizing 2 N ×N matrix whose entries are given by K =K(x ,x ) 0 ij ¯i ¯j the factorization is crucial for the overall performance of for i,j =1,...,N, x ,x ,∈X. . ¯i ¯j a kernel method. In this paper, we extend the factorization 1 Kernel matrices appear in unsupervised and supervised scheme presented in [36] in several ways: 0 statistical learning, Gaussian process, regression, and non- 7 • We present an algorithm that factorizes the ASKIT parametric statistics [8], [12], [27], [33]. Solving linear 1 approximation in O(NlogN) time instead of systems with kernel matrices is an algebraic operation that : O(Nlog2N) and we demonstrate its performance on v is required in many kernel methods. The simplest example i several datasets (see §II-C and §V). X is ridge regression in which we solve λI+K, where λ>0 • We present a hybrid level-restricted factorization is a regularization parameter that controls generalization r schemethatreducesdramaticallythefactorizationtime a accuracyandI istheidentitymatrix.Thislinearsolvecanbe by using a Krylov iterative solver on a much smaller prohibitively expensive for large N because K is typically system than the original (see §II-C). The new method dense. For example, consider the Gaussian kernel, can be used with matrices for which [36] fails. (cid:32) 1(cid:107)x −x (cid:107)2(cid:33) • We study performance on Intel’s “Knights Landing” K(x ,x )=exp − ¯i ¯j 2 , (1) (KNL) architecture. We introduce an optimized ¯i ¯j 2 h2 matrix-free kernel summation that reduces the storage requirementsofthefactorizationwithouthavingavery wherehisthekernelbandwidth.Forsmallh,K approaches significant impact on wall-clock time (see §II-D). the identity matrix whereas for large h, K approaches the rank-one constant matrix. The first regime suggests In our numerical experiments, we measure the performance sparse approximations while the second regime suggests ofthemethodonseveraldifferentdatasets.ASKIThasbeen global low-rank approximations. But for the majority of h applied to polynomial, Matern, Laplacian, and Gaussian values, K is neither sparse nor globally low-rank. Direct factorization of λI +K requires O(N3) work, whereas a 1ApproximateSkeletonizationKernelIndependentTreecode.TheASKITlibrary Krylov iterative method costs O(N2) work per iteration isavailableathttp://padas.ices.utexas.edu/libaskit. kernels in arbitrary dimensions. Due to space limitations on-diagonal blocks K and K are themselves hierarchi- ll rr we only present the factorization results for the Gaussian cal. Note that the low-rank structure is not invariant on kernel function. Also the Gaussian kernel is among the permutations, it very strongly depends on the ordering of hardest to compress in high dimensions. We examine the the columns (or rows since the matrix is symmetric). For performance of the method for different bandwidth ranges notational convenience we write K ≈K(cid:101) =D+UV, where that are relevant to learning tasks, its sensitivity to the U and V are rank s and D is also hierarchical, where use regularization parameter λ, its performance as we increase K(cid:101) to indicate the approximate kernel matrix. the number of points, and its numerical stability (see §III). Inverting hierarchical matrices. When K admits this Limitations. Not all kernel matrices admit a good hier- hierarchicallow-rankapproximation,thenwecanefficiently archical low-rank decomposition. Typically this is related to approximate K−1 using the Sherman-Morrison-Woodbury theintrinsicdimensionalityofthedatasetatdifferentscales. formula along with recursion: SoASKITandsubsequentlyourmethodcanfail.IfASKIT K(cid:101) =D+UV =D(I+WV), and W =D−1U. cancompressthematrix,thesecondpotentialpointoffailure is the choice of regularization parameter. If it is too small, K(cid:101)−1 =(I+WV)−1D−1 (3) our algorithm (as well as [36]) can become numerically =(I−W(I+VW)−1V)D−1 unstable. We can numerically detect the instability, but it =(I−WZ−1V)D−1, where Z =I+VW. is not clear how to fix it while maintaining the log-linear complexity of the algorithm. However, small regularization Recursion is used to invert D−1. After obtaining D−1 we oftenresultsinpoorlearningperformancesothiscornercase compute W = D−1U for a rank-s matrix U and factorize is not important in applications. We discuss this in more the smaller reduced system Z = I + VW ∈ Rs×s. The detail in §III and §V. Also our methods cannot be applied scheme can be easily extended to invert λI+K. tohierarchicaldecompositionsinwhichD issparseandnot To turn this formulation into an algorithm, we need (1) a just block diagonal. For such decompositions, our method method to partition K so that off-diagonal blocks have low- can be used as a preconditioner, as discussed in [36]. rank, (2) an efficient way to compute the low-rank factors Related work. Nystrom methods and their variants [7], U and V, and (3) a scheme to construct the inverse. For [13], [28], [34] can be used to build fast factorizations. the first two tasks we use ASKIT, a method we recently However, not all kernel matrices can be approximated well developed [21]–[23]. ASKIT uses geometric information by Nystrom methods [15], [18], [23], [31]. Factorization (the input points) to permute K by partitioning the points methods based on hierarchical decomposition have been recursively using a binary tree. Interactions between points studied for kernel matrices from points in two or three in a treenode correspond to diagonal blocks of K. In the dimensions [1], [2], [4], [9], but less so in high dimensions recursion,thechildren ofthenodecanbe usedtodefinethe withafewexceptions[15],[36].Earlyworksdiscussingpar- block partitioning of the parent block, similar to (2). Next, allel operations for hierarchical matrices on shared memory we summarize ASKIT features that are necessary for this system include bulk synchronous parallelization [16] and paper. Please see [23] for the complete details on ASKIT. DAG-based task parallelism [17]. Distributed factorization Partitioning the matrix. We use a ball tree [26] to and operations were discussed in [14], [32]. The difficul- partition X. Starting with the root node (which contains the ties of generalizing low-dimensional factorizations in high- entiredataset),nodesarepartitionedintotwochildren(with dimensions are discussed in [21], [36]. an equal number of points) by a splitting hyperplane. This recursive splitting terminates when a node has less than m II. METHODS points, a user-specified parameter. The root has level l = 0 and the leaves l = D = log (N/m), the depth of the tree. 2 Webeginwithasketchofhierarchicalmatricesanddirect In the following, we overload α, β to indicate both binary solvers in §II-A. We also briefly summarize the ASKIT tree nodes and the indices of the points that belong to these algorithm which we use as the basis for our new methods. nodes;|α|isthenumberofpointsinα;andl,rindicatethe We describe parallel factorization schemes in §II-B and left and right children of the node. We define X ∈Rd×N to highlight the novelty of our approach over [36]. We then be the matrix of all points and X to be the points owned α introduce our hybrid iterative/direct solver in §II-C. by tree node α, (i.e., X ={x |∀i∈α}). α ¯i Computing low rank approximations. Let α be the A. Hierarchical Matrices and Treecodes points in a leaf node. Let S = {1,...,N}\α. The skele- Broadly speaking, we consider a matrix K ∈ RN×N to tonization of a node α is a rank-s approximation of K Sα be hierarchical if it can be partitioned as using s columns of K . We refer to these columns as Sα the skeleton of α, denoted by α. Skeletonization is done (cid:20) (cid:21) (cid:20) (cid:21) (cid:20) (cid:21) (cid:101) K K K 0 0 K K = ll lr = ll + lr . (2) using the Interpolative Decomposition (ID) [11]. Using a K K 0 K K 0 rl rr rr rl pivoted rank-revealing QR factorization, the ID finds α and (cid:101) P ∈Rs×|α| such that where the off-diagonal blocks Klr and Krl can be accu- α(cid:101)α rately approximated by a low-rank factorization and the K ≈K P . (4) Sα Sα(cid:101) α(cid:101)α Algorithm II.1 [α(cid:101),Pα(cid:101)α]=Skeletonize(α) Kα−α1 ∈ Rm×m using LAPACK’s GETRF. For an internal [i(cid:101)lf,α]=isSlekaefltheetnonreitzuren(l[)α(cid:101);,[P(cid:101)r,α(cid:101)]α=] =SkIeDl(αe)t;onize(r); nworditeeαou,twKe(cid:101)αnαee=d UDαα,(WI+α,WVαα,Vaαn)daZsα−1. Using (3), we can return [α(cid:101),Pα(cid:101)[(cid:101)l(cid:101)r]] = ID([(cid:101)l(cid:101)r]); (cid:20)K(cid:101)ll (cid:21)(cid:18)I+(cid:20)Pˆl(cid:101)l (cid:21)(cid:20) K(cid:101)lr(cid:21)(cid:19), (7) K(cid:101)rr Pˆr(cid:101)r K(cid:101)rl The first s pivots from the QR define α. Using the QR, we icnanOc(odmNp2umte)Pαc(cid:101)oαm=plKexS†itα(cid:101)yKfSoαr.thTehiosvsecrh(cid:101)aellmfeachtoowriezvaetirorne.suWltes wthheer“ehwate”dneofitnaetioPˆnl)(cid:101)l.=ThK(cid:101)elr−el1foPrle(cid:101)l,aPnˆdαα(cid:101)Pˆr(cid:101)rre=quKi(cid:101)rer−sr1P“irn(cid:101)rv(enrtointigc”e can turn it to a O(dlogNm) scheme by sampling a small K(cid:101)α−α1,whichinturnrequirestraversingallthedescendantsof subsetS(cid:48)ofSandusingitinsteadofS[23].Theapproxima- α(thesubtreerootedatα)andrecursivelyapplying(7).(We introducedthisschemein[36]andresultsinO(dNlog2N) tion rank s is chosen such that σ (K )/σ (K )<τ, s+1 S(cid:48)α 1 S(cid:48)α complexity.)Butaswewillseeshortlythissubtreetraversal where τ is user-specified and σ are the singular values is not necessary.) Once we have W , we use the SMW estimated by the diagonal of the rank-revealing QR. α formula to invert (I + W V )−1. This inverse requires For a non-leaf α, we first compute the skeletons (cid:101)l and (cid:101)r W Z−1V , which in termsαofαthe block decomposition of of the children of α and then we compute the skeleton α(cid:101)⊂ α α α α, can be written as (cid:101)lID∪d(cid:101)rec=om[(cid:101)lp(cid:101)ro]siatniodnt(hAelpgroorjiethcmtioInI.1m)a.tOrinxcPeα(cid:101)th[(cid:101)le(cid:101)r]skueslientgonainzoattihoenr (cid:20)Pˆ (cid:21)(cid:20) I K Pˆ (cid:21)−1(cid:20) K (cid:21) oifeevnetrryynoofdKe(wbubtythKerowot≈)isKcomwp(αut)e+d,(cid:80)wDecaKncoPmpuwte(βth)e, l(cid:101)l Pˆr(cid:101)r K(cid:101)rlPˆl(cid:101)l (cid:101)lrI r(cid:101)r K(cid:101)rl (cid:101)lr . (8) wthhere i∈α and β arie: the sibilαings of α aln=d0itsiβ(cid:101)ancβ(cid:101)eβstors. Since α is the parent of l and r, Pˆ and Pˆ have been In ASKIT and [36], U and V of a node α for (3) are already computed. K Pˆ and K lPˆ(cid:101)l are cr(cid:101)romputed by (cid:101)lr r(cid:101)r (cid:101)rl l(cid:101)l (cid:20) (cid:21) (cid:20) (cid:21)(cid:20) (cid:21) GEMM, and the reduced system is factorized by GETRF. 0 K K P Krl 0lr ≈UαVα = Kr(cid:101)l l(cid:101)r (cid:101)ll P(cid:101)rr . (5) andWPeˆca.nWexepsloayittaha“ttePlescoipsi“ntge”lesrecloaptieodn”bfreotwmePenPˆl(cid:101)l,,PPˆr(cid:101)r, IbnetahpipsrwoxoirmkawteeiunsetwtoheeqfaucitvathleanttKfolrrm=s:KKrTl.PThuosr,PKlrKcan. and Pαr(cid:101)rα(cid:101)because is comαpα(cid:101)uted by formula in the b[o(cid:101)lx(cid:101)r]α(cid:101)belolw(cid:101)l. Here we use the second form to write l(cid:101)r (cid:101)rr l(cid:101)l (cid:101)lr D−1 P =(cid:20)K(cid:101)l−l1 (cid:21) (cid:20)Pl(cid:101)l (cid:21)P . (9) (cid:20)K0 K0lr(cid:21)≈UαVα =(cid:20)Pl(cid:101)l P (cid:21)(cid:20)K K(cid:101)lr(cid:21). (6) α αα(cid:101) K(cid:101)r−r1 Pr(cid:101)r [(cid:101)l(cid:101)r]α(cid:101) rl r(cid:101)r (cid:101)rl The calculation in the box requires just GEMM operations We will see that having the P terms on the left allows us from the children (l and r) but not all descendants. Since toLdeesviegln arenstOri(cNtiolong. NW)efapcrtoorviizdaetiodnetaaliglosritohnm.the level- DPˆα−α(cid:101)1P=K(cid:101)α−wα1iPthαα(cid:101)(9=).(IN+owWαwVeα)fi−n1dDtα−h1aPt αPˆα(cid:101),wecacnanarlespolabcee restriction feature of ASKIT. As we saw in Algorithm II.1, telαescoαpα(cid:101)ed by Pˆ and Pˆ as αα(cid:101) l(cid:101)l r(cid:101)r the skeletonization proceeds in a bottom-up traversal of the bgarollwtirnege.lAarsgewreantrda,vdeerspeenthdeintgreoen,tthheeopfrfo-bdlieamgo,ntahlebnleocceksssaarrye Pˆαα(cid:101) =(cid:18)I+(cid:20)Pˆl(cid:101)l Pˆ (cid:21)Vα(cid:19)−1(cid:20)Pˆl(cid:101)l Pˆ (cid:21)P[(cid:101)l(cid:101)r]α(cid:101). (10) r(cid:101)r r(cid:101)r rank s can increase to the extend that no compression Notice that we no longer need to solve K(cid:101)−1 and K(cid:101)−1 in takes place. To guarantee accuracy, skeletonization of α ll rr (10). Thus, no tree traversal is required. In the leaf level should terminate if α(cid:101) = (cid:101)l∪(cid:101)r. In some of our numerical (base case), Pˆ is computed directly from K−1P . experiments, instead of using this criterion, we use L to αα(cid:101) αα αα(cid:101) Giventheseformulasandtheskeletonizationcomputedin represent the level at which the skeletonization stops. AlgorithmII.1,wecomputethefactorsneededforthedirect With level restriction, the factorization described in [36] solverinapostordertraversalofthetree(AlgorithmII.2).If cannot be used. We introduce a hybrid iterative/direct αisaleafnode,wefactorizeλI+K usinganLUfactor- scheme that addresses this shortcoming in §II-C. αα ization. Otherwise, we compute K and K . Notice that B. Fast direct solver Pˆ andPˆ arecomputedinthepre(cid:101)vlriousrecu(cid:101)rrlsion;thus,we We have sketched how we compute the UV approxima- caln(cid:101)l formar(cid:101)rndfactorizethereducedsystemZ .Finally,Pˆ α αα(cid:101) tions in K(cid:101) using ASKIT. We now discuss how to use them is telescoped using (10), thus Solve(α,WαP[(cid:101)l(cid:101)r]α(cid:101),false) inthecontextof(3)tosolveλI+K directly.Forsimplicity, (Algorithm II.3) will not invoke recursion. we describe the case where λ = 0, but all the algorithms This algorithm improves on the one in [36] by removing we describe trivially generalize to the λ (cid:54)= 0 case. We first the extra subtree traversals that result in O(Nlog2N) com- consider the case in which no level restriction takes place. plexity. Instead, our algorithm exploits the nested structure We assume that all the internal nodes have been skele- of Pˆ resulting in an NlogN complexity for the factor- αα(cid:101) tonized. The factorization of K(cid:101) proceeds using a bottom- ization.Insomeofourlargestruns,thisresultedinover3× up traversal of the tree. At the leaf level, we factorize speedup without any change in the accuracy. Algorithm II.2 Factorize(α) Algorithm II.4 DistFactorize(α,q) if α is leaf then if α is at level logp then Factorize(α). LU factorization λI+K . else αα W =P and P =I. DistFactorize(c, q). α αα(cid:101) [(cid:101)l(cid:101)r]α(cid:101) 2 else {i< q} {i≥ q} Pˆαα(cid:101)FFLo=aUrcmSftaoocWltrovαiriezzw(eeαit(th,lhW)ePâαrlne(cid:101)lPdd,[u(cid:101)lPFˆ(cid:101)rcr]eaα(cid:101)(cid:101)rdc,,ftasynaosdlrtesiVmezα)eZw(uαrist)ih.inngK((8(cid:101)l1r)0.,)K. (cid:101)rl. {{R{000e}}}dSRuRee2cencecdvv,K(cid:101)lB(cid:101)r.c{xa}sPˆt{x(cid:101)r}(cid:101)l.. {{{Re22qqq}}}dRSuSee2cecnenvddK,(cid:101)rBK(cid:101)l.{cxa}sPˆt{Pˆx(cid:101)l}(cid:101)r.. . 2 (cid:101)l{x} {x}(cid:101)r {0} Bcast P and LU factorizes Z. [(cid:101)lr(cid:101)]α(cid:101) Algorithm II.3 w=Solve(α,u,do_recur) Pˆ =DistSolve(α,Pˆ P ,q,false). ¯ ¯ {x}α(cid:101) {x}(cid:101)c [(cid:101)c]α(cid:101) if α is leaf then LU solver w=(λI+K )−1u αα ¯ ¯ else Algorithm II.5 w=DistSolve(α,u,q,do_recur) if do_recur then ¯ ¯ v = [Solve(l,u ,true); Solve(r,u ,true)]. ifαisatlevellogpthenw =Solve(α,u,do_recur). ¯l ¯r ¯ Compute w=u−W Z−1V u using (8). else ¯ α α α if do_recur then u=DistSolve(c,u,q,do_recur). ¯ 2 {i< q} {i≥ q} 2 2 P𝛂𝛂~K𝛂~β Reduce K(cid:101)r{x}u. Reduce u(cid:101)l =K(cid:101)l{x}u. Pβ~βKβ~𝛂PKr~rKllr~Kl~πP βKl~lKrrl~rPβ θ~ Pθ~θK~θπ 𝛂 θ β lleevveell--12 π {{{{0000}}}} RS[Bueec(cid:101)lcna;uvds(cid:101)rt]uu=(cid:101)(cid:101)lru.. Z. −1[u(cid:101)l;u(cid:101)r]{{{u22qqqsi}}}ngSRBeec(8nca)dvs.tuu(cid:101)(cid:101)lru.. . l r w =u−Pˆ (cid:101)lu . w2=u−Pˆ (cid:101)ru . {x}(cid:101)l (cid:101)l {x}(cid:101)r (cid:101)r Pππ~Kπ~θ level-3 Distributed Local nestedconstructisenabledsuchthateachthreadwillinvoke Tree Tree parallelBLASorLAPACKroutines.Ifwehavepprocesses, Figure1ThetopfourlevelsofthetreeandthecorrespondingblocksofthematrixK˜. Thenodesbelongingtoeachprocessarehighlightedinasinglecolor.Eachprocess then above level logp of the tree, we have to communicate factorizesitsownportionofthetreeindependently.Wealsohighlightthefactorsused to compute factors, since the terms needed are distributed in thedirect solverconstruction andshow whichprocess owns whichfactor. Each processownadiagonalblockandallfactorsinthesamecolumnandthesamerow. among processes. We use the Message Passing Interface Forexample,theyellowprocessownsPββ(cid:101) andKα(cid:101)β atlevel1;similarlyitowns (MPI) library for distributed memory communication. Pβθ(cid:101)andKπ(cid:101)β atlevel0. In Figure 1, we summarize the distributed-memory algo- WethendescribehowtoapplyK(cid:101)α−α1 toavectoru¯,shown rithm. The four colors represent four different MPI ranks inAlgorithmII.3.Ifαisaleafnode,wecandirectlyinvoke and the nodes they own. The tree on the right shows l =2 an LU solver to obtain w¯ = Kα−α1u¯. Otherwise, we have treenodesareuniquelyassignedtoranks,buttreenodeswith two situations. If Algorithm II.3 is called by Factorize l > 2 are shared among ranks. To facilitate collective (Tdhou_s,rneocurercuirssiofnalissere)q,utihreend.wWehiklenoSwoluve=isWcαalPle[(cid:101)ld(cid:101)r]αt(cid:101)o. ccoommmmuunniiccaattioorn,,wheaicchheqduisatlrliybudtievdidetrsetehneodraenkcsreoaftethseapalroecnatl. solvearandomu,thenweneedtosolveK(cid:101)−1u andK(cid:101)−1u We use {i} to denote the i MPI rank in the local commu- ll ¯l rr ¯r th recursively and compute (8) with GEMV on Vα and Wα and nicator. Consider the communicator of α, which involves q an LU solve GETRS on Zα. ranks.Letcdenotethechildofαthat{i}owns.Thenc=l Therefore,thecompletealgorithmconsistsofconstructing if{i< q}.Otherwise,c=r.Foradistributednodeα,data 2 the tree, calling Algorithm II.1, then Algorithm II.2 and points {x} owned by {i} are never required by other MPI i Algorithm II.3, each called on the root of the tree. processes, and Xα =Xl∪Xr =(∪i λ, κ(λI +D) can grow even as and employ GMRES CGS refinement. If not specified, κ(λI + K) remains small. If this case happens in many KNL experiments use Cache-Quadrant configuration with levels of our factorization, then the method is not stable. OMP_PROC_BIND=spread.“T”referstothetotalruntime With narrow bandwidths where K approaches a (blocked)- in seconds, and “GFs” refers to the GFLOPS per node. diagonal matrix, σn >λ may occur. Datasets. We use real-world datasets: COVTYPE (forest Undertheframeworkofhierarchicalmatrices,thepivoting cartographic variables); SUSY and HIGGS (high-energy rowswecanchooseduringtheDfactorizationarelimitedto physics) [19]; MNIST (handwritten digit recognition) [5]; theskeletonrows.Thus,evenκ(λI+K)isnotbad,(λI+D) and MRI (brain MRI) [25]. We also use a 64D synthetic can be unstable due to the aggressive pivoting strategy if λ dataset, which is drawn from a 6D Normal distribution and is small. Our methods can detect this situation, but avoiding embeddedin64Dwithadditionalnoise.Thissetisadataset this case entirely (or fixing it) is not straightforward. withahighambientbutrelativelysmallintrinsicdimension. Accuracy metrics and parameter selection. For the linear solve, we report the relative residual IV. EXPERIMENTALSETUP (cid:15) =(cid:107)u−(λI+K˜)w(cid:107) /(cid:107)u(cid:107) . (15) We performed numerical experiments on Haswell and r 2 2 KNL architectures with four different setups to examine The parameters h and λ used in the Gaussian kernel were the accuracy and efficiency of our methods. Especially, we selectedusingcross-validation.InTableIIwereportthepa- wanttodemonstrate(1)thecomplexityimprovementagainst rameters we used. Other combinations in §V are candidates [36],(2)FLOPSefficiency,(3)scalabilityand(4)theadvan- for the cross-validation. Level restriction L is choosen such tages of our hybrid solver. We explore the task of kernel that the relative error is controlled. In Table IV, V we use ridge regression for binary supervised classification [33], L=3, and for experiments in Figure 5 we use L=5 or 7. which requires approximating the solution of (λI +K)−1 during the training step. We use the Gaussian kernel with V. EMPIRICALRESULTS bandwidth h. The model weight w is chosen by solving The experiments are labeled #1 to #39 in the tables and w =(λI+K(cid:101))−1u, where u is given (the labels). Once w is figures.Weselectrepresentativeparametercombinationsand computed, the label given by x∈/ X is sign(K(x,X)w). We comparetheruntimebetween[36]andourAlgorithmII.4in ¯ ¯ ¯ applyourmethodstotrainthismodelonreal-worlddatasets Table III. We present single node performance on Haswell employingupto3,072x86coresand4,352KNLcores.The and KNL in Table IV. In Figure 4, we present strong the percentage of correct predictions (Acc) is reported in scalability and verify the O(NlogN) complexity of our Table II, along with the optimal h and λ that were found methods (Algorithm II.4). In Table V, we compare our using holdout cross validation. Implementation and hardware. Our experiments were 3WeestimatethepeakaccordingtotheclockrateandthetheoreticalFMAthrough- conducted on Lonestar5 (two 12-core, 2.6GHz, Xeon E5- put. For 24 Haswell cores, 998 = 2×12×2.6×16. For 68 KNL cores, 3046 = 68×1.4×32. As a reference, MKL GEMM can achieve 87% on the 2690 v3 “Haswell” per node) and Stampede (68-core, Haswell node and 69% on the KNL node. We assume two VPUs can dual issue 1.4GHz, Xeon Phi 7250 “KNL” per node) clusters at DFMAs [29]. However, Intel processors may have a different frequency while fully issuingFMA,andtheclockratemaydropto1.0GHz.Thismaybethereasonwhy the Texas Advanced Computing Center. The theoretical MKLGEMMcanonlyachieve2.1TFLOPSonKNL. τ 1E-1 1E-3 1E-5 # 11 12 13 14 15 16 # dataset h log2 log log2 log log2 log Config Haswell Cache-Quad F-Quad F-SNC4 1 COVTYPE .35 <1 <1 5 3 15 8 p 1 4 1 4 1 4 2 .07 5 3 24 11 27 12 nthd 24 6 68 17 68 17 3 SUSY .50 <1 <1 1 <1 6 3 Tf 89 94 41 52 62 77 4 .05 63 23 125 37 125 37 GFf 623 592 1356 1136 895 723 5 MNIST2M 1.0 24 11 38 15 40 16 ComputeMatVecV withGEMV. 6 .10 1 <1 1 <1 4 7 Ts 0.8 0.8 0.9 0.9 0.9 1.0 7 MNIST8M 0.3 142 51 185 61 195 64 GFs 14 14 12 12 12 11 8 HIGGS 2.0 75 29 324 84 326 85 ReevaluateV withGEMM. 9 .90 216 66 317 83 323 85 Ts 4.4 5.3 7.4 4.0 7.5 5.8 10 NORMAL32M .19 18 10 46 22 84 28 GFs 158 119 40 74 39 51 ComputeMatVecV withGSKS. TLoanbelestaIrI5IuFsiancgto1r2iz8atnioondetsim(3e,0c7o2mcpoarreiss)onwitihnasdeacpotnidve. rEaxnpkesrismesnetlsectaerdebdyonτeaonnd Ts 1.1 1.8 1.1 1.2 1.2 1.3 smax.COVTYPEusedm=2,048,κ=2,048,smax=2,048.SUSYusedthe GFs 269 164 269 243 243 228 samecombination.MNISTusedκ=256.HIGGSusedm=512,κ=1,024. TableIVSinglenodeperformanceonCOVTYPE100Kwithm=s=2048(fixed NORMALusedm=512,κ=128andsmax=256. rank) and L = 3. p is the number of MPI processes, and nthd is the number hybrid (Algorithm II.6) with the direct method (Algo- of OpenMP threads per process. #11 and #12 were conducted on Haswell, and #13 to #16 where conducted on KNL with different configurations. #13 and #14 rithm II.5). In Figure 5, we report the convergence behavior use cache-quadrant. #15 uses flat-quadrant, and #16 uses flat-snc-4 configuration. ofiterativesolver(AlgorithmII.6)onλI+K(cid:101) andourhybrid #11 achieves the high factorization performance on Haswell, and #13 is the most efficientconfigurationonKNL.Ineachcolumn,wereportfactorizationtime(Tf)and factorization. Here the parameter τ indicates the relative GFLOPS(GFf)andthreedifferentsolvingtime(Ts).Thesethreesolvingschemes tolerance of approximation of the kernel matrix K, s havedifferentstoragerequirementandmemoryoperations. max is the maximum skeleton size, κ is the number of nearest neighbors used for skeletonization sampling in ASKIT, m Strong scaling Experiment Experiment is the leaf node size, and L is the level restriction. 101 IIddeeaall NNllooggN2N 32M 102100% 100% Ideal scaling 91% facCtoormizaptaiornisotinmewbitehtw[e3e6n]th(eTaOb(lNe lIoIgI)2.NW)ealcgoomritphamre[3th6e] e (sec) 16M e (sec) 93%90% 85% and our O(NlogN) algorithms using the same parameters. Tim100 4M Tim 81% 82% We only compare the case without level restriction, because 1M NORMAL 70% [36] does not support this feature. The runtime is directly 64D h0.15 101 74% #17 associated with the rank s. For example, the U, V matrices 10-1 106 107 102 103 62% in #4 are much larger than those in #3. Thus, the runtime N #core is also much longer. The overall speedup is about 2–4× Figure4O(NlogN)verification(#17)andstrongscaling(#18).#17useNORMAL with m = 512, κ = 128, a fixed s = 256 and L = 1. The blue lines are due to the logN term. Both methods construct exactly the experimentalfactorizationtime,andyellowlinesaretheoretical(ideal)time.#18use same factorization (up to roundoff errors). Although the NORMAL1Mwithκ=128,m=s=2048andL=1.Weincreasenumberof nodesupto128Haswellnodes(3,072cores)and64KNLnodes(4,352cores).The speedupisnotexactlylogN duetotheprefactorsdepending greenlinesrepresenttheidealscaling. on s and d, we can expect the asymptotic speedup to be O(dsNlogN) where d is 54 in this dataset. GEMM takes O(logN).Forexample,COVTYPEcanonlyachieve1.9×, O(sN) to store the matrix, but GSKS is matrix free with but HIGGS can achieve 3.8× because the problem size O(1)space.GSKSisonly1.2−1.6×slowerthantheGEMV is 20× larger. Both methods have NlogN complexity for approach while it is 4−7× faster than the GEMM approach. the “Solve” operation. For the experiments in Table III, the Notice that 80% of T is dominated by GETRS (triangular longest “Solve” operation (#8) is less than 2 seconds. s solver) when d is small. Since GETRS can only achieve 2 Single node performance (Table IV). We conduct a set GFLOPSonaKNLnode,theoverallGF issomewhatlow. s of single node experiments to show the FLOP rates we Scaling (Figure 4). In #17, we use 128 nodes (3,072 achieve and to test some of the memory models on KNL. cores) and increase N from 1M to 32M. We can observe On a Haswell node, #11 reaches 62% (623/998) of the that our implementation is very close to the theoretical theoretical peak. For a KNL node, #13 (Cache-Quadrant) NlogN scaling (yellow) but lower than the Nlog2N is the fastest and achieves 45% (1356/3046). Using MPI scaling (purple). In #18, we fix the data set (NORMAL on KNL (p > 1) is typically slower due to the extra 1M) and increase the number of cores. The green line is memory operations. Our implementation does not perform the ideal scaling (100%), and our implementation reaches very well on the flat memory mode (#15 and 16). The 62% efficiency on 3,072 Haswell cores and 70% on 4,352 memory requirements usually exceed 16GB and as a result KNL cores. This relatively small problem (1M) cannot U and V cannot fit into MCDRAM. We tried manually fully exploit all computing resources; thus, we can see the swapping memory between MCDRAM and DDR4 but this degradation while N is small or when the number of cores was not as efficient as using the cache memory mode. is large (∼230 points per core for 64 KNL nodes). Reducing storage (Table IV). In Table IV, we report Hybrid and direct methods comparison (Table V). three different schemes for kernel summation §II-D: GEMV, In Table V we set L = 3 (level restriction) and com- GEMM and GSKS. The first scheme takes O(sNlogN) pare Algorithm II.6 (hybrid) and Algorithm II.2 (direct) time and space. The last two schemes evaluate K u in for problems that Algorithm II.2 can be applied (i.e., the β(cid:101)α 122#901 KNLx-- ASK255I974T404SUS1TY651090fh=G8440F416.447f15λ2112T=...523s40G2F894296s 556ee--e11(cid:15)-r422 KS9P8-- 11111010000-01----086420 λκCh# 2 O0318.eV0+T72YPE 1111100000----08642 λκCh# 2 O0019..eV03+T73YPE 11110000--0121 λκCh# 3 O0010..eV00+T705Y3PE MRIh=3.5λ=10 100 150 200 250 300 100 200 300 400 500 100 200 300 400 500 222234 x-- 332091267 483647 744926577 3119...456 325190988 113ee--e11-400 93-- 111000--042 λSh U04S.010Y5 111000--042 λSh U04S.01Y5 111000--021 λSh U00S..14Y5 MNIST2Mh=1.0λ=1 1100--86 κ#3 11e+2 1100--86 κ#3 12e+3 1100--43 κ#3 13e+5 25 x 237 27 655 1.9 592 1e-13 - 10-5 200 300 400 500 600 200 300 400 500 600 200 300 400 500 600 26 - 270 47 372 2.3 489 1e-13 - 27 - 217 19 404 27.2 612 1e-3 27 100 HIGGS 100 HIGGS 100 HIGGS 10-2 λ 0.08 10-2 λ 0.008 10-2 λ 0.00008 TableVHybridanddirectmethodscomparisonwithelevelrestrictionL=3.All 10-4 h 0.9 10-4 h 0.9 10-4 h 0.9 reexppoerrtimAeSntKsIuTsebaudiladpitnigvetirmanek,sfawcittohritzoaletiroannctiemτe=(T0f).,0s0o0lv0i1ngantidmsemeamxpl=oyi2n0g4G8S.KWSe 1100--86 κ#3 14e+2 1100--86 κ#3 15e+3 1100--86 κ#3 16e+5 (Ts) and efficiency in GFLOPS. The three experiments with orange index are the 500 600 700 800 900 500 600 700 800 900 500 600 700 800 900 hybridschemeandtheremainingarethedirectfactorization.Wereporttherelative residual(cid:15)r andnumberofKryloviterations(KSP)(forthehybridmethod). 1100-02 λM N0.I0S2T 1100-02 λM N0.I0S0T2 1100-02 λM N0.I0S0T002 full factorization requires 2LsN +22Ls2 memory for level 10-4 h 0.1 10-4 h 0.1 10-4 h 0.1 10-6 κ 1e+2 10-6 κ 1e+3 10-6 κ 1e+5 restriction L), with adaptive s selection. #20, #23 and #26 10-8 #37 10-8 #38 10-8 #39 use the direct factorization on Haswell, and #19, #22, and 40 45 50 55 60 40 45 50 55 60 40 45 50 55 60 #25onKNL.TheremainingthreerunsaredoneonHaswell Figure 5 solvingλI+K(cid:101):Convergenceoftherelativeresidual(cid:15)r (verticalaxis) using the hybrid method. On Haswell, we observe that the over time (horizontal axis) in seconds. (a) (blue) is unpreconditioned GMRES; (b) factorization time Tf is about two times longer than #21, a(orreatnhgee)saimsehyabsriTdabmleethIIoId;.κAilsl ethxepecroimndeintitosnusneudmbτer=.CO1EVT−YP5E.amnd, kHIaGnGdSsmusaexd #24 and #27. In the factorization phase, both Haswell and L=5restriction.SUSYandMNISTusedL=7.Theserunsmaybe100–1000× moreexpensive(alsorunningoutofmemory)ifweweretousethedirectsolver. KNL do not perform as well as in Table IV, because using more predictable behavior. #30 is detected numerically ill- adaptive ranks s results in load imbalance. If we further conditioningofD inoursolver.Alsonoticethatin#30both increaseL,asweneedtodoinFigure5,thecostofthefull methods fail to converge. factorization can be 1000× in runtime and 30× in storage. Applying the hybrid solver to a vector is slower than VI. CONCLUSIONS applying the direct solver, due to the need of iteration. E.g., We have introduced new algorithms for approximately T in #21 is about 20× slower. Yet the overall runtime s factorizing kernel matrices. We evaluated our algorithms (T + T ) of the hybrid method is still smaller than the f s on both real-world and synthetic datasets with different direct one. When L is larger, the advantage of the hybrid parameters.Weconductedanalysisandexperimentstostudy solver will be higher. For example, in Figure 5, we report the complexity and the scalability of our methods. These results that require L = 7. Algorithm II.2 cannot be used: experiments include scaling up to 3,072 Haswell cores and the memory just for Z with s=2048 exceeds 500GB. 4,352 KNL and exhibit significant speedups over existing Convergence behavior for solving λI +K(cid:101) (Figure 5). methods. The factorization can be very fast. For example, We report the convergence rate using four different band- it only takes 10 seconds to factorize a kernel matrix with widths with two different methods: (a) unpreconditioned 32M points in 64D. Our future work will focus on further GMRES using ASKIT’s MatVec for λI + K(cid:101) (blue line) optimizationof ourimplementation. Inparticular, wewould and (b) our hybrid method Algorithm II.6 (orange line). like to introduce task parallelism in the tree traversal to Each row corresponds to a dataset with a specific h. These address the load balancing issue. While adaptive ranks or experiments resemble a cross-validation study in which we adaptive level restriction is used, each treenode may have vary λ in order to improve learning. Across columns, we different workload. In this case, scheduling is important to vary λ as [10−2,10−3,10−5]σ1(K(cid:101)), where σ1(K(cid:101)) is an avoid the critical path. Additionally, we plan to address the estimate, so that the condition number κ of λI+K(cid:101) is 102, stability issues mentioned in §III and explore other possible 103 and 105 respectively. We report the relative (to a zero variants (e.g. sparse off-diagonal blocks). initial guess) Krylov residual (cid:15) (y-axis) over time (x-axis). r The steeper the curve is, the faster the method converges. The x-axis offset represents setup costs. E.g., in #28, the offsetoftheblueline(≈140sec)isthecostofbuildingthe tree and the skeletons (spent in ASKIT). The fixed cost of (b) includes the fixed cost of (a) plus the factorization time. We see that most of the blue lines are flat when the condition number is around 1E+5, but orange lines still decrease steadily except for #30 (see §III for the stability issue). We can observe 10–1000× speedup on the “Solve” operations. Overall the hybrid scheme is faster and has REFERENCES [19] M.LICHMAN, UCI machine learning repository, 2013. [1] S.AMBIKASARANANDE.DARVE,AnO(nlogn)fastdirect [20] Z.LUETAL.,Howtoscaleupkernelmethodstobeasgood solver for partial hierarchically semi-separable matrices, asdeepneuralnets,arXivpreprintarXiv:1411.4000,(2014). Journal of Scientific Computing, 57 (2013), pp. 477–501. [21] W. B. MARCH, B. XIAO, AND G. BIROS, ASKIT: Approxi- [2] S. AMBIKASARAN, D. FOREMAN-MACKEY, L. GREEN- mate skeletonization kernel-independent treecode in high di- GARD, D. W. HOGG, ANDM. O’NEIL, Fast direct methods mensions,SIAMJournalonScientificComputing,37(2015), for Gaussian processes and the analysis of NASA Kepler pp. 1089–1110. mission data, arXiv preprint arXiv:1403.6015, (2014). [22] W. B. MARCH, B. XIAO, S. THARAKAN, C. D. YU, AND [3] S.BALAYETAL., PETSc Web page, 2014. G.BIROS,Akernel-independentFMMingeneraldimensions, in Proceedings of SC15, Austin, Texas, Nov 2015. [4] M.BEBENDORF, Hierarchical matrices, Springer, 2008. [23] , Robust treecode approximation for kernel machines, [5] C.-C.CHANGANDC.-J.LIN,Libsvm:Alibraryforsupport in Proceedings of the 21st ACM SIGKDD Conference on vector machines, ACM Transactions on Intelligent Systems Knowledge Discovery and Data Mining, Sydney, Australia, and Technology (TIST), 2 (2011), p. 27. August 2015, pp. 1–10. [6] B. DAI, B. XIE, N. HE, Y. LIANG, A. RAJ, M.-F. F. [24] W. B. MARCH, B. XIAO, C. D. YU, AND G. BIROS, An BALCAN, AND L. SONG, Scalable kernel methods via dou- algebraic parallel treecode in arbitrary dimensions, in Pro- bly stochastic gradients, in Advances in Neural Information ceedings of IPDPS15, Hyderabad, India, May 2015. Processing Systems, 2014, pp. 3041–3049. [25] MENZEETAL.,TheMultimodalBrainTumorImageSegmen- [7] A. GITTENS AND M. MAHONEY, Revisiting the Nystrom tation Benchmark (BRATS), IEEE Transactions on Medical method for improved large-scale machine learning, in Pro- Imaging, (2014), p. 33. ceedings of ICML13, 2013, pp. 567–575. [26] S. M. OMOHUNDRO, Five balltree construction algorithms, International Computer Science Institute Berkeley, 1989. [8] A. GRAY AND A. MOORE, N-body problems in statistical learning,Advancesinneuralinformationprocessingsystems, (2001), pp. 521–527. [27] C.E.RASMUSSENANDC.WILLIAMS, Gaussian Processes for Machine Learning, MIT Press, 2006. [9] W. HACKBUSCH, B. N. KHOROMSKIJ, AND E. E. TYR- [28] S.SI,C.-J.HSIEH,ANDI.DHILLON, Memory efficient ker- TYSHNIKOV, Approximate iterations for structured matrices, nel approximation, in Proceedings of The 31st International Numerische Mathematik, 109 (2008), pp. 365–383. Conference on Machine Learning, 2014, pp. 701–709. [10] W. W. HAGER, Updating the inverse of a matrix, SIAM [29] A.SODANIETAL.,Knightslanding:Second-generationintel review, 31 (1989), pp. 221–239. xeon phi product, IEEE Micro, 36 (2016), pp. 34–46. [11] N. HALKO, P.-G. MARTINSSON, AND J. TROPP, Finding [30] F.G.VANZEEANDR.A.VANDEGEIJN,Blis:Aframework structurewithrandomness:Probabilisticalgorithmsforcon- forrapidlyinstantiatingblasfunctionality,ACMTransactions structingapproximatematrixdecompositions,SIAMReview, on Mathematical Software (TOMS), 41 (2015), p. 14. 53 (2011), pp. 217–288. [31] R.WANG,Y.LI,M.W.MAHONEY,ANDE.DARVE, Struc- [12] T. HOFMANN, B. SCHO¨LKOPF, AND A. J. SMOLA, Kernel tured block basis factorization for scalable kernel matrix methodsinmachinelearning,Theannalsofstatistics,(2008), evaluation, arXiv preprint arXiv:1505.00398, (2015). pp. 1171–1220. [32] S.WANG,X.S.LI,J.XIA,Y.SITU,ANDM.V.DEHOOP, [13] C.-J.HSIEH,S.SI,ANDI.S.DHILLON, Fast prediction for Efficientscalablealgorithmsforsolvingdenselinearsystems large-scalekernelmachines,inAdvancesinNeuralInforma- with hierarchically semiseparable structures, SIAM Journal tion Processing Systems 27, 2014, pp. 3689–3697. on Scientific Computing, 35 (2013), pp. C519–C544. [14] M.IZADIKHALEGHABADIETAL., Parallel H-matrix arith- [33] L. WASSERMAN, All of Statistics: A Concise Course in metic on distributed-memory systems, (2012). Statistical Inference, Springer, 2004. [15] R. KONDOR, N. TENEVA, AND V. GARG, Multiresolu- [34] C.WILLIAMSANDM.SEEGER, Using the Nystro¨m method tion matrix factorization, in Proceedings of ICML14, 2014, tospeedupkernelmachines,inProceedingsofthe14thAn- pp. 1620–1628. nual Conference on Neural Information Processing Systems, no. EPFL-CONF-161322, 2001, pp. 682–688. [16] R.KRIEMANN,ParalleleAlgorithmenfrH-Matrizen,disser- tation, Universitt Kiel, 2004. [35] C.D.YU,J.HUANG,W.AUSTIN,B.XIAO,ANDG.BIROS, Performance optimization for the k-nearest neighbors kernel [17] , H-LU factorization on many-core systems, Preprint, on x86 architectures, in Proceedings of SC15, ACM, 2015, Max Planck Institute for Mathematics in the Sciences, 2014. pp. 7:1–7:12. [18] D. LEE, A. GRAY, AND A. MOORE, Dual-tree fast gauss [36] C. D. YU, W. B. MARCH, B. XIAO, AND G. BIROS, INV- ASKIT: a parallel fast direct solver for kernel matrices, in transforms,Advances inNeuralInformation ProcessingSys- Proceedings of the IPDPS16, Chicago, USA, May 2016. tems, 18 (2006), pp. 747–754.

An $N \log N$ Parallel Fast Direct Solver for Kernel Matrices

Similar An $N \log N$ Parallel Fast Direct Solver for Kernel Matrices