gp-shootout

Benchmark and compare several large-scale Gaussian process regression methods in 1D, 2D, and 3D, including our implementation of the equispaced Fourier method (EFGP) described in https://arxiv.org/abs/2210.10210 and https://arxiv.org/abs/2305.11065 . We focus on posterior mean prediction (kernel ridge regression). We also generate figures and tables for the first of the above papers.

Authors: Philip R Greengard, Alex H Barnett, Manas Rachh.

Installation

git clone --recurse-submodules https://github.com/flatironinstitute/gp-shootout.git

This will also download install some submodule packages (currently: RLCM, FLAM). In addition your system must also have the following required dependencies:

• MATLAB (tested on R2021a and b)

Dependencies specific to methods:

• For EFGP (and its option SSGP): FINUFFT (version 2.0 or later; please specify its location in startup.m).
• For SKI: Python 3.7, 3.8, or 3.9 (not 3.10, since MATLAB has not caught up), with python environment in which MATLAB was opened to include:
  • numpy
  • torch
  • gpytorch
  • pytorch
• For RLCM: C++ compiler (and RLCM installed as submodule)
• FLAM: a recent (2022) commit is needed (FLAM is installed as a submodule)

To test the basic installation, start MATLAB from the top-level gp-shootout directory (which will execute startup), then within MATLAB type test_all.

Advanced: to build then test all wrapped non-MATLAB methods:

1) make sure you can call python from matlab, eg via py.sys.version, then from MATLAB run test_all_nonmatlab.

Usage

If you did not start MATLAB from the top-level directory, then run startup to add required paths and apply useful settings. You may need to addpath to FINUFFT by hand if you forgot in startup.m.

Look in drivers for example scripts. You may try to run expt for a demo.

Generating figures and tables from the paper

All are run in MATLAB R2021a or R2021b unless stated. In order of appearance in paper:

• Fig. 1 {f:discr}: from paper repo equispaced_fourier_gps run discr_figs.jl in Julia. A couple of seconds runtime.
• Fig. 2 {f:materr}: drivers/fig_materr.m, takes around 1 minute.
• Fig. 3 {fig:cond1} and Fig. 4 {fig:cond_heatmap}: results/philip/cond_number/cond_number.m. Uses 1d data from efgp_tables. Takes ~5 seconds without exact 1e4 condition number calculation.
• Table 2 {table3}, ie the EFGP only N-scaling tests: paper_results/efgp_tables Data is subsampled from MAT files of size N=1e7 for each N. Time: about 1.5 hours for the full table, mostly due to reference solutions. Excluding Matern-3D examples, runtime is considerably faster (~10 minutes).
• Fig. 5 {acc_v_time}: paper_results/time_v_accuracy: input data is of size N=1e5 and is the same as the efgp tables of that size. Reference solutions can take a long time. For example, Matern 3d reference takes ~3 hours.
• Fig. 6 {fig:co2}: paper_results/c02: performs GP regression on CO2 data set described here -- https://www.tandfonline.com/doi/full/10.1080/01621459.2017.1419136. Time: ~15 seconds for l=5, 50, including reference solutions.
• Table 3 {tbig1} and Table 4 {tbig2}: paper_results/big_example. To generate the data, run gen_data.m, running this would require approximately 64 GB of memory. Set the data directory where you wish to store the data. To generate individual rows of the table, run big_example.m by setting appropriate values for iNvals, isig, itol, and alg. For table 3, iNvals = [1,2,3,4,5] and isig = 1, and for table 4, isig = iNvals. Options for alg are EFGP and FLAM. Set the directory to read the data from and the directory for storing the results. The memory requirements for storing all results is approximately 56 GB. Followed by this, run big_example_postproc.m, set values for iNvals and isig, and the directories where the data is stored and where the results are stored. Finally, the tex table for these results can be generated by running nicelatextable.m. For convenience, we have included the results generated for the paper in the results folder inside the big_example directory.

To do

• Octave compatibility (apart from Py engine?)
• understand empirical error breakdown as sigma->0 at large N
• understand CG num iters growing like 1/sqrt(tol) sometimes
• insert switch from EFGP to SSGP w/ same xi quad nodes
• explore EFGP accel via padding fftn to powers of only 2,3,5
• Move EFGP to a different repository and include it as a submodule

Done (CHANGELOG)

• rationalize interface for Philip dim-specific codes, move getL() into kernels
• add Matern nu=3/2, 5/2 in all 3 dims and to tester - a Matern kernel with nu as an input.
• CO2 dataset, 2D, N=1e6
• GPytorch (SKI) implementation
• Heaton et al. 2019 comparison datasets (two of), 2D, N=1e5
• posterior variances at target points for gp_naive
• testing routines for comparing methods to gp_naive for small problems
• rationalized dims 1,2,3 EFGP, handle arb shifted & scaled x data
• covar output to naive_gp?
• datasets -> data/*
• other methods -> algs/* (Python via system calls from matlab)
• 1d option for dense linear solve (vs. iterative)
• RLCM wrapped (SE ker only) via binary tmp file IO, tested in d=1,2,3.
• switched matlab pyenv engine setting to OutOfProcess to prevent crashes
• SKI fixed mem alloc problem
• protect from annoying FLAM chol fail if tol too large
• add Matern to RLCM (it's SE only so far)
• FLAMGP in 3D
• renamed m to mtot in efgp routines
• other top-level fig-generating driver scripts
• Update FLAMGP proxy setting to updated version proposed by Ken Ho