Multinomial logit Hessian matrix of the log-likelihood
params : array-like
hess : ndarray, (J*K, J*K)
where equals 1 if j = l and 0 otherwise.
The actual Hessian matrix has J**2 * K x K elements. Our Hessian is reshaped to be square (J*K, J*K) so that the solvers can use it.
This implementation does not take advantage of the symmetry of the Hessian and could probably be refactored for speed.