[Paper Review] Compression, inversion, and approximate PCA of dense kernel matrices at near-linear computational complexity

Florian Schäfer, T. J. Sullivan+1 more|arXiv (Cornell University)|Jun 7, 2017

Sparse and Compressive Sensing Techniques4 citations

TL;DR

This paper proposes a near-linear complexity algorithm for compressing, inverting, and performing approximate PCA on dense kernel matrices arising from smooth covariance functions. By identifying a sparse sparsity pattern via spatial configuration and applying zero-fill-in incomplete Cholesky factorization, it achieves ε-approximate matrix compression and inversion in O(N log²N log²d(N/ε)) time and O(N log N logd(N/ε)) space, with provable exponential accuracy and optimal convergence rates for PCA.

ABSTRACT

Dense kernel matrices $\\Theta \\in \\mathbb{R}^{N \ imes N}$ obtained from point evaluations of a covariance function $G$ at locations $\\{ x_{i} \\}_{1 \\leq i \\leq N} \\subset \\mathbb{R}^{d}$ arise in statistics, machine learning, and numerical analysis. For covariance functions that are Green's functions of elliptic boundary value problems and homogeneously-distributed sampling points, we show how to identify a subset $S \\subset \\{ 1 , \\dots , N \\}^2$, with $\\# S = O ( N \\log (N) \\log^{d} ( N /\\epsilon ) )$, such that the zero fill-in incomplete Cholesky factorisation of the sparse matrix $\\Theta_{ij} 1_{( i, j ) \\in S}$ is an $\\epsilon$-approximation of $\\Theta$. This factorisation can provably be obtained in complexity $O ( N \\log( N ) \\log^{d}( N /\\epsilon) )$ in space and $O ( N \\log^{2}( N ) \\log^{2d}( N /\\epsilon) )$ in time, improving upon the state of the art for general elliptic operators; we further present numerical evidence that $d$ can be taken to be the intrinsic dimension of the data set rather than that of the ambient space. The algorithm only needs to know the spatial configuration of the $x_{i}$ and does not require an analytic representation of $G$. Furthermore, this factorization straightforwardly provides an approximate sparse PCA with optimal rate of convergence in the operator norm. Hence, by using only subsampling and the incomplete Cholesky factorization, we obtain, at nearly linear complexity, the compression, inversion and approximate PCA of a large class of covariance matrices. By inverting the order of the Cholesky factorization we also obtain a solver for elliptic PDE with complexity $O ( N \\log^{d}( N /\\epsilon) )$ in space and $O ( N \\log^{2d}( N /\\epsilon) )$ in time, improving upon the state of the art for general elliptic operators.

Motivation & Objective

• To address the O(N³) computational bottleneck in dense kernel matrix operations common in statistics, machine learning, and PDEs.
• To develop a method that enables near-linear complexity inversion and compression of dense kernel matrices without requiring analytic forms of the kernel.
• To provide a provably accurate, data-adaptive algorithm that automatically exploits intrinsic low-dimensional structure in the data.
• To achieve optimal convergence rates for approximate PCA and fast solvers for elliptic PDEs using only subsampling and incomplete Cholesky factorization.
• To generalize operator-adapted wavelets and gamblet transforms by removing the need for vanishing moment conditions and local averaging.

Proposed method

• Identifies a sparse sparsity pattern S ⊂ {1,…,N}² with size O(N log N logd(N/ε)) based solely on the spatial configuration of sampling points.
• Applies zero-fill-in incomplete Cholesky factorization to the sparse matrix Θ_ij 1_{(i,j)∈S}, which provably approximates the full dense kernel matrix Θ within ε in the operator norm.
• Leverages theoretical results on operator-adapted wavelets and the screening effect to justify the sparsity pattern selection without requiring analytic expressions of the kernel G.
• Uses the Cholesky factor L to directly compute an approximate sparse PCA with optimal convergence rate in the operator norm.
• Reverses the elimination order of the Cholesky factorization to construct a fast direct solver for elliptic PDEs with complexity O(N log²d(N/ε)) in time and O(N logd(N/ε)) in space.
• Employs a hierarchical, data-driven approach that adapts to the intrinsic dimension d of the data set, not the ambient dimension.

Experimental results

Research questions

• RQ1 Can dense kernel matrices from smooth covariance functions be compressed and inverted in near-linear time without requiring analytic forms of the kernel?
• RQ2 Does a sparse Cholesky factorization of a carefully selected sparsity pattern yield an ε-approximation of the full dense kernel matrix with provable accuracy?
• RQ3 Can the incomplete Cholesky factorization of such a pattern provide an approximate PCA with optimal convergence rate in the operator norm?
• RQ4 Can the same factorization be reversed to yield a fast direct solver for elliptic PDEs with improved complexity over existing methods?
• RQ5 Is the algorithm robust to high ambient dimensionality while automatically exploiting the intrinsic dimension of the data?

Key findings

• The algorithm achieves ε-approximate compression and inversion of dense kernel matrices in O(N log²N log²d(N/ε)) time and O(N log N logd(N/ε)) space, significantly improving upon the O(N³) standard Cholesky factorization.
• The sparse Cholesky factorization of the selected sparsity pattern S provides an approximate PCA with optimal convergence rate in the operator norm.
• The method does not require knowledge of the analytic form of the kernel G, only the spatial configuration of the sampling points.
• Numerical evidence suggests that the exponent d in the complexity can be replaced by the intrinsic dimension of the data set, not the ambient dimension.
• The algorithm enables a fast direct solver for elliptic PDEs with O(N log²d(N/ε)) time complexity and O(N logd(N/ε)) space complexity, advancing the state of the art.
• The approach generalizes the gamblet transform by removing the need for vanishing moment conditions and local averaging, enabling broader applicability.