def METHOD_NAME(self, x1: Tensor, x2: Tensor, diag: bool = False, last_dim_is_batch: bool = False, **kwargs) -> Tensor:
    if last_dim_is_batch:
        x1 = x1.transpose(-1, -2).unsqueeze(-1)
        x2 = x2.transpose(-1, -2).unsqueeze(-1)
    num_dims = x1.size(-1)
    if not hasattr(self, "randn_weights"):
        self._init_weights(num_dims, self.num_samples)
    x1_eq_x2 = torch.equal(x1, x2)
    z1 = self._featurize(x1, normalize=False)
    if not x1_eq_x2:
        z2 = self._featurize(x2, normalize=False)
    else:
        z2 = z1
    D = float(self.num_samples)
    if diag:
        return (z1 * z2).sum(-1) / D
    if x1_eq_x2:
        # Exploit low rank structure, if there are fewer features than data points
        if z1.size(-1) < z2.size(-2):
            return LowRankRootLinearOperator(z1 / math.sqrt(D))
        else:
            return RootLinearOperator(z1 / math.sqrt(D))
    else:
        return MatmulLinearOperator(z1 / D, z2.transpose(-1, -2))