<p>Just to add to Mark's answer above, when using nat grads in non-conjugate models it can take a bit of tuning to get the best performance, and instability is potentially a problem. As Mark points out, the large steps that provide potentially faster convergence can also lead to the parameters ending up in in bad regions of the parameter space. When the variational approximation is good (i.e. the true and approximate posterior are close) then there is good reason to expect that the nat grad will perform well, but unfortunately there is no silver bullet in the general case. See <a href="https://arxiv.org/abs/1903.02984" rel="nofollow noreferrer">https://arxiv.org/abs/1903.02984</a> for some intuition.</p>
Exploding gradient for gpflow SVGP
<p>When optimizing a SVGP with Poisson Likelihood for a big data set I see what I think are exploding gradients.
After a few epochs I see a spiky drop of the ELBO, which then very slowly recovers after getting rid of all progress made before.
Roughly 21 iterations correspond to an Epoch.</p>

<p><img src="https://i.imgur.com/VbZMKqT.png" alt="ELBO"></p>

<p>This spike (at least the second one) resulted in a complete shift of the parameters (for vectors of parameters I just plotted the norm to see changes):
<img src="https://i.imgur.com/o7y4ySh.png" alt="Parameters"></p>

<p>How can I deal with that? My first approach would be to clip the gradient, but that seems to require digging around the gpflow code.</p>

<p><b>My Setup:</b></p>

<p>Training works via Natural Gradients for the variational parameters and ADAM for the rest, with a slowly (linearly) increasing schedule for the Natural Gradient Gamma.</p>

<p>The batch and inducing point sizes are as large as possible for my setup
(both 2^12, with the data set consisting of ~88k samples). I include 1e-5 jitter and initialize the inducing points with kmeans.</p>

<p>I use a combined kernel, consisting of a combination of RBF, Matern52, a periodic and a linear kernel on a total of 95 features (a lot of them due to a one-hot encoding), all learnable.
The lengthscales are transformed with gpflow.transforms.</p>

<pre><code>    with gpflow.defer_build():
        k1 = Matern52(input_dim=len(kernel_idxs["coords"]), active_dims=kernel_idxs["coords"], ARD=False)
        k2 = Periodic(input_dim=len(kernel_idxs["wday"]), active_dims=kernel_idxs["wday"])
        k3 = Linear(input_dim=len(kernel_idxs["onehot"]), active_dims=kernel_idxs["onehot"], ARD=True)
        k4 = RBF(input_dim=len(kernel_idxs["rest"]), active_dims=kernel_idxs["rest"], ARD=True)
        #
        k1.lengthscales.transform = gpflow.transforms.Exp()
        k2.lengthscales.transform = gpflow.transforms.Exp()
        k3.variance.transform = gpflow.transforms.Exp()
        k4.lengthscales.transform = gpflow.transforms.Exp()

        m = gpflow.models.SVGP(X, Y, k1 + k2 + k3 + k4, gpflow.likelihoods.Poisson(), Z,
                               mean_function=gpflow.mean_functions.Constant(c=np.ones(1)),
                               minibatch_size=MB_SIZE, name=NAME)
        m.mean_function.set_trainable(False)

    m.compile()
</code></pre>

<p><b>UPDATE: Using only ADAM</b>
Following the suggestion by Mark, I switched to ADAM only,
which helped me get rid of that sudden explosion. However, I still initialized with an epoch of natgrad only, which seems to save a lot of time.</p>

<p><img src="https://i.imgur.com/zS4wh9c.png" alt="Comparison of ADAM+NATGRAD(grey) vs ADAM only(red)"></p>

<p>In addition, the variational parameters seem to change a lot less abrupt (in terms of their norm at least). I guess they'll converge way slower now, but at least it's stable.
<img src="https://i.imgur.com/pMvv2LS.png" alt="Comparison of ADAM+NATGRAD(grey) vs ADAM only(red)"></p>