Gpytorch: Non-Square Hadamard Multiplication Fails [Bug]

[wjmaddox]

🐛 Bug

Looks like hadamard matrix multiplication is calling a root decomposition which is not well-defined for non-square lazy tensors. The fix is to check to see if the two matrices are non-square here.

To reproduce

* Code snippet to reproduce *

import torch
import gpytorch

a = gpytorch.lazify(torch.randn(30, 10))
b = gpytorch.lazify(torch.randn(10, 10))

ab = a.matmul(b)
# is not a matmul lazy tensor

(ab * ab) # fails

(ab.evaluate() * ab.evaluate()) # does not fail

* Stack trace/error message *

~/Documents/GitHub/gpytorch/gpytorch/lazy/lazy_tensor.py in _mul_matrix(self, other)
    508             left_lazy_tensor = self if self._root_decomposition_size() < other._root_decomposition_size() else other
    509             right_lazy_tensor = other if left_lazy_tensor is self else self
--> 510             return MulLazyTensor(left_lazy_tensor.root_decomposition(), right_lazy_tensor.root_decomposition())

RuntimeError: root_decomposition only operates on (batches of) square (symmetric) LazyTensors. Got a MatmulLazyTensor of size torch.Size([30, 10]).

Expected Behavior

Should return a Hadamard? or matmul lazy tensor rather than failing.

System information

Please complete the following information:

• gpytorch 1.2.0
• torch 1.7.0
• mac os

Additional context

I'd like to circumvent passing ab into a RootLazyTensor to get the diagonal out and so would like to directly do the Hadamard product here and reduce. (e.g. RootLazyTensor(ab).diag()).

[jacobrgardner]

@wjmaddox So, this is kind of intended. LazyTensor was always really intended to be PositiveDefiniteLazyTensor, and it's only recently with the addition of TriangularLazyTensor that we've really even considered matrices that arent spd.

This assumption is kind of littered throughout the package, and this is an example of it. Another example is that MatmulLazyTensor doesn't override _solve, and therefore will obviously return wrong results on non spd matrices. Realistically, MatmulLazyTensor(AB)._solve should do like B._solve(A._solve(...)) or something.

[wjmaddox]

Ahh, that would make sense. I'll close because of the easy workaround but separating out generic lazy tensors from PSD lazy tensors is probably something that would be useful in the future.