Hessian matrix of ML model

The docs state that hessian(f, x) ... where x is a real number or an array, ... which I interpret as if it should be a scalar array.
Doing something like

function f2(params)
	loss = Flux.mse(params[:, 2:end] * x .+ params[:, 1:1], y)
end
initial_params = [model[1].b model[1].W]
Zygote.hessian(f, initial_params)

seems to work for me. I think it is something with mutating arrays that Zygote does not like, have not read into it much but here is some discussion around the topic.