# \[ANN\] Vecchia.jl: linear-cost differentiable approximations for Gaussian log-likelihoods

**URL:** <https://discourse.julialang.org/t/ann-vecchia-jl-linear-cost-differentiable-approximations-for-gaussian-log-likelihoods/71915>\
**Category:** Package Announcements\
**Created:** [November 22, 2021, 8:56pm UTC](https://discourse.julialang.org/t/ann-vecchia-jl-linear-cost-differentiable-approximations-for-gaussian-log-likelihoods/71915 "2021-11-22T20:56:40Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [November 22, 2021, 8:56pm UTC](https://discourse.julialang.org/t/ann-vecchia-jl-linear-cost-differentiable-approximations-for-gaussian-log-likelihoods/71915/1 "2021-11-22T20:56:40Z")

</div>

Hi all. I’ve just registered [Vecchia.jl](https://github.com/cgeoga/Vecchia.jl), which provides functionality for a very convenient class of approximations for the Gaussian log-likelihood (as of now, only mean-zero, although for some technical reasons this is a bit less of a limitation than you might think). My preferred reference for those who like to read papers is [Stein/Chi/Welty 2004](https://rss.onlinelibrary.wiley.com/doi/abs/10.1046/j.1369-7412.2003.05512.x), but the idea can be described very simply. Any joint likelihood for a multivariate \textbf{z} can of course be decomposed into a product of conditional likelihoods as

\ell\_{\theta}(\textbf{z}) = \ell\_{\theta}(z\_1) \prod\_{j=2}^n \ell\_{\theta}(z\_j \; | \; \textbf{z}\_{1:(j-1)}).

The point of Vecchia approximations is that this quantity can in some circumstances be very well-approximated with

\ell\_{\theta}(\textbf{z}) \approx \ell\_{\theta}(z\_1) \prod\_{j=2}^n \ell\_{\theta}(z\_j \; | \; \textbf{z}\_{\sigma(j)}),

where \sigma(j) \subseteq [j-1] and is typically O(1) in size. In this setting and from this formulation, it is reasonably direct to see that you will evaluate O(n) many O(1)-sized conditional Gaussian likelihoods, giving a very attractive linear complexity in the data size. There are plenty of dials to turn, like taking “chunks” \textbf{z}\_j instead of scalar measurements z\_j, but for the sake of brevity I won’t discuss all that in the initial announcement post here.

As far as why this is sometimes a sensible thing to do, in the simplest case you could imagine a time series for \textbf{z}, and there is an obvious intuition that perhaps conditioning on the recent past (say, \sigma(j) = \{j-3, j-2, j-1\}) makes the the present and the far-past nearly independent. For those who would like a reference or more concrete discussion of when exactly that might be reasonably close to the truth, let me know and I can provide plenty of references or discussion here. But for the moment I won’t introduce any more technical or mathematical language unless somebody asks.

The package is pretty easy to use. The README gives a complete demo, but here is a slightly shortened demo:

```julia
using LinearAlgebra, StaticArrays, Vecchia

# VERY IMPORTANT FOR MULTITHREADING, since this is many small BLAS/LAPACK calls:
BLAS.set_num_threads(1)

# Covariance function, in this case Matern(v=3/2):
kfn(x,y,p) = p[1]*exp(-norm(x-y)/p[2])*(1.0+norm(x-y)/p[2])

# Locations for fake measurements, in this case 2048 of them, and fake data 
# (data NOT from the correction distribution, this is just a maximally simple demo):
const pts = [SVector{2, Float64}(randn(2)) for _ in 1:2048]
const dat = randn(length(pts))

# Create the VecchiaConfig:
const chunksize = 64
const num_conditioning_chunks = 3
const vecc = Vecchia.kdtreeconfig(dat, pts, chunksize, num_conditioning_chunks, kfn)

# Now you can evaluate the likelihood:
const sample_p = ones(2)
Vecchia.nll(vecc, sample_p)

```

Thanks to the beautiful Folds ecosystem of @tkf, this code thread-parallelizes very well. And as yet another bonus, it is fully AD-compatible, so you can using `ForwardDiff.hessian` and so on to perform higher-order optimization, which is a must for GP estimation in my opinion. The example files provide a full demonstration of directly using Ipopt to perform parameter estimation. I acknowledge that I’m an Ipopt cultist, but I really don’t think it gets any better than that.

Not demonstrated here but shown in the README, the package also provides functionality for the induced sparse approximation to the precision matrix. This depends even more on @tkf’s ecosystem and now also on @elrod’s `LoopVectorization`, but with all those fancy tools assembling the sparse approximation is only about 8x the time cost of the direct evaluation of the many small likelihoods. All things considered, I think that’s pretty cool.

Also, I’m releasing this on the MIT license instead of my preferred GPLv2 license that I’ve always used in the past, so if that was ever a deterrent to using code I’ve released, no need for that concern here. I’d really love to see more of the spatial statistics community adopt Julia, because I’m sick of writing R bindings so that R users can access all of the incredible functionality of Julia. So crossing my fingers that this helps…

Once the three day waiting period is up, you should be able to `]add Vecchia`. If you try it and have any issues or comments or questions, please do reach out. This is my first time registering a package and trying to deal with things like compat, so please also tell me if I’ve done something stupid there.

---

<div class="post-metadata">

**Author:** ![jkbest2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jkbest2/32/7350_2.png) [@jkbest2](https://discourse.julialang.org/u/jkbest2)\
**Post date:** [November 22, 2021, 10:06pm UTC](https://discourse.julialang.org/t/ann-vecchia-jl-linear-cost-differentiable-approximations-for-gaussian-log-likelihoods/71915/2 "2021-11-22T22:06:03Z")

</div>

This looks great, I’m excited to check it out! @ElOceanografo I think you’ll find this interesting as well.

---

<div class="post-metadata">

**Author:** ![ElOceanografo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eloceanografo/32/624_2.png) [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)\
**Post date:** [November 23, 2021, 1:52pm UTC](https://discourse.julialang.org/t/ann-vecchia-jl-linear-cost-differentiable-approximations-for-gaussian-log-likelihoods/71915/3 "2021-11-23T13:52:07Z")

</div>

Yes! This looks great @cgeoga, can’t wait to play around with it.

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [November 23, 2021, 2:16pm UTC](https://discourse.julialang.org/t/ann-vecchia-jl-linear-cost-differentiable-approximations-for-gaussian-log-likelihoods/71915/4 "2021-11-23T14:16:22Z")

</div>

Excellent! I’ll be very curious to hear from you and/or @jkbest2 if you do use it.

I should also clarify to anybody reading this thread that adding support for mean functions would be just a few lines of code. Since all the derivatives are AD-based it would be sufficient to evaluate for `data - meanfun.(locations)` instead of just `data`. I primarily haven’t added it because it’s not obvious what the “best” interface for it is. But if you or anybody else would be more interested if it had mean function stuff built in, I can easily plop that in and add that to the demo.

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [November 30, 2021, 4:03pm UTC](https://discourse.julialang.org/t/ann-vecchia-jl-linear-cost-differentiable-approximations-for-gaussian-log-likelihoods/71915/5 "2021-11-30T16:03:56Z")

</div>

Hey all–just an update that I’ve released a new version that now also offers the option of `LoopVectorization`-powered SIMD for the likelihood function. Particularly if you’re going to want sequential execution instead of threaded, the SIMD easily buys you a factor of two or three speedup. As of now it’s not completely stable for autodiff, but even with it just falls back to `@inbounds` that doesn’t come at a speed cost in the sequential case (unsure how exactly about many-core systems). But in some applications people just want a fast likelihood (right?), and so maybe even in its slightly incomplete form this would be of interest.

The usage is pretty simple once you’ve created the `VecchiaConfig` object like in the above code example or README:

```julia
# re-write the kernel function to accept scalar components of the two spatial
# coordinates.
function kfn_scalar(x1, x2, y1, y2, p)
  nrm = sqrt((x1-y1)^2 + (x2-y2)^2)
  p[1]*exp(-nrm/p[2])*(1+nrm/p[2])
end

# change some internals so that looping over points is LV-approved.
const vecc_s = Vecchia.scalarize(vecc, kfn_scalar) 

# Enjoy the AVX:
Vecchia.nll(vecc_s, sample_p)

```

Huge thanks to @elrod for help with the `@generated` functions that make this extension just as dimension-agnostic as the standard implementation that uses `SVectors` for coordinates.

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [June 3, 2022, 2:34am UTC](https://discourse.julialang.org/t/ann-vecchia-jl-linear-cost-differentiable-approximations-for-gaussian-log-likelihoods/71915/6 "2022-06-03T02:34:26Z")

</div>

Hey all–just another quick update. I’ve added a basic implementation of the “reverse” Cholesky factor for the sparse precision matrix. In a word, the typical Vecchia likelihood is written as a bunch of small likelihoods (see top post). But this can also be written in terms of a sparse precision matrix \Omega \, ``\approx" \Sigma^{-1} such that

\log \ell\_{\text{Vecchia}} (\theta) = \frac{1}{2} \left( -\log |\Omega(\theta)| + \textbf{z}^T \Omega(\theta) \textbf{z} \right)

for (suitably re-permuted!!) data \textbf{z}. The functionality to compute this \Omega has been in the package for a while as `Vecchia.precisionmatrix(vecc, params)` (see above for information about the `vecc` configuration object).

Further, though, there are some really nice ways to directly compute a symmetric factor of \Omega so that \Omega = U U^T _and U is upper triangular_. The code for creating U is now available in the package as `Vecchia.rchol(vecc, params)`, and is a decent bit faster than building \Omega. Not only is this pretty convenient for likelihoods, but it also makes things like simulations much easier. Probably the most important benefit of this, though, is that the U-based likelihood method does not require any sparse matrix factorizations (or any `SparseMatrixCSC` structures at all), and so the AD compatibility is excellent. The constructor for U does a lot of mutation, so I can’t promise that it will work for all AD systems. But with `ForwardDiff.jl` at least it’s great. Finally, I’m hoping that this helps in making this package composable with other things in the GP ecosystem, as you can use this function and just get a matrix factor out and then do whatever other crazy MCMC thing or whatever you want with it.

By default, `U=Vecchia.rchol(vecc, params)` gives you back a struct that has custom application methods (`U*x`, `U'*x`, and `U*U'*x` with some in-place variants) and a log-determinant. If you use big conditioning sets, these custom methods can be a lot faster than the sparse methods, I would guess because working with the dense chunks means getting to use BLAS. But if you want the actual sparse matrix, you can just call `sparse(U)` and you’re good to go.

As a quick reminder: this is all built on top of the beautiful `JuliaFolds` ecosystem. So you can really enjoy some good multi-threading speedup. **BUT** , remember to `BLAS.set_num_threads(1)` if you want to use multiple Julia threads (as in, `julia -t5` or whatever)!

Happy to provide more information and/or references if anybody is interested, but this post is getting long enough for now.

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [December 16, 2025, 8:12pm UTC](https://discourse.julialang.org/t/ann-vecchia-jl-linear-cost-differentiable-approximations-for-gaussian-log-likelihoods/71915/7 "2025-12-16T20:12:51Z")

</div>

# **Vecchia.jl version 0.11**

Hi all. I’ve just gone through and done a **big** rewrite of this package to try and take user-friendly design more seriously. I ended up removing a bit of functionality and paring down to try and focus on a smaller core of well-polished and performant routines. **If you are considering using this package or do use/try it, I would greatly appreciate your thoughts and feedback**.

Here are the two major use cases that I’ve tried to clean up and make as fun and easy to use as all of my favorite packages in the ecosystem are.

## GP fitting and interpolation

The setting here is that you have:

- locations `pts::Vector{SVector{D,Float64}}`
- measurements `data[j] = f(pts[j])`, and you think f is a (mean-zero) Gaussian process (and `data` may be a matrix, with each column representing an iid copy)
- a covariance function `kernel(pts[j], pts[k], params)`, which you think describes the covariance of `data[j]` and `data[k]`. Typically, `params` is unknown and needs to be estimated from data.

Then you specify your approximation with:

```julia
using Vecchia
approx = VecchiaApproximation(pts, kernel, data)

```

which has theoretically-motivated defaults that should perform well in a reasonably wide variety of settings. The docstrings and README go over options that you can play with, and I’ll continue to add new options (potentially via extensions to try and keep the base package light) (feedback and requests very welcome!).

You can then fit your model by loading in some extensions and bringing your own `NLPModels.jl`-complaint solver. For example, here is the estimation using the new `Uno` solver:

```julia
using NLPModels, ForwardDiff, UnoSolver
solver = NLPModelsSolver(uno; :preset=>"filtersqp")
mle = vecchia_estimate(approx, init_kernel_params, solver)

```

Finally, you can then predict at new locations `pred_pts` with

```julia
preds = predict(approx, pred_pts, mle)

```

## Fast preconditioners for kernel matrices

If you have some matrix `S` with `S[j,k] = kernel(pts[j], pts[k])` that is positive-definite, then you can make a \mathcal{O}(n)-cost preconditioner that is fast to assemble and (in at least some cases) works exceptionally well. Here is a demo of doing that using `cg` from the amazing `Krylov.jl`:

```julia
using StaticArrays, Vecchia, Krylov
using BesselK # for the matern kernel

# just an example of random locations
pts = rand(SVector{2,Float64}, 5_000)

# kernel matrix params:
params = [1.0, 0.1, 1.75]

# a dense matrix you don't want to factorize
S = [matern(x, y, params) for x in pts, y in pts]

# make the VecchiaApproximation, but no need to pass in data
# (because you don't need it). In my experience, it can be worthwhile
# to crank up the number of conditioning points beyond what works
# well for likelihoods though.
approx = VecchiaApproximation(pts, matern; ordering=KNNConditioning(30))
pre = rchol_preconditioner(approx, params)

# enjoy rapid convergence:
v = collect(1.0:length(pts)) # some RHS 
res = cg(S, v; M=pre, verbose=1) # ~15-20 iterations

```

I think this variety of preconditioner is a nice and complementary tool to other standard options like low-rank-plus-diagonal/Nystrom and stuff like that. Particularly if `kernel` is continuous but not smooth at the origin and you get slow spectral decay, this preconditioner may work much better than low rank-based alternatives.

## not implemented (yet): mean functions.

I have really been feeling some decision paralysis on the best way to implement mean functions. I basically always fit mean-zero GP models, so this isn’t something I’ve really had to deal with in my own work. But obviously this is something that needs to be implemented and offered for this package to be more broadly usable and feature-complete. If anybody has thoughts or suggestions, I would be _very_ interested and appreciative to hear them.

Thanks for reading!

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [April 2, 2026, 11:15pm UTC](https://discourse.julialang.org/t/ann-vecchia-jl-linear-cost-differentiable-approximations-for-gaussian-log-likelihoods/71915/8 "2026-04-02T23:15:08Z")

</div>

# Vecchia.jl version 0.13.[XXX]

Hi all. In continuing the big re-write and polishing, I’ve just released version 0.13 with the following new features:

- Parametric mean functions, including nonlinear ones (see example files for a demonstration)
- Carefully designed linear-cost prediction, including UQ and conditional simulation
- Expected Fisher information matrices as faster Hessian proxies for fitting mean-zero models (with the biggest selling point being that they require only the first derivative of the kernel function with respect to parameters), which is still designed with AD so that you can just plug in your favorite kernel and everything should just work.

Here is an example for fitting a mean-zero model with an expected Fisher matrix Hessian approximation and doing some prediction with confidence intervals (with is linear cost in the number of prediction points and can easily scale to millions+ of points).

```julia-auto
using NLPModels, ForwardDiff, UnoSolver
using sequentialknn_jll # an optional extension for blazing fast conditioning set design

# this script just simulates some data (pts_train, data_train, pts_pred,
# data_pred). In a real application, you would naturally bring your own
# locations and data.
include("example_setup.jl")

# specify the approximation for the data. This uses the Mat\'ern covariance
# function that looks like 
#
# matern(location_1, location_2, params)
#
# provided by BesselK.jl. But you can easily write your own kernel that takes
# two locations and a vector of parameters.
approx = VecchiaApproximation(pts_train, matern, data_train)

# compute the approximate mle that identifies the "best" kernel parameters.
solver = NLPModelsSolver(uno; preset="filtersqp", TR_radius=0.1)
mle = vecchia_estimate(approx, ones(3), solver;
                          expected_fisher=true, # a faster Hessian proxy
                          box_lower=[1e-8, 1e-8, 0.25], 
                          box_upper=[10.0, 10.0, 5.0])

# now predict at the un-observed locations.
preds = predict(approx, pts_pred, mle)
cmean = conditional_mean(preds)
cvars = conditional_variances(preds)

# summarize the first few results:
using Printf
@printf "\n\n ****A few predictions**** \n"
@printf "---------------------------\n"
@printf "True value Prediction\n"
@printf "---------------------------\n"
for j in 1:5
  @printf " % 01.3f % 01.3f ± %1.2f \n" data_pred[j] cmean[j] sqrt(cvars[j])*1.96
end

```

And the output of that print statement looks like this:

```julia-auto
****A few predictions****
---------------------------
True value Prediction
---------------------------
   2.382 2.389 ± 0.10 
  -7.371 -7.364 ± 0.06 
   2.000 2.056 ± 0.05 
  -3.833 -3.860 ± 0.08 
  -4.480 -4.721 ± 0.28 

```

Naturally those confidence intervals are so reasonable looking because I’m working with data simulated with the exact model that I’m fitting with, and so if you don’t take the probabilistic part of your model seriously then GPs and this whole framework becomes just one of many linear-cost function interpolators and the confidence interval functionality becomes much more dubious.

# Feedback

From anybody who tries this package would be great! It would be cool to try and tag a 1.0 release, but I would really appreciate hearing from at least a few people about the experience of using the code so that I don’t commit to patterns that are annoying.