Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

**BREAKING**: makes GP evaluation explicitly require RowVecs/ColVecs wrapper; calling GP with AbstractMatrix now errors #260

Closed
wants to merge 20 commits into from
Closed
Show file tree
Hide file tree
Changes from 13 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
@@ -1,7 +1,7 @@
name = "AbstractGPs"
uuid = "99985d1d-32ba-4be9-9821-2ec096f28918"
authors = ["JuliaGaussianProcesses Team"]
version = "0.5.4"
version = "0.6.0"

[deps]
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
Expand Down
2 changes: 1 addition & 1 deletion docs/Project.toml
Expand Up @@ -4,5 +4,5 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[compat]
AbstractGPs = "0.4, 0.5"
AbstractGPs = "0.4, 0.5, 0.6"
Documenter = "0.27"
10 changes: 10 additions & 0 deletions docs/src/api.md
Expand Up @@ -48,6 +48,16 @@ fx = f(x, Σ)
fx = AbstractGPs.FiniteGP(f, x, Σ)
```

!!! note

When considering higher-dimensional input features, these are often
provided as a matrix `X`. Evaluating `f(X)` interprets these as a vector of
vectors, where by default each _column_ of `X` is seen to represent an
input point. The implementation is handled by KernelFunctions.jl; see its
[documentation on Input
Types](https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable/api/#Vector-Valued-Inputs)
for more details.

st-- marked this conversation as resolved.
Show resolved Hide resolved
The `FiniteGP` has two API levels.
The [Primary Public API](@ref) should be supported by all `FiniteGP`s, while the [Secondary Public API](@ref) will only be supported by a subset.
Use only the primary API when possible.
Expand Down
77 changes: 40 additions & 37 deletions docs/src/concrete_features.md
@@ -1,26 +1,6 @@
# Features

## Plotting
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These changes seem unrelated, would have been cleaner to make them in a separate PR.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

true. just done that


### Plots.jl

We provide functions for plotting samples and predictions of Gaussian processes with [Plots.jl](https://github.com/JuliaPlots/Plots.jl). You can see some examples in the [One-dimensional regression](@ref) tutorial.

```@docs
AbstractGPs.RecipesBase.plot(::AbstractVector, ::AbstractGPs.FiniteGP)
AbstractGPs.RecipesBase.plot(::AbstractGPs.FiniteGP)
AbstractGPs.RecipesBase.plot(::AbstractVector, ::AbstractGPs.AbstractGP)
sampleplot
```

### Makie.jl

You can use the Julia package [AbstractGPsMakie.jl](https://github.com/JuliaGaussianProcesses/AbstractGPsMakie.jl) to plot Gaussian processes with [Makie.jl](https://github.com/JuliaPlots/Makie.jl).

![posterior animation](https://juliagaussianprocesses.github.io/AbstractGPsMakie.jl/stable/posterior_animation.mp4)


## Setup
### Setup

```julia
using AbstractGPs, Random
Expand All @@ -34,85 +14,108 @@ x = randn(rng, 10)
y = randn(rng, 10)
```

### Finite dimensional projection
## Finite dimensional projection
Look at the finite-dimensional projection of `f` at `x`, under zero-mean observation noise with variance `0.1`.
```julia
fx = f(x, 0.1)
```

### Sample from GP from the prior at `x` under noise.
## Sample from GP from the prior at `x` under noise.
```julia
y_sampled = rand(rng, fx)
```

### Compute the log marginal probability of `y`.
## Compute the log marginal probability of `y`.
```julia
logpdf(fx, y)
```

### Construct the posterior process implied by conditioning `f` at `x` on `y`.
## Construct the posterior process implied by conditioning `f` at `x` on `y`.
```julia
f_posterior = posterior(fx, y)
```

### A posterior process follows the `AbstractGP` interface, so the same functions which work on the posterior as on the prior.
## A posterior process follows the `AbstractGP` interface, so the same functions which work on the posterior as on the prior.

```julia
rand(rng, f_posterior(x))
logpdf(f_posterior(x), y)
```

### Compute the VFE approximation to the log marginal probability of `y`.
## Compute the VFE approximation to the log marginal probability of `y`.
Here, z is a set of pseudo-points.
```julia
z = randn(rng, 4)
```

#### Evidence Lower BOund (ELBO)
### Evidence Lower BOund (ELBO)
We provide a ready implentation of elbo w.r.t to the pseudo points. We can perform Variational Inference on pseudo-points by maximizing the ELBO term w.r.t pseudo-points `z` and any kernel parameters. For more information, see [examples](https://github.com/JuliaGaussianProcesses/AbstractGPs.jl/tree/master/examples).
```julia
elbo(VFE(f(z)), fx, y)
```

#### Construct the approximate posterior process implied by the VFE approximation.
### Construct the approximate posterior process implied by the VFE approximation.
The optimal pseudo-points obtained above can be used to create a approximate/sparse posterior. This can be used like a regular posterior in many cases.
```julia
f_approx_posterior = posterior(VFE(f(z)), fx, y)
```

#### An approximate posterior process is yet another `AbstractGP`, so you can do things with it like
### An approximate posterior process is yet another `AbstractGP`, so you can do things with it like
```julia
marginals(f_approx_posterior(x))
```

### Sequential Conditioning
## Sequential Conditioning
Sequential conditioning allows you to compute your posterior in an online fashion. We do this in an efficient manner by updating the cholesky factorisation of the covariance matrix and avoiding recomputing it from original covariance matrix.

```julia
# Define GP prior
f = GP(SqExponentialKernel())
```

#### Exact Posterior
### Exact Posterior

Generate posterior with the first batch of data by conditioning the prior on them:
```julia
# Generate posterior with the first batch of data on the prior f1.
p_fx = posterior(f(x[1:3], 0.1), y[1:3])
```

# Generate posterior with the second batch of data considering posterior p_fx1 as the prior.
Generate posterior with the second batch of data, considering the previous posterior `p_fx` as the prior:
```julia
p_p_fx = posterior(p_fx(x[4:10], 0.1), y[4:10])
```

#### Approximate Posterior
##### Adding observations in an sequential fashion
### Approximate Posterior
#### Adding observations in a sequential fashion
```julia
Z1 = rand(rng, 4)
Z2 = rand(rng, 3)
Z = vcat(Z1, Z2)
p_fx1 = posterior(VFE(f(Z)), f(x[1:7], 0.1), y[1:7])
u_p_fx1 = update_posterior(p_fx1, f(x[8:10], 0.1), y[8:10])
```
##### Adding pseudo-points in an sequential fashion

#### Adding pseudo-points in a sequential fashion
```julia
p_fx2 = posterior(VFE(f(Z1)), f(x, 0.1), y)
u_p_fx2 = update_posterior(p_fx2, f(Z2))
```

## Plotting

### Plots.jl

We provide functions for plotting samples and predictions of Gaussian processes with [Plots.jl](https://github.com/JuliaPlots/Plots.jl). You can see some examples in the [One-dimensional regression](@ref) tutorial.

```@docs
AbstractGPs.RecipesBase.plot(::AbstractVector, ::AbstractGPs.FiniteGP)
AbstractGPs.RecipesBase.plot(::AbstractGPs.FiniteGP)
AbstractGPs.RecipesBase.plot(::AbstractVector, ::AbstractGPs.AbstractGP)
sampleplot
```

### Makie.jl

You can use the Julia package [AbstractGPsMakie.jl](https://github.com/JuliaGaussianProcesses/AbstractGPsMakie.jl) to plot Gaussian processes with [Makie.jl](https://github.com/JuliaPlots/Makie.jl).

![posterior animation](https://juliagaussianprocesses.github.io/AbstractGPsMakie.jl/stable/posterior_animation.mp4)
2 changes: 1 addition & 1 deletion examples/intro-1d/Project.toml
Expand Up @@ -13,7 +13,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"

[compat]
AbstractGPs = "0.4, 0.5"
AbstractGPs = "0.4, 0.5, 0.6"
AdvancedHMC = "0.2, 0.3"
Distributions = "0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25"
DynamicHMC = "2.2, 3.1"
Expand Down
2 changes: 1 addition & 1 deletion examples/mauna-loa/Project.toml
Expand Up @@ -9,7 +9,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

[compat]
AbstractGPs = "0.5"
AbstractGPs = "0.5, 0.6"
CSV = "0.9"
DataFrames = "1"
Literate = "2"
Expand Down
2 changes: 2 additions & 0 deletions examples/test/Project.toml
@@ -0,0 +1,2 @@
[deps]
Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
17 changes: 13 additions & 4 deletions src/finite_gp_projection.jl
Expand Up @@ -20,10 +20,19 @@ function FiniteGP(f::AbstractGP, x::AbstractVector, σ²::Real=default_σ²)
return FiniteGP(f, x, Fill(σ², length(x)))
end

function FiniteGP(
f::AbstractGP, X::AbstractMatrix, σ²=default_σ²; obsdim::Int=KernelFunctions.defaultobs
)
return FiniteGP(f, KernelFunctions.vec_of_vecs(X; obsdim=obsdim), σ²)
function FiniteGP(f::AbstractGP, X::AbstractMatrix, σ²=default_σ²)
nrows, ncols = size(X)
throw(
ArgumentError(
"""projecting a GP on a finite number of points only accepts vectors (e.g. a vector of vectors for multi-dimensional features), but was called with a $nrows × $ncols matrix.
- If this is supposed to represent $ncols data points of dimension $nrows,
wrap it in ColVecs (e.g. `f(ColVecs(X), σ²)`).
- If this is supposed to represent $nrows data points of dimension $ncols,
wrap it in RowVecs (e.g. `f(RowVecs(X), σ²)`).
For more details, see https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable/api/#Input-Types
willtebbutt marked this conversation as resolved.
Show resolved Hide resolved
""",
),
)
end

## conversions
Expand Down
6 changes: 5 additions & 1 deletion src/sparse_approximations.jl
Expand Up @@ -281,7 +281,11 @@ function _compute_intermediates(fx::FiniteGP, y::AbstractVector{<:Real}, fz::Fin
end

function consistency_check(fx, y)
@assert length(fx) == length(y)
return length(fx) == length(y) || throw(
DimensionMismatch(
"length(fx) = $(length(fx)) and length(y) = $(length(y)) do not match"
),
)
end

function tr_Cf_invΣy(f::FiniteGP, Σy::Diagonal)
Expand Down