diff --git a/src/Rings/mpoly.jl b/src/Rings/mpoly.jl
index bc6e8fbf7965..029c0788371a 100644
--- a/src/Rings/mpoly.jl
+++ b/src/Rings/mpoly.jl
@@ -1294,3 +1294,20 @@ end
 
 # assure compatibility with generic code for MPolyQuos:
 lift(f::MPolyRingElem) = f
+
+### Hessian matrix
+function hessian_matrix(f::MPolyElem)
+  R = parent(f)
+  n = nvars(R)
+  df = jacobi_matrix(f)
+  result = zero(MatrixSpace(R, n, n))
+  for i in 1:n
+    for j in i:n
+      result[i, j] = result[j, i] = derivative(df[i], j)
+    end
+  end
+  return result
+end
+
+hessian(f::MPolyElem) = det(hessian_matrix(f))
+
diff --git a/src/exports.jl b/src/exports.jl
index 37e8a6ab4d7c..52bc98d6290b 100644
--- a/src/exports.jl
+++ b/src/exports.jl
@@ -620,6 +620,8 @@ export helper_images
 export helper_kappa
 export helper_ring
 export hermitian_form
+export hessian
+export hessian_matrix
 export hilbert_basis
 export hilbert_function
 export hilbert_polynomial
diff --git a/test/Rings/mpoly.jl b/test/Rings/mpoly.jl
index 03b599b55d91..d85187593116 100644
--- a/test/Rings/mpoly.jl
+++ b/test/Rings/mpoly.jl
@@ -316,3 +316,11 @@ end
   equidimensional_decomposition_weak(I)
 end
 
+@testset "Hessian matrix" begin
+  R, (x, y, z) = QQ[:x, :y, :z]
+  f = x^2 + x*y^2 - z^3
+  H = hessian_matrix(f)
+  @test H == R[2 2*y 0; 2*y 2*x 0; 0 0 -6*z]
+  @test hessian(f) == -24*x*z + 24*y^2*z
+end
+