Avilay_Parekh:

torch.mm(C.transpose(0, 1), y1.view(-1, 1)).view(4, 4)


This needs to be
torch.mm(C.transpose(0, 1), torch.ones_like(y1).view(-1, 1)).view(4, 4)

When you want have a product that you want to backpropagate through, you replace the factor w.r.t. which you want to differentiate by the (appropriately expanded and summed up) output gradient, not by the output itself.
My favourite way of looking at derivatives of (multi)linear functions is in terms of Einstein summation notation, let’s take torch.nn.functional.bilinear as an example: It can be written for (inefficient) einsum as out = torch.einsum('bi,kij,bj->bk', left, weight, right). Now as torch.expand and torch.sum are dual to each other for taking derivatives and by the product differentiation rule, you have that
weight.grad = torch.einsum('bi,bk,bj->kij', left, out.grad, right) etc., so you swap in grad.out for one of the arguments and make the exchange between the right hand side and that factor in the equation, and you get the right gradient formula.
But now I got carried away, …
Best regards
Thomas