# Natural Gradients in Stochastic Variational Inference (SVI) for Gaussian Process Regression

Currently, I've hard times in understanding the natural gradients update in SVI method for Gaussian Process. I'm learning the SVI method for Gaussian Process through Gaussian Process for Big Data paper by James Hensman.

The paper said that Stochastic variational inference works by taking steps in the direction of the approximate natural gradient $$\tilde{g}(\theta) = G(\theta)^{-1}\frac{\partial L_3}{\partial \theta} = \frac{\partial L_3}{\partial \eta}$$. Therefore, the step in the direction of the natural gradient will be $$\theta_{(t+1)}=\theta_{(t)} + l \frac{\partial L_3}{\partial\eta}$$ where $$\theta$$ is canonical parameters and $$\eta$$ is expectation parameters of variational distribution $$q(u) = N(m,S)$$. Hence, if we substitute each of the natural gradients and canonical parameters, we will get: $$\theta_{2(t+1)} = -\frac{1}{2}S^{-1}_{(t+1)}=-\frac{1}{2}S^{-1}_{(t)} +l(-\frac{1}{2}\Lambda + \frac{1}{2}S^{-1}_{(t)})$$ and $$\theta_{1(t+1)} = S^{-1}_{(t+1)}m_{(t+1)} = S^{-1}_{(t)}m_{(t)} + l (\beta K^{-1}_{mm} K_{mn} y - S^{-1}_{(t)}m_{(t)})$$

Please note that, $$\Lambda = K^{-1}_{mm} + \Sigma_{i}\Lambda_i$$. As we can see that if we take $$l$$ which is the length step equals to 1 then we will get the same solution of variational distribution in Titsias' Sparse GP [2009] paper and it's the true posterior.

What I'm really confused about is how do we apply that optimization? Because we can just use length step equals to 1 and recover the same solution of Titsias' Sparse GP and automatically reach the true posterior. And also what makes this SVI cheap in computational cost? because if take step length equals to 1 we will have the same equation like in Titsias [2009].

Will really appreciate if someone can also give a simple example code of this SVI GP natural gradient update (I can read R and Python code). Thanks.

Here's the example of my code in R

library(dplyr)
library(tidyr)

set.seed(2)
train_pred <- seq(-5,5,length.out = 100)
e1 <- rnorm(length(train_pred), mean = 0, sd = 0.05)
train_resp <- sin(train_pred) + e1
df_sp <- data.frame(x = train_pred, y = train_resp)

kernel <- function(x1, x2, variance, lengthscale){
kern <- matrix(rep(0, times = length(x1)*length(x2)), ncol = length(x2), nrow = length(x1))
for(i in 1:length(x1)){
for(j in 1:length(x2)){
kern[i,j] = variance*exp(-((x1[i]-x2[j])^2)/(2*(lengthscale^2)))
}
}
return(kern)
}

induce_idx <- c(2,35,65,90,95)

induced_sample <- df_sp[induce_idx,]

minibatch <- function(df, batchsize = 2){
dummy_df <- df
dummy_df$$idx <- 1:nrow(dummy_df) dummy_id <- dummy_df$$idx
list_batch <- list()
i <- 0
while(length(dummy_id) > 0){
i <- i + 1
dummy <- sample(dummy_id, size = batchsize)
list_batch[[i]] <- df_sp %>% filter(id %in% dummy)
dummy_id <- dummy_id[!(dummy_id %in% dummy)]
}
return(list_batch)
}

list_mini <- minibatch(df_sp)
k_mm <- kernel(sample_df$$x1, sample_df$$x1, variance = 1, lengthscale = 1)
k_mn <- kernel(sample_df$$x1, df_sp$$x1, variance = 1, lengthscale = 1)
mean_svi <- list(matrix(0, nrow = nrow(k_mm), ncol = 1))
inverse_s_svi <- list(matrix(0, nrow = nrow(k_mm), ncol = nrow(k_mm)))
batchsize <- 2
steplength <- 0.02
sigma <- 1
counter <- 1
history <- c()
#The Update of Parameters using Natural Gradients
for(i in 1:4){
for(i in 1:(nrow(df_sp)/2)){
counter <- counter + 1
batch_df <- list_mini[[i]]
k_mb <- kernel(sample_df$$x1, batch_df$$x1, variance = 1, lengthscale = 1)
Lambda <- solve(k_mm) + (solve(k_mm) %*% k_mb %*% t(k_mb) %*% solve(k_mm))
inverse_s_svi[[counter]] <- inverse_s_svi[[counter - 1]] + steplength * (Lambda - inverse_s_svi[[counter - 1]])
old_s_svi <- inverse_s_svi[[counter - 1]]
mean_svi[[counter]] <- solve(inverse_s_svi[[counter]] + diag(1E-07, nrow = nrow(inverse_s_svi[[counter]]))) %*%
(old_s_svi %*% mean_svi[[counter - 1]] + steplength *
((sigma^(-2) * solve(k_mm) %*% k_mn %*% df_sp\$y1) - (old_s_svi %*% mean_svi[[counter - 1]])))
history <- c(history, sum((mean_svi[[counter]] - mean_svi[[counter - 1]])^2))
}
}
df_dummy <- data.frame(h = history, id = 1:length(history))
ggplot(df_dummy, aes(x = id, y = history)) + geom_line()


The parameter's turn out converge with small number of epoch. However, it doesn't converge into the correct value (compared to Sparse GP).