Most efficient way to calculate hessian of cost function in neural network

Question

I am coding a MLP network and I would like to implement the levenberg-marquardt algorithm. With levenberg-marquardt, the weights' update after each iteration is given by this formula:

W(t+1) = W(t) - (H(t)+ l(t)*I)^-1 * J

// W(t) is the matrix of weight at instant t
// H(t) is the `hessian` of the cost function
// l(t) is a learning rate
// J is the gradient of the cost function.

But I can't find an algorithm to calculate (or to have a acceptable estimation of) the hessian. How can I do that?