diff --git a/src/qibo/backends/abstract.py b/src/qibo/backends/abstract.py
index fe7d09195c..0d0918738a 100644
--- a/src/qibo/backends/abstract.py
+++ b/src/qibo/backends/abstract.py
@@ -358,7 +358,7 @@ def calculate_overlap(self, state1, state2):  # pragma: no cover
 
     @abc.abstractmethod
     def calculate_overlap_density_matrix(self, state1, state2):  # pragma: no cover
-        """Calculate norm of two density matrices."""
+        """Calculate overlap of two density matrices."""
         raise_error(NotImplementedError)
 
     @abc.abstractmethod
diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py
index dc48d71348..7dfbdb461d 100644
--- a/src/qibo/backends/numpy.py
+++ b/src/qibo/backends/numpy.py
@@ -687,7 +687,9 @@ def calculate_overlap(self, state1, state2):
         return self.np.abs(self.np.sum(self.np.conj(state1) * state2))
 
     def calculate_overlap_density_matrix(self, state1, state2):
-        raise_error(NotImplementedError)
+        state1 = self.cast(state1)
+        state2 = self.cast(state2)
+        return self.np.trace(self.np.transpose(self.np.conj(state1)) @ state2)
 
     def calculate_eigenvalues(self, matrix, k=6):
         if self.issparse(matrix):
diff --git a/src/qibo/callbacks.py b/src/qibo/callbacks.py
index f68e8d2312..2d96f795ba 100644
--- a/src/qibo/callbacks.py
+++ b/src/qibo/callbacks.py
@@ -206,7 +206,7 @@ def apply(self, backend, state):
         self.append(overlap)
         return overlap
 
-    def apply_density_matrix(self, backend, state):  # pragma: no cover
+    def apply_density_matrix(self, backend, state):
         overlap = backend.calculate_overlap_density_matrix(self.state, state)
         self.append(overlap)
         return overlap
diff --git a/tests/test_callbacks.py b/tests/test_callbacks.py
index d2e58b2cc8..16a5eabf6d 100644
--- a/tests/test_callbacks.py
+++ b/tests/test_callbacks.py
@@ -343,19 +343,29 @@ def test_norm(backend, density_matrix, seed):
     backend.assert_allclose(final_norm, target_norm)
 
 
+@pytest.mark.parametrize("seed", list(range(1, 5 + 1)))
 @pytest.mark.parametrize("density_matrix", [False, True])
-def test_overlap(backend, density_matrix):
-    state0 = np.random.random(4) + 1j * np.random.random(4)
-    state1 = np.random.random(4) + 1j * np.random.random(4)
+@pytest.mark.parametrize("nqubits", list(range(2, 6 + 1, 2)))
+def test_overlap(backend, nqubits, density_matrix, seed):
+    dims = 2**nqubits
+    if density_matrix:
+        state0 = random_density_matrix(dims, seed=seed, backend=backend)
+        state1 = random_density_matrix(dims, seed=seed + 1, backend=backend)
+    else:
+        state0 = random_statevector(dims, seed=seed, backend=backend)
+        state1 = random_statevector(dims, seed=seed + 1, backend=backend)
+
     overlap = callbacks.Overlap(state0)
-    overlap.nqubits = 2
+    overlap.nqubits = nqubits
+
     if density_matrix:
-        with pytest.raises(NotImplementedError):
-            overlap.apply_density_matrix(backend, state1)
+        final_overlap = overlap.apply_density_matrix(backend, state1)
+        target_overlap = np.trace(np.transpose(np.conj(state0)) @ state1)
     else:
         final_overlap = overlap.apply(backend, state1)
-        target_overlap = np.abs((state0.conj() * state1).sum())
-        backend.assert_allclose(final_overlap, target_overlap)
+        target_overlap = np.abs(np.sum(np.conj(state0) * state1))
+
+    backend.assert_allclose(final_overlap, target_overlap)
 
 
 @pytest.mark.parametrize("density_matrix", [False, True])