If it’s not just a feature lacking temporarily, what alternatives are out there? I assume torch.solve works in most cases but that doesn’t leverage the positive definitiveness of a matrix:

l = torch.cholesky(x)
x_inv_u = torch.cholesky_solve(u, l)

Have you found any work around for this? I’m working with GPs and for numerical stability reasons the consensus is you need the cholesky_solve …which isn’t supported for .backward calls as you’ve pointed out.