How to design a Hessian Algorithm based on Gradient Algorithm?

Question

I have a question which I think should be pretty simple. I have written a subroutine in Fortran which will compute the Gradient at a point of the form

subroutine get_gradient(xyz, Grad)

where xyz is a 3x1 array and Grad is also a 3x1 array. The gradient method is a fourth order accurate method where the step is chosen inside the algorithm. Based on this, is there a simple way to use it to then turn around and generate the Hessian at xyz?

for example something like

for i = 1,3
  for j = 1,3
    txyz = xyz
    xyz(i) = xyz(i) + h
    xyz(j) = xyz(j) + h
    call get_gradient(xyz, Grad)
    Hess(i,:) = Hess(i,:) + Grad
    Hess(:,j) = Hess(:,j) + Grad
  end
end

I know the previous loop will not give the Hessian, but the idea is there - use the gradient method to construct the Hessian. Is there a simple algorithm which will do this?

Answer 1

Hessian is in simple terms the Jacobian of the gradient. So if you compute the Jacobian as in: http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant

form the gradient, you would obtain Hessian Matrix. While it may not be the best way to compute Hessian, it might be sufficient to do it.

Answer 2

As it turns out the answer is quite simple. In this example, f1, f2, f3, f4 are all 3x1 arrays.

  do ixyz = 1,3
    ! x+2h
    txyz = xyz
    txyz(ixyz) = txyz(ixyz) + two*h
    call get_grad(txyz, f1)
    ! x+h
    txyz = xyz
    txyz(ixyz) = txyz(ixyz) + h
    call get_grad(txyz, f2)
    ! x-h
    txyz = xyz
    txyz(ixyz) = txyz(ixyz) - h
    call get_grad(txyz, f3)
    ! x-2h
    txyz = xyz
    txyz(ixyz) = txyz(ixyz) - two*h
    call get_grad(txyz, f4)
    HessianRho(ixyz,:) = HessianRho(ixyz,:) + (-f1 + 8.0_rkind*(f2-f3) + f4)/(h*12.0_rkind)
    HessianRho(:,ixyz) = HessianRho(:,ixyz) + (-f1 + 8.0_rkind*(f2-f3) + f4)/(h*12.0_rkind)
  enddo
  HessianRho = HessianRho/two

This will generate the Hessian in a symmetric manner (which is why we do (ixyz,:) as well as (:,ixyz) and then divide by 2)