I want to accomplish a similar task and I would propose using a different method than Olden's method. I would be happy to learn the flaws with my proposed method. Let me mention a couple reasons not to use Olden's method. First, I think (but don't know) that it may be difficult to track the combinations of weights in a network with many layers. Second, the weights may not "accurately" reflect what is happening because of the activation functions.
I propose that a better way to accomplish this task is to add a small delta to each feature one at a time and use the predict command to observe the effect on the outcome variable. This is analogous to calculating the marginal effect of each variable (or the coefficients in a classical model). Note, that this could be time intensive on a large sample or with many features (though you could probably use a representative subsample for the prediction instead). Also note that this could be different from the true "importance" of the feature because e.g. the feature could have positive and negative effects that are important but combine to cancel. If that were a concern, then we would want to do something richer like dropping features from the model and observing the effect on the outcome. Here is my marginal effects code that accomplishes what I propose:
def marginal_effects(Data, model, delta):
inputcount = Data.shape[1]
baselevel = model.predict(Data).mean()
identity = np.identity(inputcount)
effects = np.zeros([inputcount, 1])
for i in range(inputcount):
addedterm = np.zeros([inputcount,inputcount])
addedterm[i,i] = delta
multiplier = identity + addedterm
deltaData = Data.dot(multiplier)
effects[i] = model.predict(deltaData).mean() - baselevel
return effects * (1/delta)
Then I use the function with a model (with a modest number of features) trained on the Data and a small delta like 0.01. The code is designed for a model with a single output, but could be easily generalized to multiple outputs.
Correction:
If you wanted to do Olden's Methdod (and I suppose it would be better if you had very many features) then it would probably just be a matter of matrix multiplying the weights from each layer in sequence. If you had 1,000 inputs and 100 neurons, then your weights should form a 1,000 x 100 matrix. If that 100 neuron layer was connected to a 100 neuron layer, then there should be a 100 x 100 matrix of weights and the two matrices could be multiplied to yield a 1,000 x 100 matrix. And so on, until the final matrix (assuming a single output) is 100 X 1 and multiplying it yields a 1,000 x 1 matrix with the "overall effect" of each input.