def h_update(hk, sk, yk):
  """Nocedal & Wright (1999) Eq. 9.2
     Pre-condition: hk must be positive definite"""
  yts = yk.dot(sk)
  assert yts > 0
  rho = 1 / yts
  v = -rho * yk.matrix_outer_product(sk)
  v.matrix_diagonal_add_in_place(1)
  hl = v.matrix_transpose_multiply(hk).matrix_multiply(v)
  hl += rho * sk.matrix_outer_product(sk)
  return hl