Thanks for your warm suggestion, @rasbt.
If I understand what you mean correctly, your method needs to first compute all layers’ representations (or intermediate outputs); and then uses those related outputs to obtain cost for each particular layer; finally, computes the gradients of each layer parameters and updates them respectively.
It almost solved my problem. However, as I mentioned in my problem, the hidden representations (outputs) don’t have any closed-form expressions, and meanwhile their values depend on spatially adjacent layers’ representations (outputs). Therefore, it seems like that those representations cannot be directly computed in a forward function.
In order to compute the representations (outputs) of all layers, it has to perform an iterative optimization process. So in this case, the problem becomes intractable.