generalise partial gauss-newton but dont use it as default

bayesnewton/likelihoods.py

@@ -171,7 +171,7 @@ def gauss_newton(self, y, f):
         log_target = -0.5 * V.T @ V + self.log_normaliser(f)
         jacobian = -J.T @ V + ZJ
         # hessian_approx = -J.T @ J + self.log_normaliser_hessian(f)
-        hessian_approx = -J.T @ J - np.diag((ZJ**2).flatten())
+        hessian_approx = -J.T @ J - ZJ @ ZJ.T
         # second_order_term = -H.T * V
         return log_target, jacobian, hessian_approx  # , second_order_term
 
@@ -1241,7 +1241,7 @@ def conditional_moments(self, f):
         return conditional_expectation, conditional_covariance
 
 
-class HeteroscedasticNoise(MultiLatentLikelihood, PartialGaussNewtonMixin):
+class HeteroscedasticNoise(MultiLatentLikelihood, GaussNewtonMixin):
     """
     The Heteroscedastic Noise likelihood:
         p(y|f1,f2) = N(y|f1,link(f2)^2)