Yes, you can basically just replace all numpy methods with the PyTorch ones:
device = 'cuda:0'
query_vecs = torch.randn(123, 32, device=device)
reference_vecs = torch.randn(456, 32, device=device)

dist = torch.cat([torch.sqrt(torch.sum((q - reference_vecs)**2, dim=1)) for q in query_vecs])