So, I now am using (apparently) correct values for normalize but still the predictions are very off.
network = Network()
network.cuda()    

criterion = nn.MSELoss()
optimizer = optim.Adam(network.parameters(), lr=0.0001)

loss_min = np.inf
num_epochs = 1

start_time = time.time()
for epoch in range(1,num_epochs+1):
    
    loss_train = 0
    loss_test = 0
    running_loss = 0
    
    
    network.train()
    print('size of train loader is: ', len(train_loader))

    for step in range(1,len(train_loader)+1):

        
        batch = next(iter(train_loader))
        images, landmarks = batch['image'], batch['landmarks']
        #RuntimeError: Given groups=1, weight of size [64, 3, 7, 7], expected input[64, 600, 800, 3] to have 3 channels, but got 600 channels instead
        #using permute below to fix the above error
        images = images.permute(0,3,1,2)
        
        images = images.cuda()
    
        landmarks = landmarks.view(landmarks.size(0),-1).cuda() 
    
        norm = transforms.Normalize([0.3809, 0.3810, 0.3810], [0.1127, 0.1129, 0.1130]) 
        for image in images:
            image = image.float()
            ##image = to_tensor(image) #TypeError: pic should be PIL Image or ndarray. Got <class 'torch.Tensor'>
            image = norm(image)
        
        print(type(images))
        ##landmarks = torchvision.transforms.Normalize(landmarks) #Do I need to normalize the target?
        
        predictions = network(images)
        
        # clear all the gradients before calculating them
        optimizer.zero_grad()
        
        print('predictions are: ', predictions.float())
        print('landmarks are: ', landmarks.float())
        # find the loss for the current step
        loss_train_step = criterion(predictions.float(), landmarks.float())
        
        
        loss_train_step = loss_train_step.to(torch.float32)
        print("loss_train_step before backward: ", loss_train_step)
        
        # calculate the gradients
        loss_train_step.backward()
        
        # update the parameters
        optimizer.step()
        
        print("loss_train_step after backward: ", loss_train_step)

        
        loss_train += loss_train_step.item()
        
        print("loss_train: ", loss_train)
        running_loss = loss_train/step
        print('step: ', step)
        print('running loss: ', running_loss)
        
        print_overwrite(step, len(train_loader), running_loss, 'train')
        
    network.eval() 
    with torch.no_grad():
        
        for step in range(1,len(test_loader)+1):
            
            batch = next(iter(train_loader))
            images, landmarks = batch['image'], batch['landmarks']
            images = images.permute(0,3,1,2)
            images = images.cuda()
            landmarks = landmarks.view(landmarks.size(0),-1).cuda()
        
            predictions = network(images)

            # find the loss for the current step
            loss_test_step = criterion(predictions, landmarks)

            loss_test += loss_test_step.item()
            running_loss = loss_test/step

            print_overwrite(step, len(test_loader), running_loss, 'Validation')
    
    loss_train /= len(train_loader)
    loss_test /= len(test_loader)
    
    print('\n--------------------------------------------------')
    print('Epoch: {}  Train Loss: {:.4f} Valid Loss: {:.4f}'.format(epoch, loss_train, loss_test))
    print('--------------------------------------------------')
    
    if loss_test < loss_min:
        loss_min = loss_test
        torch.save(network.state_dict(), '../moth_landmarks.pth') 
        print("\nMinimum Valid Loss of {:.4f} at epoch {}/{}".format(loss_min, epoch, num_epochs))
        print('Model Saved\n')
     
print('Training Complete')
print("Total Elapsed Time : {} s".format(time.time()-start_time))

Here’s the run results for the first step:
size of train loader is:  90
<class 'torch.Tensor'>
predictions are:  tensor([[ 0.5567,  0.1901,  0.1691, -0.5159,  0.7102,  0.1237,  0.1823, -0.0213],
        [ 0.5388,  0.1358,  0.1788, -0.5070,  0.7536,  0.1243,  0.1586, -0.0193],
        [ 0.5594,  0.1783,  0.1923, -0.5533,  0.7156,  0.1126,  0.1507, -0.0277],
        [ 0.5570,  0.1751,  0.1592, -0.5033,  0.7256,  0.1115,  0.1789,  0.0019],
        [ 0.5629,  0.1907,  0.1947, -0.4955,  0.6791,  0.1034,  0.1385,  0.0105],
        [ 0.5977,  0.1855,  0.1474, -0.5349,  0.7066,  0.0978,  0.1709,  0.0122],
        [ 0.5796,  0.1818,  0.2099, -0.5082,  0.7344,  0.1145,  0.1653, -0.0095],
        [ 0.5708,  0.1742,  0.1952, -0.5718,  0.7394,  0.1266,  0.1485, -0.0274]],
       device='cuda:0', grad_fn=<AddmmBackward>)
landmarks are:  tensor([[496.1400, 238.9700, 684.3768, 325.7215, 307.3357, 262.1889, 469.2918,
         323.8934],
        [489.7400, 240.3700, 708.0000, 253.0000, 327.0000, 331.0000, 485.0000,
         331.0000],
        [491.0200, 243.2600, 700.0000, 285.0000, 349.0000, 301.0000, 406.9148,
         349.3146],
        [496.0200, 244.2800, 587.0000, 115.0000, 336.0000, 147.0000, 492.0000,
         331.0000],
        [492.0800, 243.5100, 565.4799, 160.8031, 272.0000, 245.0000, 462.0000,
         344.0000],
        [501.0800, 241.7700, 720.0000, 286.0000, 304.0000, 310.0000, 513.1892,
         286.2780],
        [496.5700, 246.5700, 699.0000, 300.0000, 384.0000, 338.0000, 504.0000,
         326.0000],
        [494.0100, 239.8300, 539.0000, 150.0000, 345.0000, 116.0000, 441.0000,
         345.0000]], device='cuda:0')
loss_train_step before backward:  tensor(162903.9375, device='cuda:0', grad_fn=<MseLossBackward>)
loss_train_step after backward:  tensor(162903.9375, device='cuda:0', grad_fn=<MseLossBackward>)
loss_train:  162903.9375
step:  1
running loss:  162903.9375

I am not sure where I should divide by 255 in my code. Also, seems landmarks are also should be fed to to_tensor (ToTensor) since they are shown as the raw value.