Hi,
apologies for being less clear.
I’d do something like the following (the probability of the log going wrong is not that large, given that the mean is 10 standard deviations from 0):
class PCEN(nn.Module):
    def __init__(self,in_features,bias=False):
        super(PCEN,self).__init__()

        self.log_alpha = nn.Parameter((torch.randn(in_features)*0.1+1.0).log_())
        self.log_delta = nn.Parameter((torch.randn(in_features)*0.1+1.0).log_())
        self.log_r     = nn.Parameter((torch.randn(in_features)*0.1+1.0).log_())
        self.eps = 0.00001

    def forward(self,x,smoother):
        alpha = self.log_alpha.exp().expand_as(x)
        delta = self.log_delta.exp().expand_as(x)
        r     = self.log_r.exp().expand_as(x)

        pcen = (x/(self.eps + smoother)**alpha + delta)**r - delta**r
        return pcen

#40 dimensional filterbank energies
pcen = PCEN(4)

#dummy data, energy is non-negative
feats = torch.randn(10,4).exp_()
smoother = torch.randn(10,4).exp_()

feats = Variable(feats)
smoother = Variable(smoother)

pcen_feats = pcen(feats,smoother)

This way, the backprop will just compute correct adjustments to the log parameters. My understanding is that they did it similarly.
I also took the liberty to generate the dummy data in pytorch directly and to make it positive with exp_. The fractional powers don’t really mix well with negative numbers (that is why your code got NaNs) and we all prefer positive energy. 
Best regards
Thomas