I’d like to compute the Hessian for array values and to get
some array with dimension (2,2) but instead I got this
julia> hessian(x -> sum(x.^3), [1 2; 3 4]) # uses linear indexing of x
4×4 Array{$Int,2}:
6 0 0 0
0 18 0 0
0 0 12 0
0 0 0 24
If I use the 2 x 10 array then I already have 20 x 20 matrix! Can I put the result in the other dimension?
Because my neural network will be very slow as I use hessian in the loss function
You have 20 inputs, so the Hessian is 20x20 by definition. If you want a 2x2 matrix, you are thinking of something other than a Hessian.
The result you got is correct. sum computes the sum of all elements, so your function x -> sum(x.^3) is indeed defined as \mathbb{R}^4\rightarrow\mathbb{R} and hence the hessian is a 4\times 4 matrix. If you feed to your function a 2\times 10 matrix, the function is \mathbb{R}^{20}\rightarrow\mathbb{R} and hence the hessian will be 20\times 20. In general if your function is \mathbb{R}^n\rightarrow\mathbb{R}, the hessian will be a n\times n matrix.
Yes, sorry.
I mean that I’d like to have the result with three dimensions (2,2,2) for 2x2 array and
the hessian with shape (2,2,10) for (2,10) array
The hessian is by definition a n\times n matrix. Are you sure you don’t mean something else? 
This makes me think you actually have a \mathbb{R}^2\rightarrow\mathbb{R} function and you want to evaluate the hessian at n distinct points (e.g if your array is (2,10) at 10 points) is this what you want to do?
A possible way I can think of is the following (there might be better ones)
julia> f(x) = sum(x.^3)
f (generic function with 1 methods)
julia> pnts = [[1,2], [3,4], [5,6], [7,8]]
4-element Vector{Vector{Int64}}:
[1, 2]
[3, 4]
[5, 6]
[7, 8]
julia> ForwardDiff.hessian.(f, pnts)
4-element Vector{Matrix{Int64}}:
[6 0; 0 12]
[18 0; 0 24]
[30 0; 0 36]
[42 0; 0 48]
Yes, that’s exactly what I wanted! Thanks