I am testing the result. It turns out that for each time to calculate the loss and backpropagation, the repeating usage of some components will generate the error. If naive regulate retain_graph = True, the leaf tensor will have a problem in the next calculation. Therefore, my solution is to initialize the repeating parts each time. The resulting codes are in the followings for others might have the same problem:
for epoch in range(EPOCHS):
    sum_loss = 0
    target = torch.randn(26)
    
    for step in range(STEPS): 
        A0_B = torch.FloatTensor([[s_ for _ in range(8)] for _ in range(8)])
        qstep = torch.zeros(8, dtype=torch.float32)
        one_hot_location = F.one_hot(matrix_location,num_classes=8).view(-1).float()
        one_hot_location.requires_grad=True
        lamda_vector = Poisson(one_hot_location.float())
        for p in range(0, matrix_location.shape[0]):
            i_index = int(matrix_positions[p][0])
            j_index = int(matrix_positions[p][1])
            lamda = lamda_vector[p]*1.0
            epion = torch.rand(1).detach()
            A0_B[i_index][j_index] = lamda * (1 + 0.5 * epion)
        for index in range(len(A0_B)):
            transition_timing = (A0_B[index].clone() + external_inputs[step]) * 1.0
            # qval[step][index] = torch.matmul(sm(gm(sm(tempd.clone()))), tempd.clone())
            qstep[index] = torch.matmul(soft_max(gumbel_softmax(soft_max(transition_timing))), transition_timing)  
        qval[step] = qstep
        qlog = torch.log(qstep)
        # Calculate the backpropagation
        loss = criterion(qlog,torch.log(states[step]))
        sum_loss += loss.item()
        Distribution_Optimizer.zero_grad()
        loss.backward()  
#        print(one_hot_location.grad)
        Distribution_Optimizer.step()
#        print(((oldvals==Poisson.InputLayer.weight).sum().float())/124800)
#    writer.add_scalar('sum_loss', sum_loss,epoch)
    loss_list.append(sum_loss)
    print('Epoch: %.f, Loss: %.4f' %((epoch + 1),sum_loss))