I think that maybe I fixed the issue:
import torch
import torch.nn as nn
from torch.autograd import grad
r = torch.randn(3,10) #
d = torch.sum(r**2,dim=0, keepdim=True) + torch.sum(r**2,dim=0,keepdim=True).t() - 2 * r.t() @ r
d.requires_grad_(True)
y_true = torch.randn(10)

F = nn.Linear(10,10) #Dummy network
output = F(d)

E = torch.sum(output)
# E.backward(retain_graph=True)
y = -  grad(E, d, create_graph=True)[0].requires_grad_(True)
ysum = torch.sum(y,dim=0)

loss = torch.sum(ysum - y_true)

loss.backward() #Now I wish to do the standard loss backward propagation with an optimizer to update the network parameters
#optimizer.step()
print("Done")

This seems to run, I haven’t confirmed that it actually does what I intend but I think it does.