Soo some hints.
Variable was deprecated couple of years ago 
You don’t need to acess .data That just computes in-place modifications of the tensor.
You can define an optimizer and to pass tensors:
Some corrections would be
import torch
N=15
# load into pytorch
D = torch.rand(N, N).float()
H = torch.eye(N).float() - \
    1 / N * torch.ones(D.shape).float()
B = 0.5 * torch.matmul(torch.matmul(H, D ** 2), H)
pts = torch.rand(N, 3).float()

# declare constants
stop_loss = 1e-2
step_size = stop_loss / 3

# emulate simulate SVD
U = torch.rand(N, N, requires_grad=True)
s = torch.diag(torch.sort(torch.rand(N), descending=True).values).requires_grad_()
Vh = torch.rand(N, N, requires_grad=True)

optimizer = torch.optim.SGD([U,s,Vh], lr=step_size)
# find embedding
embed = torch.matmul(torch.sqrt(s), Vh)
X_hat = embed.T[:, :3]  # select the first 3 coordinates --> x, y, z

for i in range(100000):

    # calculate loss1: how close is our SVD function?
    delta1 = (torch.matmul(U.T, U) - torch.eye(N)) + \
             (torch.matmul(V.T, V) - torch.eye(N)) + \
             (torch.matmul(torch.matmul(U, s), Vh.T) - B)

    L1 = torch.norm(delta1, p=2)

    # calculate loss2: how close is our MDS approximation?
    delta2 = torch.nan_to_num(X_hat - pts, nan=0.0)

    L2 = torch.norm(delta2, p=2)

    # Backprop
    loss = L1 + L2
    loss.backward()

    # update
    optimizer.zero_grad()
    optimizer.step()

    if i % 1000 == 0:
        print('Loss is %s at iteration %i' % (loss, i))

    if abs(loss) < stop_loss:
        break

If you want to do your own optimization you should zeroe the gradients each step.
lastly, I don’t know whether this function is backpropagable or not torch.nan_to_num(X_hat - pts, nan=0.0)    But I cannot imagine why you would have NaNs
Just to answer your questions:
    U.data.data.zero_()
    s.data.data.zero_()
    V.data.data.zero_()

This is making zero the values of U,s and V each iteration. Thus the result will be static cos you are “resetting” the tensor. And yes, you can combine L=L1+L2