Thanks So much @ptrblck
I have written a minimal, executable code snippet to reproduce the issue.
Here is the code, you will see that the parameter is not updating. You can run it directly, without the need for additional libraries.
import numpy as np
import torch
import torch.nn as nn
from torch.autograd import Variable
class Torch_stat_model(nn.Module):
def __init__(self, initial_landmarks):  # landmark size (batch_size nx3)
    super(Torch_stat_model, self).__init__()
    self.initial_land = initial_landmarks

def forward(self, shape_params):  # shape-param size (batch_size nx1)
    land_numpy = (shape_params * self.initial_land).cpu().detach().numpy()
    return land_numpy

class Decoder_stat_model(nn.Module):
def __init__(self, initial_landmarks, scale=1.0):
    super(Decoder_stat_model, self).__init__()
    self.scale = scale
    self.model = Torch_stat_model(initial_landmarks)

def forward(self, theta):
    land_detached = self.model(theta)
    land_tensor = Variable(torch.from_numpy(land_detached).to(dtype=torch.float32).to(device=theta.device),
                           requires_grad=theta.requires_grad)
    return land_tensor * self.scale

device = torch.device(‘cuda’)
torch.manual_seed(0)
batch = 3
initial_landmark = torch.rand((batch, 48, 3)).to(device=device)
decoder = decoder = Decoder_stat_model(initial_landmark).to(device=device)
tparam = []
for i in range(batch):
alpha = np.random.normal(np.zeros(1), 2.0, 48)
tparam.append(alpha)
shape_parameters = torch.tensor(tparam, dtype=torch.float32,
device=device)
theta = torch.unsqueeze(shape_parameters, 2)
print(theta.shape)
target_landmarks = decoder(theta)
param = []
for i in range(batch):
alpha = np.random.normal(np.zeros(1), 1.5 + 1.0 / (i + 1), 48)
param.append(alpha)
shape_parameters_estimated = torch.zeros([batch, 48, 1], dtype=torch.float32, device=device, requires_grad=True)
loss = nn.MSELoss()
report_interval = 1
learning_rate = 1
optimizer = torch.optim.Adagrad([
shape_parameters_estimated],
lr=learning_rate)
max_iter = 5
for iteration in range(max_iter):
theta_estimated = shape_parameters_estimated
fit_landmarks = decoder(theta_estimated)
print(theta_estimated[0][0])

loss_lm = loss(target_landmarks, fit_landmarks)

loss_value = loss_lm

loss_value.backward()
optimizer.step()
optimizer.zero_grad()

# visualization
if iteration % report_interval == 0:
    print(f'epoch {iteration + 1}: loss={loss_value:.8f}')

You mentioned writing a custom `autograd.Function, how would I include that in the Decoder_stat_model class to attach the numpy array to the computation graph?