diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py
index f280a05640..7248b2348e 100644
--- a/src/qibo/backends/numpy.py
+++ b/src/qibo/backends/numpy.py
@@ -400,6 +400,17 @@ def execute_circuit(self, circuit, initial_state=None, nshots=1000):
                 return self.execute_circuit(initial_state + circuit, None, nshots)
         elif initial_state is not None:
             initial_state = self.cast(initial_state)
+            valid_shape = (
+                2 * (2**circuit.nqubits,)
+                if circuit.density_matrix
+                else (2**circuit.nqubits,)
+            )
+            if tuple(initial_state.shape) != valid_shape:
+                raise_error(
+                    ValueError,
+                    f"Given initial state has shape {initial_state.shape} instead of "
+                    f"the expected {valid_shape}.",
+                )
 
         if circuit.repeated_execution:
             if circuit.measurements or circuit.has_collapse:
diff --git a/tests/test_models_circuit_execution.py b/tests/test_models_circuit_execution.py
index bf493bc893..0158cf1001 100644
--- a/tests/test_models_circuit_execution.py
+++ b/tests/test_models_circuit_execution.py
@@ -2,6 +2,7 @@
 import pytest
 
 from qibo import Circuit, gates
+from qibo.quantum_info import random_density_matrix
 
 
 def test_eager_execute(backend, accelerators):
@@ -91,8 +92,6 @@ def test_final_state_property(backend):
 
 
 def test_density_matrix_circuit(backend):
-    from qibo.quantum_info import random_density_matrix
-
     initial_rho = random_density_matrix(2**3, backend=backend)
 
     c = Circuit(3, density_matrix=True)
@@ -146,3 +145,14 @@ def test_initial_state_error(backend):
 
     with pytest.raises(ValueError):
         backend.execute_circuit(c, c1)
+
+
+@pytest.mark.parametrize("density_matrix", [False, True])
+def test_initial_state_shape_error(backend, density_matrix):
+    nqubits = 2
+    c = Circuit(nqubits, density_matrix=density_matrix)
+    c.add(gates.X(i) for i in range(nqubits))
+
+    initial_state = random_density_matrix(2, backend=backend)
+    with pytest.raises(ValueError):
+        backend.execute_circuit(c, initial_state=initial_state)