docs: def _build_natgrad_step_op(self, model, q_mu_param, q_sqrt_param, xi_transform):
          """
          Implements equation 10 from
  
          @inproceedings{salimbeni18,
              title={Natural Gradients in Practice: Non-Conjugate  Variational Inference in Gaussian Process Models},
              author={Salimbeni, Hugh and Eleftheriadis, Stefanos and Hensman, James},
              booktitle={AISTATS},
              year={2018}
  
          In addition, for convenience with the rest of GPflow, this code computes ∂L/∂η using
          the chain rule:
  
          ∂L/∂η = (∂[q_μ, q_sqrt] / ∂η)(∂L / ∂[q_μ, q_sqrt])
  
          In total there are three derivative calculations:
          natgrad L w.r.t ξ  = (∂ξ / ∂nat ) [ (∂[q_μ, q_sqrt] / ∂η)(∂L / ∂[q_μ, q_sqrt]) ]^T
  
          Note that if ξ = nat or [q_μ, q_sqrt] some of these calculations are the identity.
  
          """
          objective = model.objective
          q_mu, q_sqrt = q_mu_param.constrained_tensor, q_sqrt_param.constrained_tensor
  
          # the three parameterizations as functions of [q_mu, q_sqrt]
          etas = meanvarsqrt_to_expectation(q_mu, q_sqrt)
          nats = meanvarsqrt_to_natural(q_mu, q_sqrt)
          xis = xi_transform.meanvarsqrt_to_xi(q_mu, q_sqrt)
  
          # we need these to calculate the relevant gradients
          _meanvarsqrt = expectation_to_meanvarsqrt(*etas)
  
          ## three derivatives
          # 1) the oridinary gpflow gradient
          dL_d_mean, dL_d_varsqrt = tf.gradients(objective, [q_mu, q_sqrt])
  
          # 2) the chain rule to get ∂L/∂η, where eta are the expectation parameters
          dL_detas = tf.gradients(_meanvarsqrt, etas, grad_ys=[dL_d_mean, dL_d_varsqrt])
  
          # 3) the forward mode gradient to calculate (∂ξ / ∂nat)(∂L / ∂η)^T,
          if isinstance(xi_transform, XiNat):
              nat_dL_xis = dL_detas
          else:
              _xis = xi_transform.naturals_to_xi(*nats)
              # this line should be removed if the placeholder_with_default problem is rectified
              _xis = [tf.reshape(_xis[0], q_mu_param.shape), tf.reshape(_xis[1], q_sqrt_param.shape)]
              nat_dL_xis = self._forward_gradients(_xis, nats, dL_detas)
  
          # perform natural gradient descent on the ξ parameters
          xis_new = [xi - self.gamma * nat_dL_xi for xi, nat_dL_xi in zip(xis, nat_dL_xis)]
  
          # transform back to the model parameters [q_μ, q_sqrt]
          mean_new, varsqrt_new = xi_transform.xi_to_meanvarsqrt(*xis_new)
  
          # these are the tensorflow variables to assign
          q_mu_u = q_mu_param.unconstrained_tensor
          q_sqrt_u = q_sqrt_param.unconstrained_tensor
  
          # so the transform to work for LowerTriangular
          mean_new.set_shape(q_mu_param.shape)
          varsqrt_new.set_shape(q_sqrt_param.shape)
  
          q_mu_assign = tf.assign(q_mu_u, q_mu_param.transform.backward_tensor(mean_new))
          q_sqrt_assign = tf.assign(q_sqrt_u, q_sqrt_param.transform.backward_tensor(varsqrt_new))
          return q_mu_assign, q_sqrt_assign