# AGP with linear predictors

**URL:** <https://discourse.julialang.org/t/agp-with-linear-predictors/59630>\
**Category:** Probabilistic Programming\
**Tags:** gaussian-process\
**Created:** [April 20, 2021, 12:43am UTC](https://discourse.julialang.org/t/agp-with-linear-predictors/59630 "2021-04-20T00:43:54Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![jzr](https://avatars.discourse-cdn.com/v4/letter/j/eb9ed0/32.png) [@jzr](https://discourse.julialang.org/u/jzr)\
**Post date:** [April 20, 2021, 12:43am UTC](https://discourse.julialang.org/t/agp-with-linear-predictors/59630/1 "2021-04-20T00:43:55Z")

</div>

Continuing my experiments from [another thread](https://discourse.julialang.org/t/faster-gp-fitting-on-million-observations/59519/11) using the suggestion of @theogf, I have some code using AGP. With a single predictor I can fit:

> ****
>
> ![image|25%](https://global.discourse-cdn.com/julialang/original/3X/a/b/ab6c8a5e8bf0b4381c0725042c646ff5e8d74a84.png)
> 
> ```julia
> using AugmentedGaussianProcesses
> const AGP = AugmentedGaussianProcesses
> using Distributions
> using CairoMakie
> using ThreadsX
> using DataFrames
> using StatsModels
> using GLM
> 
> N = 1000
> σ = 1.0
> df = let
> a = randn(N)
> b = randn(N)
> c = randn(N)
> ϵ = σ * randn(N)
> x = [a b c]
> t = 5rand(N)
> β = [.1, .2, .3]
> y = x * β .+ sin.(t) .+ ϵ
> sort(DataFrame(;a,b,c,t,y), :t)
> end
> 
> lm(@formula(y ~ a + b + c + sin(t)), df)
> 
> fig = Figure(resolution=(1000,2000))
> Axis(fig[1,1])
> scatter!(df.t, df.y, markersize=2)
> lines!(df.t, sin.(df.t))
> fig
> 
> ##
> 
> Ms = [4, 8, 16,];
> 
> kernel = SqExponentialKernel() + PeriodicKernel()
> 
> @time models = ThreadsX.map(Ms) do num_inducing
> m = SVGP(
> reshape(df.t, (nrow(df),1)),
> df.y, # First arguments are the input and output
> kernel, # Kernel
> GaussianLikelihood(σ), # Likelihood used
> AnalyticVI(), # Inference usede to solve the problem
> num_inducing; # Number of inducing points used
> optimiser=false, # Keep kernel parameters fixed
> Zoptimiser=false, # Keep inducing points locations fixed
> )
> @time train!(m, 100) # Train the model for 100 iterations
> m
> end
> 
> ##
> 
> @time for (i, (m, model)) in enumerate(zip(Ms, models))
> Axis(fig[i,1], title=string(m))
> n_grid = 40
> mins = minimum(df.t)
> maxs = maximum(df.t)
> x_grid = range(mins, maxs; length=n_grid)
> y_grid, sig_y_grid = proba_y(model, reshape(x_grid, :, 1))
> band!(x_grid, y_grid .- 2 * sqrt.(sig_y_grid), y_grid .+ 2 * sqrt.(sig_y_grid))
> scatter!(df.t, df.y, markersize=4, strokewidth=0, color=:black)
> scatter!.(vec(model.f[1].Z), Ref(zeros(dim(model.f[1]))); color=:orange, markersize=4, strokewidth=0)
> end
> 
> fig
> 
> ```

I would like to know if it’s possible to fit with multiple predictors. I don’t see docs about it in AGP, so I just guessed 🙂 . Currently I get

```julia
ERROR: LoadError: DimensionMismatch("Dimensions of the inputs are not matching : a = 2, b = 2, r = 1")
Stacktrace:
 [1] _evaluate(::KernelFunctions.Sinus{Float64}, ::SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true}, ::SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true})

```

> ****
>
> ```julia
> using AugmentedGaussianProcesses
> const AGP = AugmentedGaussianProcesses
> using Distributions
> using CairoMakie
> using ThreadsX
> using DataFrames
> using StatsModels
> using GLM
> 
> N = 1000
> σ = 1.0
> df = let
> a = randn(N)
> b = randn(N)
> c = randn(N)
> ϵ = σ * randn(N)
> x = [a b c]
> t = 5rand(N)
> β = [.1, .2, .3]
> y = x * β .+ sin.(t) .+ ϵ
> sort(DataFrame(;a,b,c,t,y), :t)
> end
> 
> lm(@formula(y ~ a + b + c + sin(t)), df)
> 
> fig = Figure(resolution=(1000,2000))
> Axis(fig[1,1])
> scatter!(df.t, df.y, markersize=2)
> lines!(df.t, sin.(df.t))
> fig
> 
> ##
> 
> Ms = [4, 8, 16,];
> 
> kernel = SqExponentialKernel() * PeriodicKernel()
> X = Matrix(df[:, [:t, :a]])
> @time models = map(Ms) do num_inducing
> m = SVGP(
> X,
> df.y, # First arguments are the input and output
> kernel, # Kernel
> GaussianLikelihood(σ), # Likelihood used
> AnalyticVI(), # Inference usede to solve the problem
> num_inducing; # Number of inducing points used
> optimiser=false, # Keep kernel parameters fixed
> Zoptimiser=false, # Keep inducing points locations fixed
> )
> @time train!(m, 100) # Train the model for 100 iterations
> m
> end
> 
> ##
> 
> @time for (i, (m, model)) in enumerate(zip(Ms, models))
> Axis(fig[i,1], title=string(m))
> n_grid = 40
> mins = minimum(df.t)
> maxs = maximum(df.t)
> x_grid = range(mins, maxs; length=n_grid)
> y_grid, sig_y_grid = proba_y(model, reshape(x_grid, :, 1))
> band!(x_grid, y_grid .- 2 * sqrt.(sig_y_grid), y_grid .+ 2 * sqrt.(sig_y_grid))
> scatter!(df.t, df.y, markersize=4, strokewidth=0, color=:black)
> scatter!.(vec(model.f[1].Z), Ref(zeros(dim(model.f[1]))); color=:orange, markersize=4, strokewidth=0)
> end
> 
> fig
> 
> ```

I think I want something like y = \beta\_1 a + \beta\_2 b + \beta\_3 c + \text{GP}(t). If I understand right, it’s possible to express the linear terms as GPs too, using linear kernels.

---

<div class="post-metadata">

**Author:** ![theogf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/theogf/32/1987_2.png) [@theogf](https://discourse.julialang.org/u/theogf)\
**Post date:** [April 20, 2021, 10:09am UTC](https://discourse.julialang.org/t/agp-with-linear-predictors/59630/2 "2021-04-20T10:09:34Z")

</div>

> `ERROR: LoadError: DimensionMismatch("Dimensions of the inputs are not matching : a = 2, b = 2, r = 1")`

This error message (which could be definitely improved), tells you that the `r` parameter you gave to your `PeriodicKernel` should be 2 dimensional to match the size of your inputs.

> [@jzr](#):
>
> I think I want something like y = \beta\_1 a + \beta\_2 b + \beta\_3 c + \text{GP}(t) . If I understand right, it’s possible to express the linear terms as GPs too, using linear kernels.

Yes that should definitely be possible. You can compose kernels by using the [`KernelTensorProduct`](https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable/kernels/#KernelFunctions.KernelTensorProduct) to use different kernels for each dimension.

---

<div class="post-metadata">

**Author:** ![jzr](https://avatars.discourse-cdn.com/v4/letter/j/eb9ed0/32.png) [@jzr](https://discourse.julialang.org/u/jzr)\
**Post date:** [April 20, 2021, 8:24pm UTC](https://discourse.julialang.org/t/agp-with-linear-predictors/59630/3 "2021-04-20T20:24:53Z")

</div>

Is this close?

```julia
ERROR: LoadError: MethodError: no method matching slices(::OptimIP{SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true},ColVecs{Float64,Array{Float64,2},SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true}},KmeansIP{SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true},ColVecs{Float64,Array{Float64,2},SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true}}},Nothing})

```

```julia
using AugmentedGaussianProcesses
using KernelFunctions

N = 1000
σ = 1.0

df = let
a = randn(N)
b = randn(N)
c = randn(N)
ϵ = σ * randn(N)
x = [a b c]
t = 5rand(N)
β = [.1, .2, .3]
y = x * β .+ sin.(t) .+ ϵ
(;a,b,c,t,y)
end

kernel = LinearKernel() ⊗ SqExponentialKernel()

X = [df.a df.t]

m = SVGP(
        X,
        df.y, # First arguments are the input and output
        kernel, # Kernel
        GaussianLikelihood(σ), # Likelihood used
        AnalyticVI(), # Inference usede to solve the problem
        4; # Number of inducing points used
        optimiser=false, # Keep kernel parameters fixed
        Zoptimiser=false, # Keep inducing points locations fixed
)
@time train!(m, 100) # Train the model for 100 iterations

```

---

<div class="post-metadata">

**Author:** ![theogf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/theogf/32/1987_2.png) [@theogf](https://discourse.julialang.org/u/theogf)\
**Post date:** [April 21, 2021, 4:50pm UTC](https://discourse.julialang.org/t/agp-with-linear-predictors/59630/4 "2021-04-21T16:50:14Z")

</div>

Oh I see! That’s a useful bug 🙂 It’s because of the implementation of InducingPoints.jl. I think that should be easily solvable.