def METHOD_NAME(self, x1, x2):
    k_ux1 = to_dense(self.base_kernel(x1, self.inducing_points))
    if torch.equal(x1, x2):
        covar = LowRankRootLinearOperator(k_ux1.matmul(self._inducing_inv_root))
        # Diagonal correction for predictive posterior
        if not self.training and settings.sgpr_diagonal_correction.on():
            correction = (self.base_kernel(x1, x2, diag=True) - covar.diagonal(dim1=-1, dim2=-2)).clamp(0, math.inf)
            covar = LowRankRootAddedDiagLinearOperator(covar, DiagLinearOperator(correction))
    else:
        k_ux2 = to_dense(self.base_kernel(x2, self.inducing_points))
        covar = MatmulLinearOperator(
            k_ux1.matmul(self._inducing_inv_root), k_ux2.matmul(self._inducing_inv_root).transpose(-1, -2)
        )
    return covar