You could pre-compute it somewhere outside of the function like you are already doing it with tZ and B:
reduce[idx] = 0
...
idx = torch.eye(n).byte().cuda()

Alternatively you could multiply reduce with a mask.
As a side note, I’m not sure, if it’s a good idea to use a Python keyword as a variable name.
I’ve changed a bit and get a peak GPU utilization (GTX 1070) of approx. 62%.
def distMatrix(m):
    n = m.size(0)
    d = m.size(1)
    x = m.unsqueeze(1).expand(n, n, d)
    y = m.unsqueeze(0).expand(n, n, d) 
    return torch.sqrt(torch.pow(x - y, 2).sum(2) + 1e-4)

def lossIdx(tY):
    d = -distMatrix(tZ)+B
    sigmoidD = torch.sigmoid(d)
    
    #calculating cost
    r = torch.add(torch.mul(tY, torch.log(sigmoidD)), torch.mul((1-tY), torch.log(1-sigmoidD)))
    #remove diagonal
    r = r * idx
    return -r.sum()

device = 'cuda:0'

Y = np.asarray(np.random.randn(3000, 3000))

n = np.shape(Y)[0]
k = 2
idx = 1 - torch.eye(n).to(device)

Z = np.random.rand(n,k)

tZ = torch.tensor(Z, dtype=torch.float, requires_grad=True, device=device)
B = torch.tensor([0.], requires_grad=True, device=device)
tY = torch.from_numpy(Y.astype(np.float32))
tY.requires_grad_(True)
tY = tY.to(device)

optimizer = torch.optim.Adam([tZ, B], lr = 1e-3)

steps = 100
torch.cuda.synchronize()
t0 = time.time()
for i in range(steps):
    optimizer.zero_grad()
    l = lossIdx(tY)#.cuda()
    l.backward(retain_graph=False)
    optimizer.step()

torch.cuda.synchronize()
print(time.time() - t0)