Skip to content

Commit

Permalink
perf(fock): Euler decomposition
Browse files Browse the repository at this point in the history
`FockSpace.get_linear_fock_operator` has been refactored to avoid embedding
into Fock space. For clarity, the Euler decomposition has been re-implemented
in `euler` in `decompositions`.
  • Loading branch information
Kolarovszki committed Feb 13, 2024
1 parent 9e98868 commit 270f32c
Show file tree
Hide file tree
Showing 6 changed files with 109 additions and 168 deletions.
43 changes: 32 additions & 11 deletions piquasso/_backends/fock/general/calculations.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,8 @@

from piquasso.instructions import gates

from piquasso._math.decompositions import euler

from piquasso.api.instruction import Instruction
from piquasso.api.result import Result

Expand All @@ -41,18 +43,22 @@ def passive_linear(
state._calculator, state._config
)

_apply_passive_linear(state, operator, instruction.modes)

return Result(state=state)


def _apply_passive_linear(state, interferometer, modes):
fock_operator = state._space.get_passive_fock_operator(
operator,
modes=instruction.modes,
interferometer,
modes=modes,
d=state._space.d,
)

state._density_matrix = (
fock_operator @ state._density_matrix @ fock_operator.conjugate().transpose()
)

return Result(state=state)


def particle_number_measurement(
state: FockState, instruction: Instruction, shots: int
Expand Down Expand Up @@ -290,16 +296,31 @@ def squeezing(state: FockState, instruction: Instruction, shots: int) -> Result:
def linear(
state: FockState, instruction: gates._ActiveLinearGate, shots: int
) -> Result:
operator = state._space.get_linear_fock_operator(
modes=instruction.modes,
passive_block=instruction._get_passive_block(state._calculator, state._config),
active_block=instruction._get_active_block(state._calculator, state._config),
)
calculator = state._calculator
modes = instruction.modes

state._density_matrix = (
operator @ state._density_matrix @ operator.conjugate().transpose()
np = calculator.np

passive_block = instruction._get_passive_block(state._calculator, state._config)
active_block = instruction._get_active_block(state._calculator, state._config)

symplectic = calculator.block(
[
[passive_block, active_block],
[np.conj(active_block), np.conj(passive_block)],
],
)

unitary_last, squeezings, unitary_first = euler(symplectic, calculator)

_apply_passive_linear(state, unitary_first, modes)

for mode, r in zip(instruction.modes, squeezings):
matrix = state._space.get_single_mode_squeezing_operator(r=r, phi=0.0)
_apply_active_gate_matrix_to_state(state, matrix, mode)

_apply_passive_linear(state, unitary_last, modes)

state.normalize()

return Result(state=state)
Expand Down
31 changes: 26 additions & 5 deletions piquasso/_backends/fock/pure/calculations/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,11 +20,15 @@
import random
import numpy as np

from .passive_linear import _apply_passive_linear

from ...calculations import calculate_state_index_matrix_list

from ..state import PureFockState
from ..batch_state import BatchPureFockState

from piquasso._math.decompositions import euler

from piquasso.instructions import gates

from piquasso.api.result import Result
Expand Down Expand Up @@ -318,13 +322,30 @@ def cubic_phase(state: PureFockState, instruction: Instruction, shots: int) -> R
def linear(
state: PureFockState, instruction: gates._ActiveLinearGate, shots: int
) -> Result:
operator = state._space.get_linear_fock_operator(
modes=instruction.modes,
passive_block=instruction._get_passive_block(state._calculator, state._config),
active_block=instruction._get_active_block(state._calculator, state._config),
calculator = state._calculator
modes = instruction.modes

np = calculator.np

passive_block = instruction._get_passive_block(state._calculator, state._config)
active_block = instruction._get_active_block(state._calculator, state._config)

symplectic = calculator.block(
[
[passive_block, active_block],
[np.conj(active_block), np.conj(passive_block)],
],
)

state._state_vector = operator @ state._state_vector
unitary_last, squeezings, unitary_first = euler(symplectic, calculator)

_apply_passive_linear(state, unitary_first, modes, calculator)

for mode, r in zip(instruction.modes, squeezings):
matrix = state._space.get_single_mode_squeezing_operator(r=r, phi=0.0)
_apply_active_gate_matrix_to_state(state, matrix, mode)

_apply_passive_linear(state, unitary_last, modes, calculator)

state.normalize()

Expand Down
12 changes: 7 additions & 5 deletions piquasso/_backends/fock/pure/calculations/passive_linear.py
Original file line number Diff line number Diff line change
Expand Up @@ -44,17 +44,19 @@ def passive_linear(
state._calculator, state._config
).astype(state._config.complex_dtype)

_apply_passive_linear(state, interferometer, instruction.modes, calculator)

return Result(state=state)


def _apply_passive_linear(state, interferometer, modes, calculator):
subspace = state._get_subspace(dim=len(interferometer))

subspace_transformations = _get_interferometer_on_fock_space(
interferometer, subspace, calculator
)

_apply_passive_gate_matrix_to_state(
state, subspace_transformations, instruction.modes
)

return Result(state=state)
_apply_passive_gate_matrix_to_state(state, subspace_transformations, modes)


def _get_interferometer_on_fock_space(interferometer, space, calculator):
Expand Down
41 changes: 38 additions & 3 deletions piquasso/_math/decompositions.py
Original file line number Diff line number Diff line change
Expand Up @@ -84,9 +84,11 @@ def _rotation_to_positive_above_diagonals(block_diagonal_matrix):

return block_diag(
*[
identity
if block_diagonal_matrix[2 * index, 2 * index + 1] > 0
else rotation
(
identity
if block_diagonal_matrix[2 * index, 2 * index + 1] > 0
else rotation
)
for index in range(d)
]
)
Expand Down Expand Up @@ -260,3 +262,36 @@ def mean_photon_number_gradient(scaling: float) -> float:
)

return result.root


def _complex_to_real(X):
Id = np.identity(len(X) // 2)

Check warning on line 268 in piquasso/_math/decompositions.py

View check run for this annotation

Codecov / codecov/patch

piquasso/_math/decompositions.py#L268

Added line #L268 was not covered by tests

W = 1 / np.sqrt(2) * np.block([[Id, 1j * Id], [Id, -1j * Id]])

Check warning on line 270 in piquasso/_math/decompositions.py

View check run for this annotation

Codecov / codecov/patch

piquasso/_math/decompositions.py#L270

Added line #L270 was not covered by tests

return W.conj().T @ X @ W

Check warning on line 272 in piquasso/_math/decompositions.py

View check run for this annotation

Codecov / codecov/patch

piquasso/_math/decompositions.py#L272

Added line #L272 was not covered by tests


def euler(symplectic, calculator):
np = calculator.np
d = len(symplectic) // 2

identity = np.identity(d)
zeros = np.zeros(shape=(d, d), dtype=complex)

U_orig, R = calculator.polar(symplectic, side="left")

K = calculator.block(
[
[identity, zeros],
[zeros, -identity],
],
)

H_active = 1j * K @ calculator.logm(R)

Z = 1j * H_active[:d, d:]

D, U = takagi(Z, calculator)

return U, D, np.conj(U).T @ U_orig[:d, :d]
143 changes: 0 additions & 143 deletions piquasso/_math/fock.py
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,6 @@
from piquasso.api.config import Config
from piquasso._math.indices import get_operator_index
from piquasso._math.combinatorics import partitions
from piquasso._math.decompositions import takagi
from piquasso._math.gradients import (
create_single_mode_displacement_gradient,
create_single_mode_squeezing_gradient,
Expand Down Expand Up @@ -240,148 +239,6 @@ def get_single_mode_cubic_phase_operator(
position = (annih.T + annih) * np.sqrt(hbar / 2)
return calculator.expm(1j * calculator.powm(position, 3) * (gamma / (3 * hbar)))

def embed_matrix(
self,
matrix: np.ndarray,
modes: Tuple[int, ...],
auxiliary_modes: Tuple[int, ...],
) -> np.ndarray:
indices = []
updates = []

for embedded_index, operator_basis in self.operator_basis_diagonal_on_modes(
modes=auxiliary_modes
):
index = (
operator_basis.ket.on_modes(modes=modes),
operator_basis.bra.on_modes(modes=modes),
)

indices.append(embedded_index)
updates.append(matrix[index][0])

return self._calculator.scatter(indices, updates, shape=(self.cardinality,) * 2)

def get_linear_fock_operator(
self,
*,
modes: Tuple[int, ...],
active_block: np.ndarray,
passive_block: np.ndarray,
) -> np.ndarray:
r"""The matrix of the symplectic transformation in Fock space.
Any symplectic transformation (in complex representation) can be written as
.. math::
S = \begin{bmatrix}
P & A \\
A^* & P^*
\end{bmatrix}.
As a first step, this symplectic matrix is polar decomposed:
.. math::
S = R U,
where :math:`R` is a hermitian matrix and :math:`U` is a unitary one. The polar
decomposition of a symplectic matrix is also a symplectic matrix, therefore
:math:`U` (being unitary) corresponds to a passive transformation, and :math:`R`
to an active one.
The symplectic matrix :math:`R` has the form
.. math::
\exp \left ( i K \begin{bmatrix}
0 & Z \\
Z^* & 0
\end{bmatrix}
\right ),
where :math:`Z` is a (complex) symmetric matrix. This can be decomposed via
Takagi decomposition:
.. math::
Z = U D U^T,
where :math:`U` is a unitary matrix and :math:`D` is a diagonal matrix. The
diagonal entries in :math:`D` correspond to squeezing amplitudes, and :math:`U`
corresponds to an interferometer.
Args:
modes (Tuple[int, ...]):
The modes on which the transformation should be applied.
active_block (np.ndarray): Active part of the symplectic transformation.
passive_block (np.ndarray): Passive part of the symplectic transformation.
Returns:
np.ndarray:
The resulting transformation, which could be applied to the state.
"""
np = self._calculator.np
fallback_np = self._calculator.fallback_np

d = len(modes)
identity = np.identity(d)
zeros = np.zeros_like(identity)

K = self._calculator.block(
[
[identity, zeros],
[zeros, -identity],
],
)

symplectic = self._calculator.block(
[
[passive_block, active_block],
[np.conj(active_block), np.conj(passive_block)],
],
)

U, R = self._calculator.polar(symplectic, side="left")

H_active = 1j * K @ self._calculator.logm(R)
H_passive = 1j * K @ self._calculator.logm(U)

singular_values, unitary = takagi(1j * H_active[:d, d:], self._calculator)

transformation = np.identity(self.cardinality, dtype=self.config.complex_dtype)

fock_operator = self.get_passive_fock_operator(
np.conj(unitary).T @ H_passive[:d, :d],
modes=modes,
d=self.d,
)

transformation = fock_operator @ transformation

for index, mode in enumerate(modes):
operator = self.get_single_mode_squeezing_operator(
r=singular_values[index],
phi=0.0,
)

squeezing_matrix = self.embed_matrix(
operator,
modes=(mode,),
auxiliary_modes=tuple(
fallback_np.delete(fallback_np.arange(self.d), (mode,))
),
)

transformation = squeezing_matrix @ transformation

fock_operator = self.get_passive_fock_operator(
unitary,
modes=modes,
d=self.d,
)

transformation = fock_operator @ transformation

return transformation

@property
def cardinality(self) -> int:
return cutoff_cardinality(cutoff=self.cutoff, d=self.d)
Expand Down
7 changes: 6 additions & 1 deletion tests/backends/test_backend_equivalence.py
Original file line number Diff line number Diff line change
Expand Up @@ -412,7 +412,12 @@ def test_fock_probabilities_with_two_mode_squeezing(SimulatorClass):

probabilities = state.fock_probabilities

assert all(probability >= 0 for probability in probabilities)
# NOTE: It should not happen that the probabilities are negative, but in some cases
# (possibly due to floating point errors) they turn out to be slightly negative.
assert all(
probability >= 0.0 or np.isclose(probability, 0.0)
for probability in probabilities
)
assert sum(probabilities) <= 1.0 or np.isclose(sum(probabilities), 1.0)

assert is_proportional(
Expand Down

0 comments on commit 270f32c

Please sign in to comment.