FIX!(fermionic/gaussian): Passive linear fix

The time evolution with passive linear gates in the experimental fermionic
Gaussian simulator has been unintentionally implemented with time evolving
backwards. This issue got fixed in this patch.

piquasso/fermionic/gaussian/calculations.py

@@ -113,12 +113,12 @@ def passive_linear_gate(
     The transformation is equivalent to
 
     .. math::
-        \Gamma^{a^\dagger a} &\mapsto U \Gamma^{a^\dagger a} U^\dagger, \\\\
-        \Gamma^{a^\dagger a^\dagger} &\mapsto U \Gamma^{a^\dagger a^\dagger} U^T,
+        \Gamma^{a^\dagger a} &\mapsto \overline{U} \Gamma^{a^\dagger a} U^T, \\\\
+        \Gamma^{a^\dagger a^\dagger} &\mapsto \overline{U} \Gamma^{a^\dagger a^\dagger} U^\dagger,
 
     where :math:`U` is the unitary corresponding to the passive linear gate in the Dirac
     basis.
-    """
+    """  # noqa: E501
 
     unitary = instruction._get_passive_block(state._connector, state._config)
     modes = instruction.modes
@@ -135,13 +135,17 @@ def passive_linear_gate(
     select_rows = fallback_np.ix_(modes, all_modes)
 
     state._D = connector.assign(
-        state._D, select_columns, state._D[select_columns] @ unitary.conj().T
+        state._D, select_columns, state._D[select_columns] @ unitary.T
+    )
+    state._D = connector.assign(
+        state._D, select_rows, unitary.conj() @ state._D[modes, :]
     )
-    state._D = connector.assign(state._D, select_rows, unitary @ state._D[modes, :])
 
     state._E = connector.assign(
-        state._E, select_columns, state._E[select_columns] @ unitary.T
+        state._E, select_columns, state._E[select_columns] @ unitary.T.conj()
+    )
+    state._E = connector.assign(
+        state._E, select_rows, unitary.conj() @ state._E[select_rows]
     )
-    state._E = connector.assign(state._E, select_rows, unitary @ state._E[select_rows])
 
     return Result(state=state)

tests/fermionic/gaussian/test_gates.py

@@ -42,7 +42,7 @@ def test_GaussianHamiltonian_and_Interferometer_equivalence(
 
     unitary = expm(-2j * A.conj())
 
-    hamiltonian = np.block([[-A.conj(), zeros], [zeros, A]])
+    hamiltonian = np.block([[A, zeros], [zeros, -A.conj()]])
 
     simulator = pq.fermionic.GaussianSimulator(d=d, connector=connector)
 
@@ -115,12 +115,12 @@ def test_Interferometer_on_state_vector(connector):
 
     expected_covariance_matrix = np.array(
         [
-            [0.0, -0.932682, 0.156675, 0.30532, 0.106053, -0.033002],
-            [0.932682, 0.0, -0.30532, 0.156675, 0.033002, 0.106053],
-            [-0.156675, 0.30532, 0.0, 0.749426, 0.557812, 0.097148],
-            [-0.30532, -0.156675, -0.749426, 0.0, -0.097148, 0.557812],
-            [-0.106053, -0.033002, -0.557812, 0.097148, 0.0, -0.816744],
-            [0.033002, -0.106053, -0.097148, -0.557812, 0.816744, 0.0],
+            [0.0, -0.932682, -0.156675, 0.30532, -0.106053, -0.033002],
+            [0.932682, 0.0, -0.30532, -0.156675, 0.033002, -0.106053],
+            [0.156675, 0.30532, 0.0, 0.749426, -0.557812, 0.097148],
+            [-0.30532, 0.156675, -0.749426, 0.0, -0.097148, -0.557812],
+            [0.106053, -0.033002, 0.557812, 0.097148, 0.0, -0.816744],
+            [0.033002, 0.106053, -0.097148, 0.557812, 0.816744, 0.0],
         ]
     )
 
@@ -158,12 +158,12 @@ def test_Interferometer(connector):
 
     expected_covariance_matrix = np.array(
         [
-            [-0.0, 0.504253, 0.254254, 0.435585, 0.69889, 0.053882],
-            [-0.504253, 0.0, -0.444235, 0.643754, -0.259509, 0.258109],
-            [-0.254254, 0.444235, 0.0, -0.502168, -0.061072, 0.69434],
-            [-0.435585, -0.643754, 0.502168, -0.0, 0.260561, 0.275285],
-            [-0.69889, 0.259509, 0.061072, -0.260561, -0.0, -0.610399],
-            [-0.053882, -0.258109, -0.69434, -0.275285, 0.610399, 0.0],
+            [0.0, 0.469966, -0.19859, 0.780816, 0.228407, 0.279019],
+            [-0.469966, 0.0, -0.66337, 0.08203, -0.493825, -0.297456],
+            [0.19859, 0.66337, 0.0, -0.392254, -0.478024, 0.371661],
+            [-0.780816, -0.08203, 0.392254, 0.0, -0.069497, 0.474242],
+            [-0.228407, 0.493825, 0.478024, 0.069497, 0.0, -0.686026],
+            [-0.279019, 0.297456, -0.371661, -0.474242, 0.686026, 0.0],
         ]
     )
 
@@ -198,12 +198,12 @@ def test_Interferometer_on_subsystem(connector):
 
     expected_covariance_matrix = np.array(
         [
-            [0.0, 0.208302, 0.116836, 0.656998, 0.654506, -0.287986],
-            [-0.208302, -0.0, -0.968133, 0.023488, 0.129378, -0.045151],
-            [-0.116836, 0.968133, -0.0, -0.147157, -0.165207, -0.010927],
-            [-0.656998, -0.023488, 0.147157, -0.0, 0.281808, 0.683177],
-            [-0.654506, -0.129378, 0.165207, -0.281808, 0.0, -0.669458],
-            [0.287986, 0.045151, 0.010927, -0.683177, 0.669458, 0.0],
+            [0.0, 0.16447, 0.252142, 0.588858, 0.365978, -0.654737],
+            [-0.16447, 0.0, -0.951866, 0.162249, 0.145479, -0.139325],
+            [-0.252142, 0.951866, -0.0, -0.147157, -0.040228, 0.084272],
+            [-0.588858, -0.162249, 0.147157, -0.0, 0.672212, 0.391659],
+            [-0.365978, -0.145479, 0.040228, -0.672212, 0.0, -0.625626],
+            [0.654737, 0.139325, -0.084272, -0.391659, 0.625626, 0.0],
         ]
     )
 

tests/fermionic/gaussian/test_state.py

@@ -1458,23 +1458,136 @@ def test_density_matrix_Interferometer_2_by_2_simple(connector):
         ]
     )
 
-    with pq.Program() as program:
+    with pq.Program() as preparation:
         pq.Q(0, 1) | pq.StateVector([0, 1])
+
+    with pq.Program() as program:
+        pq.Q(0, 1) | preparation
         pq.Q(0, 1) | pq.Interferometer(U)
 
     simulator = pq.fermionic.GaussianSimulator(d=2, connector=connector)
-    state = simulator.execute(program).state
+    initial_state = simulator.execute(preparation).state
+    final_state = simulator.execute(program).state
 
-    initial_state_vector = np.array([0, 1, 0, 0])  # Lexicographic ordering
+    # The basis is ordered as [0, 0], [0, 1], [1, 0], [1, 1], but the unitary acting
+    # on the one-particle Hilbert space is understood in the basis [1, 0], [0, 1],
+    # therefore we have to flip it.
+    indices = [0, 2, 1, 3]
+
+    initial_state_density_matrix = initial_state.density_matrix[
+        np.ix_(indices, indices)
+    ]
+    final_state_density_matrix = final_state.density_matrix[np.ix_(indices, indices)]
+
+    fock_space_unitary = block_diag(np.array([1.0]), U, np.array([1.0]))
+
+    assert np.allclose(
+        final_state_density_matrix,
+        fock_space_unitary @ initial_state_density_matrix @ fock_space_unitary.T.conj(),
+    )
+
+
+@pytest.mark.monkey
+@for_all_connectors
+def test_density_matrix_Interferometer_2_by_2_random(
+    connector, generate_unitary_matrix
+):
+    U = generate_unitary_matrix(2)
+
+    with pq.Program() as preparation:
+        pq.Q(0, 1) | pq.StateVector([0, 1])
+
+    with pq.Program() as program:
+        pq.Q(0, 1) | preparation
+        pq.Q(0, 1) | pq.Interferometer(U)
+
+    simulator = pq.fermionic.GaussianSimulator(d=2, connector=connector)
+    initial_state = simulator.execute(preparation).state
+    final_state = simulator.execute(program).state
 
     # The basis is ordered as [0, 0], [0, 1], [1, 0], [1, 1], but the unitary acting
     # on the one-particle Hilbert space is understood in the basis [1, 0], [0, 1],
     # therefore we have to flip it.
-    U_big = block_diag(np.array([1.0]), U[::-1, ::-1], np.array([1.0]))
+    indices = [0, 2, 1, 3]
 
-    evolved_state_vector = U_big @ initial_state_vector
+    initial_state_density_matrix = initial_state.density_matrix[
+        np.ix_(indices, indices)
+    ]
+    final_state_density_matrix = final_state.density_matrix[np.ix_(indices, indices)]
+
+    fock_space_unitary = block_diag(np.array([1.0]), U, np.array([1.0]))
+
+    assert np.allclose(
+        final_state_density_matrix,
+        fock_space_unitary @ initial_state_density_matrix @ fock_space_unitary.T.conj(),
+    )
+
+
+@for_all_connectors
+def test_density_matrix_StateVector_ordering(connector):
+    d = 3
+
+    state_vectors = [
+        [0, 0, 0],
+        [0, 0, 1],
+        [0, 1, 0],
+        [0, 1, 1],
+        [1, 0, 0],
+        [1, 0, 1],
+        [1, 1, 0],
+        [1, 1, 1],
+    ]
+
+    simulator = pq.fermionic.GaussianSimulator(d=d, connector=connector)
+
+    for i in range(len(state_vectors)):
+        with pq.Program() as preparation:
+            pq.Q(0, 1, 2) | pq.StateVector(state_vectors[i])
+
+        state = simulator.execute(preparation).state
+
+        density_matrix = state.density_matrix
+
+        assert np.isclose(density_matrix[i, i], 1.0)
+
+
+@pytest.mark.monkey
+@for_all_connectors
+def test_density_matrix_Interferometer_3_by_3_random(
+    connector, generate_unitary_matrix
+):
+    d = 3
+
+    U = generate_unitary_matrix(d)
+
+    with pq.Program() as preparation:
+        pq.Q(0, 1, 2) | pq.StateVector([1, 0, 0])
+
+    with pq.Program() as program:
+        pq.Q(0, 1, 2) | preparation
+        pq.Q(0, 1, 2) | pq.Interferometer(U)
+
+    simulator = pq.fermionic.GaussianSimulator(d=d, connector=connector)
+    initial_state = simulator.execute(preparation).state
+    final_state = simulator.execute(program).state
+
+    # The basis is ordered as lexicographically, but the unitary acting
+    # on the one-particle Hilbert space is understood in a different basis,
+    # therefore we have to flip it.
+    indices = [0, 4, 2, 1, 6, 5, 3, 7]
+
+    initial_state_density_matrix = initial_state.density_matrix[
+        np.ix_(indices, indices)
+    ]
+    final_state_density_matrix = final_state.density_matrix[np.ix_(indices, indices)]
+
+    does_not_matter = np.random.rand(d, d)
+
+    fock_space_unitary = block_diag(
+        np.array([1.0]), U, does_not_matter, np.array([1.0])
+    )
 
     assert np.allclose(
-        state.density_matrix,
-        np.outer(evolved_state_vector.conj(), evolved_state_vector),
+        final_state_density_matrix,
+        fock_space_unitary @ initial_state_density_matrix @ fock_space_unitary.T.conj(),
     )