PaperlayerView exact revision on arXiv

Source-linked AI summary

Fast Direct Methods for Gaussian Processes

Sivaram Ambikasaran, Daniel Foreman-Mackey, Leslie Greengard, David W. Hogg, Michael O'Neil

arXiv:1403.6015v2math.NAastro-ph.IMmath.ST

TL;DR

Large-scale Gaussian-process inference is limited by O(n^3) dense covariance inversion and determinant evaluation. The paper uses HODLR factorizations with block low-rank updates to obtain faster direct calculations, including O(n log^2 n) inversion and O(n log n) determinants. The approach supports direct high-dimensional probabilities and related marginalization and hyperparameter tasks under suitable kernel conditions.

  • Problem

    Dense covariance matrices make inversion and determinant evaluation O(n^3), limiting large-scale Gaussian-process inference and exact likelihood calculations.

  • Method

    The method hierarchically factors covariance matrices into products of block low-rank updates to the identity using fast low-rank approximations of off-diagonal blocks.

  • Results

    The factorization requires O(n log^2 n) work for inversion and enables determinant evaluation in O(n log n) operations.

  • Takeaways & Limitations

    The algorithms provide fast, accurate, nearly optimal direct operations for determinants, inverses, and matrix-vector products in Gaussian-process covariance problems.

  • Takeaways & Limitations

    Performance degrades for data in dimensions d > 3 and for oscillatory covariance functions, partly because off-diagonal ranks grow with dimensionality.

Abstract

from arXiv · show

A number of problems in probability and statistics can be addressed using the multivariate normal (Gaussian) distribution. In the one-dimensional case, computing the probability for a given mean and variance simply requires the evaluation of the corresponding Gaussian density. In the $n$-dimensional setting, however, it requires the inversion of an $n \times n$ covariance matrix, $C$, as well as the evaluation of its determinant, $\det(C)$. In many cases, such as regression using Gaussian processes, the covariance matrix is of the form $C = σ^2 I + K$, where $K$ is computed using a specified covariance kernel which depends on the data and additional parameters (hyperparameters). The matrix $C$ is typically dense, causing standard direct methods for inversion and determinant evaluation to require $\mathcal O(n^3)$ work. This cost is prohibitive for large-scale modeling. Here, we show that for the most commonly used covariance functions, the matrix $C$ can be hierarchically factored into a product of block low-rank updates of the identity matrix, yielding an $\mathcal O (n\log^2 n) $ algorithm for inversion. More importantly, we show that this factorization enables the evaluation of the determinant $\det(C)$, permitting the direct calculation of probabilities in high dimensions under fairly broad assumptions on the kernel defining $K$. Our fast algorithm brings many problems in marginalization and the adaptation of hyperparameters within practical reach using a single CPU core. The combination of nearly optimal scaling in terms of problem size with high-performance computing resources will permit the modeling of previously intractable problems. We illustrate the performance of the scheme on standard covariance kernels.

I. INTRODUCTION

Gaussian-process inference becomes difficult at scale because dense covariance matrices make inversion and determinant evaluation O(n^3). The paper introduces hierarchical factorizations to reduce these costs for broad covariance kernels.

  • Dense n×n covariance matrices require O(n^3) work for both inversion and determinant evaluation, making high-dimensional Gaussian probabilities costly.The covariance matrix is typically symmetric, positive-definite, and dense.
  • The proposed HODLR factorization yields O(n log^2 n) inversion and O(n log n) determinant evaluation under suitable conditions.It factors the covariance matrix into block low-rank updates of the identity.
  • Gaussian processes are attractive for flexible continuous-function modeling, but computational complexity limits their practical use on large datasets.
  • Gaussian-process applications require covariance-matrix determinants for likelihood and evidence evaluations, including optimization and MCMC loops.
  • The methods target regression, classification, marginalization, and hyperparameter adaptation with large-scale Gaussian-process covariance matrices.

A. Prediction

Gaussian-process prediction computes a conditional distribution for a new input using covariance-based formulas. For large datasets, the central bottleneck is repeatedly applying or inverting the covariance matrix.

  • Prediction evaluates the conditional distribution of a new function value given observed inputs, outputs, and the new input.
  • The covariance vector between a new input and observed data reduces prediction and confidence estimation to matrix-vector multiplications.
  • Direct inversion of σ^2I + K scales as O(n^3), motivating the paper’s rapid covariance-matrix inversion algorithm.
  • The covariance function may depend on hyperparameters that represent physical quantities or support robustness and uncertainty quantification.
  • Marginalizing r nuisance parameters requires an r-dimensional numerical integral with roughly O(q^r) quadrature nodes.

C. Adaptation of hyperparameters

Hyperparameter fitting and marginalization require repeated covariance inversions and determinant evaluations. Existing accelerations often approximate the covariance or struggle to compute accurate determinants directly.

  • Hyperparameter inference uses nonlinear optimization because every covariance-matrix entry depends nonlinearly on the hyperparameters.
  • Each likelihood or evidence evaluation requires C^-1 and det(C), so repeated recalculation is central to optimization cost.
  • Low-rank covariance approximations can use Sherman-Morrison-Woodbury and Sylvester’s theorem in O(p^2n), but constructing general approximations typically costs at least O(pn^2).
  • State-space, sparsity, low-rank, and iterative approaches accelerate Gaussian-process calculations by introducing approximations at different mathematical levels.
  • Accurate determinant evaluation remains difficult because many accelerated methods still require O(n^3) work or rely on approximate estimators.

III. HIERARCHICAL MATRICES

HODLR matrices represent dense covariance matrices through recursively organized low-rank off-diagonal blocks. This structure supports hierarchical factorizations built from block-diagonal low-rank updates.

  • III. HIERARCHICAL MATRICES: Hierarchical matrix methods apply to covariance matrices and other dense operators that become data-sparse after suitable row and column ordering.
  • III. HIERARCHICAL MATRICES: A kd-tree ordering followed by a global reordering such as a Z-curve aligns data points with matrix rows and columns.
  • III. HIERARCHICAL MATRICES: HODLR matrices recursively represent off-diagonal blocks with low-rank factors across hierarchical scales.
  • III. HIERARCHICAL MATRICES: The HODLR representation can be manipulated into a product containing one block-diagonally dense factor and block-diagonal low-rank updates to the identity.
  • III. HIERARCHICAL MATRICES: The two-level factorization is an illustrative special case before the arbitrary-level HODLR construction is developed.

A. Fast low-rank approximation of off-diagonal blocks

The method rapidly approximates off-diagonal blocks and recursively factors HODLR matrices into block-diagonal low-rank updates, enabling efficient inversion.

  • Fast low-rank approximation of off-diagonal blocks: Fast low-rank factorizations of off-diagonal blocks are required to keep the overall HODLR method computationally tractable.The target cost for obtaining a rank-r factorization of an n × n matrix is at most O(r^2n), rather than direct SVD costs of O(mnr).
  • Fast low-rank approximation of off-diagonal blocks: Smooth covariance kernels support analytic approximations, including interpolation and multipole or Taylor expansions, for efficiently compressing matrix blocks.Polynomial interpolation can achieve near spectral accuracy when matrix entries come from a smooth function.
  • Fast low-rank approximation of off-diagonal blocks: When no a-priori matrix structure is available, linear-algebraic methods such as partial-pivoted LU, randomized algorithms, and rank-revealing factorizations provide alternatives.The paper describes partial-pivoted LU for off-diagonal blocks as constructing sums of rank-one matrices from columns and rows of the matrix.
  • Fast low-rank approximation of off-diagonal blocks: The partial-pivoted LU approach relies on sub-sampling for error estimation and may lose scalability or accuracy for highly oscillatory kernels or matrices with small-scale structure.For such cases, the paper identifies analytic methods as preferable when they provide more efficient, high-accuracy approximations.
  • HODLR matrix factorization: An n × n, κ-level HODLR matrix is recursively factored into κ + 1 block-diagonal matrices whose diagonal blocks are low-rank updates.The construction uses dense diagonal submatrices at the finest level and recursively reduces dense inversions to smaller blocks, applying Sherman-Morrison-Woodbury updates.
  • HODLR matrix factorization: O(n log n) inverse application follows because each factor is block diagonal or a block-diagonal low-rank update to the identity matrix.Recursive factorization across log n levels avoids repeatedly inverting large dense blocks, while Sherman-Morrison-Woodbury handles the low-rank updates.
  • HODLR matrix factorization: Nonuniformly distributed data points can slightly reduce factorization performance, while spatial tree structures may improve off-diagonal compression.The paper specifically mentions kd-tree-style constructions based on spatial considerations as a possible alternative to point-count partitioning.

C. Computational complexity

The HODLR factorization is constructed by approximating off-diagonal blocks at multiple levels and recursively factoring the matrix. Its factorization, inversion, and low-rank construction have near-linear complexity in n up to logarithmic factors.

  • Construction: HODLR construction first computes low-rank approximations for all off-diagonal blocks, then recursively factors the matrix into roughly O(log n) pieces.The factorization uses a multilevel hierarchy with user-defined finest-level block size pmax and precision ϵ.
  • Construction: O(rn log n) flops are required to construct all off-diagonal low-rank factorizations using cross approximation.Here r is the ϵ-rank of each sub-matrix.
  • Factorization: O(n log^2 n) work is required for the complete HODLR factorization when p and r are much smaller than n.The cost accumulates across the approximately log2 n hierarchy levels.
  • Inversion: O(n log n) time is required to apply the inverse after factorization, using O(log n) low-rank-update factors and the Sherman-Morrison-Woodbury formula.Each factor can be inverted in O(n) time because it is a low-rank update to the identity.

IV. DETERMINANT COMPUTATION

After HODLR factorization, determinant evaluation reduces to determinants of block-diagonal factors and their low-rank perturbations. Sylvester’s determinant theorem makes the overall computation O(n log n).

  • Determinant evaluation: Sylvester’s determinant theorem enables determinant evaluation at a cost of O(n log n) operations after HODLR factorization.This provides an efficient direct route to evaluating high-dimensional probabilities.
  • Determinant evaluation: The determinant of a block-diagonal factor is computed as the product of its diagonal-block determinants.The determinant is also multiplicative across the factors in the HODLR product.
  • Determinant evaluation: Each diagonal block is a low-rank perturbation of the identity, so Sylvester’s theorem computes its determinant efficiently.For a rank-p update to the identity, the determinant computation costs O(p^2n).
  • Complexity: Each factor determinant costs O(n), yielding total determinant complexity O(κn) with κ ∼ log n.The factors correspond to the levels of the HODLR hierarchy.

V. NUMERICAL RESULTS

The numerical study benchmarks HODLR assembly, factorization, inversion, and determinant computation for several covariance kernels and data dimensions. It also highlights the method’s applicability beyond cases where the covariance matrix itself is low rank.

  • Experimental setup: The experiments measure assembly, factorization, inversion, and determinant computation for n-dimensional datasets in one, two, and three dimensions.Gaussian and multiquadric covariance kernels are evaluated, with additional one-dimensional benchmarks for other kernels.
  • Kernel applicability: The hierarchical approach applies even when the covariance matrix is full rank, because its off-diagonal blocks can still have small ranks.The exponential covariance function is given as an example where K has full rank.
  • Gaussian-kernel scaling: Figure 5 compares covariance-factorization time for Gaussian kernels in one, two, and three dimensions against conventional direct calculation.The direct calculation is independent of dimension but scales as O(n^3).

A. Gaussian covariance

For Gaussian covariance kernels, the study evaluates points in one, two, and three dimensions and compares HODLR scaling with direct calculation.

  • Gaussian covariance: The Gaussian covariance kernel is evaluated at point pairs in one, two, and three dimensions.Results are aggregated in Table II and scaling is compared with direct calculation in Figure 5.
  • Gaussian covariance: Figure 5 compares HODLR scaling across embedded dimensions with conventional direct calculation.The comparison concerns the time required to factorize the covariance matrix.

B. Multiquadric covariance matrices

The paper reports timing results for multiquadric, exponential, inverse multiquadric, and biharmonic covariance matrices. Multiquadric scaling is virtually identical to the Gaussian case, while Figure 6 compares scaling across kernels.

  • Multiquadric covariance matrices: Multiquadric covariance functions are a frequently used class of radial basis functions.The paper presents analogous numerical results in one, two, and three dimensions.
  • Multiquadric covariance matrices: Multiquadric scaling is virtually identical to the Gaussian case, so the corresponding plot is omitted.The reported numerical results cover one, two, and three dimensions.
  • Multiquadric covariance matrices: Exponential covariance functions are evaluated in one dimension, with scaling compared across kernels in Figure 6.The comparison includes exponential, inverse multiquadric, and biharmonic covariance kernels.
  • Multiquadric covariance matrices: Timing results are reported in Tables V and VI for inverse multiquadric and biharmonic covariance functions.The inverse multiquadric and biharmonic kernels are frequently used in radial basis function interpolation and geostatistics.

E. Scaling in high dimensions

The high-dimensional experiments examine how the algorithm scales when data occupy unscaled or scaled hypercubes. Scaling saturates near d ≈ 10 for fixed sample size, while sparse high-dimensional data retain low off-diagonal ranks.

  • Unscaled data: For unscaled data in [−3, 3]^d, computational speed and accuracy increase as dimension d grows after some point.These results are reported in Table VII, although the paper calls this comparison unfair.
  • Scaled data: For scaled data in [−3/d, 3/d]^d, scaling saturates once d ≈ 10 because n is not increased with dimension.With fixed n, the points become sparse and the off-diagonal block ranks remain the same.
  • Curse of dimensionality: At high dimensions, 5000 points undersample spatial structure, producing off-diagonal ranks close to full and behavior equivalent to dense linear algebra.The paper identifies this as a manifestation of the curse of dimensionality.
  • Regression performance: The regression experiment generates 1024 synthetic data points and compares accelerated and exact regression across factorization precision ϵ.The experiment uses randomly selected x_j values in (−3, 3) and evaluates regression at the same points.
  • Regression performance: The absolute RMSE difference between dense and accelerated regression is proportional to factorization precision ϵ.Figure 8 reports this difference as a function of ϵ.

VI. CONCLUSIONS

The paper presents HODLR-based direct algorithms for covariance determinants, inverses, and matrix-vector products in Gaussian-process problems. The method avoids prior compromises, but performance degrades for dimensions above three and oscillatory covariance functions.

  • Conclusions: The HODLR-based algorithms compute covariance determinants, inverses, and matrix-vector products using fast, accurate, nearly optimal hierarchical direct linear algebra.The method targets covariance matrices encountered in Gaussian-process models and also applies to related classification and prediction problems.
  • Conclusions: The approach avoids sacrificing covariance-kernel fidelity, constructing global low-rank approximations, or paying dense full-rank computational costs.These were the compromises associated with previous acceleration attempts.
  • Limitations: The method’s performance degrades for covariance functions on R^d when d > 3 and for oscillatory covariance functions.The paper attributes part of the loss to larger off-diagonal ranks and part to compression difficulties.
  • Limitations: The authors identify extensions to these high-dimensional and oscillatory cases as current research.They also investigate high-dimensional anisotropic quadratures for marginalization and moment computation.
Loading 1403.6015v2…