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.