how to build a customer MultiOutput kernels ?

[linkzhao]

hi, i am new guy. in my project, I want to usea customer MultiOutput kernels, I tried, but failed. my kernel form refer to the attachment. can you help me code it by gpytorch when you are free? thank you very much!

[linkzhao]

I show my code as follow, please help me to check it ,thanks @gpleiss @arthus701

class CurlFreeKernel(gpytorch.kernels.Kernel):
def init(self, num_tasks, length_prior=None, length_constraint=None , **kwargs):

    super(CurlFreeKernel, self).__init__(**kwargs)

    # register the raw parameter
    self.register_parameter(
        name='raw_length', parameter=torch.nn.Parameter(torch.zeros(*self.batch_shape, 1, 1))
    )

    # set the parameter constraint to be positive, when nothing is specified
    if length_constraint is None:
        length_constraint = Positive()

    # register the constraint
    self.register_constraint("raw_length", length_constraint)

    self.num_tasks = num_tasks

# now set up the 'actual' paramter
@property
def length(self):
    # when accessing the parameter, apply the constraint transform
    return self.raw_length_constraint.transform(self.raw_length)

@length.setter
def length(self, value):
    return self._set_length(value)

def _set_length(self, value):
    if not torch.is_tensor(value):
        value = torch.as_tensor(value).to(self.raw_length)
    # when setting the paramater, transform the actual value to a raw one by applying the inverse transform
    self.initialize(raw_length=self.raw_length_constraint.inverse_transform(value))

def forward(self, x1, x2, diag=False, last_dim_is_batch=False, **params):

    x1_ = x1.div(self.length)
    x2_ = x2.div(self.length)

    I = torch.eye(3, 3, dtype=torch.float)

    num_x1=x1_.size(-2)
    num_x2=x2_.size(-2)
    # print("num_x1 ",num_x1)
    # print("num_x2 ",num_x2)
    with torch.no_grad():
        for i in range(num_x1):               
            for j in range(num_x2):
                diff=x1_[i]-x2_[j]                    
                dist_martri=( I-torch.matmul(diff, diff.T) )*(-0.5*torch.norm( diff,2)).exp()
                if j<1:
                    temp=dist_martri
                else:
                    temp = torch.cat((temp, dist_martri), 1)
            if i<1:
                result=temp
            else:
                result = torch.cat((result, temp), 0)
            #print("result:  ",result.shape)

    result= lazify(result)
    #print("result.numpy: ", result.numpy())
    return result

def num_outputs_per_input(self, x1, x2):
    """
    Given `n` data points `x1` and `m` datapoints `x2`, this multitask
    kernel returns an `(n*num_tasks) x (m*num_tasks)` covariance matrix.
    """
    return self.num_tasks

[linkzhao]

Should the matrix “diff ” be transposed

[arthus701]

No. The transpose here corresponds to an outer product, which is realized via broadcasting.

[linkzhao]

I tried an example and it was obvious that there were some problems via broadcasting:

[arthus701]

You're right. I updated the code, but didn't test it further. If you find the solution, I will update the kernel again to have the right code here...

[linkzhao]

So far, I have not found a solution. I will sudy the underlying implementation of gptorch. Maybe you can update the code? and i will help you test it

[gpleiss]

I think there may have been a few typos in @arthus701 's suggestion. Here's my attempt.
Again, this is un tested. I've commented what all the shapes should be.

A debugging tip: make sure the sizes of the object actually match the sizes in the comments.

from gpytorch.lazy import KroneckerLazyProduct

def forward(self, x1, x2, **params):
    # THIS IS JUST A CRUDE SKETCH! DON'T TAKE ANYTHING HERE FOR GRANTED!
    n = x1.size(-2)
    m = x2.size(-2)
    d = x1.size(-1)
    
    x1 = x1 / self.lengthscale
    x2 = x2 / self.lengthscale
    
    # Diff tensor: ... x n x m x d
    diff = x1.unsqueeze(-2) - x2.unsqueeze(-3)
    
    # Outer product term
    outer = diff.unsqueeze(-1) - diff.unsqueeze(-2)  # ... x n x m x d x d
    outer = outer.transpose(-2, -3).contiguous  # ... x n x d x m x d

    # Exp term
    norm = diff.pow(2).sum(dim=-1)  # ... x n x m
    exp_part = (-0.5 * norm).exp()
    exp_part = exp_part[..., :, None, :, None]  # ... x n x 1 x m x 1
    
    # Combine terms
    eye = torch.eye(d, dtype=x1.dtype, device=x1.device)  # ... x d x d
    eye = eye[..., None, :, None, :]  # ... x 1 x d x 1 x d
    kernel = (eye - outer) * exp_part
    
    kernel = kernel * self.outputscale / self.lengthscale.pow(2)
    return kernel

Since your kernel is multi-output, you will also need to add this method to the kernel:

    def num_outputs_per_input(self, x1, x2):
        return x1.size(-1)