class_docstr: This function takes in
  
      - the kernel interpolation term K_ZZ^{-1/2} k_ZX
      - the natural parameters of the variational distribution
  
  and returns
  
      - the predictive distribution mean/covariance
      - the inducing KL divergence KL( q(u) || p(u))
  
  However, the gradients will be with respect to the **cannonical parameters**
  of the variational distribution, rather than the **natural parameters**.
  This corresponds to performing natural gradient descent on the variational distribution.
code: def forward(
      ctx: FunctionCtx,
      interp_term: torch.Tensor,
      natural_vec: torch.Tensor,
      natural_mat: torch.Tensor,
  ) -> Tuple[torch.Tensor, torch.Tensor, torch.Tensor]:
      # Compute precision
      prec = natural_mat.mul(-2.0)
      diag = prec.diagonal(dim1=-1, dim2=-2).unsqueeze(-1)
  
      # Make sure that interp_term and natural_vec are the same batch shape
      batch_shape = torch.broadcast_shapes(
          interp_term.shape[:-2], natural_vec.shape[:-1]
      )
      expanded_interp_term = interp_term.expand(
          *batch_shape, *interp_term.shape[-2:]
      )
      expanded_natural_vec = natural_vec.expand(
          *batch_shape, natural_vec.size(-1)
      )
  
      # Compute necessary solves with the precision. We need
      # m = expec_vec = S * natural_vec
      # S K^{-1/2} k
      solves = linear_cg(
          prec.matmul,
          torch.cat(
              [
                  expanded_natural_vec.unsqueeze(-1),
                  expanded_interp_term,
              ],
              dim=-1,
          ),
          n_tridiag=0,
          max_iter=settings.max_cg_iterations.value(),
          tolerance=min(
              settings.eval_cg_tolerance.value(),
              settings.cg_tolerance.value(),
          ),
          max_tridiag_iter=settings.max_lanczos_quadrature_iterations.value(),
          preconditioner=lambda x: x / diag,
      )
      expec_vec = solves[..., 0]
      s_times_interp_term = solves[..., 1:]
  
      # Compute the interpolated mean
      # k^T K^{-1/2} m
      interp_mean = (
          s_times_interp_term.transpose(-1, -2)
          @ natural_vec.unsqueeze(-1)
      ).squeeze(-1)
  
      # Compute the interpolated variance
      # k^T K^{-1/2} S K^{-1/2} k = k^T K^{-1/2} (expec_mat - expec_vec expec_vec^T) K^{-1/2} k
      interp_var = (s_times_interp_term * interp_term).sum(dim=-2)
  
      # Let's not bother actually computing the KL-div in the foward pass
      # 1/2 ( -log | S | + tr(S) + m^T m - len(m) )
      # = 1/2 ( -log | expec_mat - expec_vec expec_vec^T | + tr(expec_mat) - len(m) )
      kl_div = torch.zeros_like(interp_mean[..., 0])
  
      # We're done!
      ctx.save_for_backward(
          interp_term,
          s_times_interp_term,
          interp_mean,
          natural_vec,
          expec_vec,
          prec,
      )
      return interp_mean, interp_var, kl_div
  
imports: from typing import Optional, Tuple
  import torch
  from linear_operator import to_linear_operator
  from linear_operator.operators import DiagLinearOperator, LinearOperator, MatmulLinearOperator, SumLinearOperator
  from linear_operator.utils import linear_cg
  from torch import Tensor
  from torch.autograd.function import FunctionCtx
  from ..None import settings
  from ..distributions import Delta, Distribution, MultivariateNormal
  from ..module import Module
  from ..utils.memoize import cached
  from ._variational_strategy import _VariationalStrategy
  from .natural_variational_distribution import NaturalVariationalDistribution
label_desc: Use image node_id 1 for calling the _NgdInterpTerms obj's underlying member method code with example usage: obj.forward(ctx, interp_term, natural_vec, natural_mat) and returns: interp_mean,  interp_var,  kl_div