From 0b6909ddcc42e951ea35a1390ac55e9d85ae47a5 Mon Sep 17 00:00:00 2001 From: Kemal Bidzhiev Date: Wed, 25 Mar 2026 16:06:09 +0100 Subject: [PATCH 01/13] bug + docs --- emu_mps/hamiltonian.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/emu_mps/hamiltonian.py b/emu_mps/hamiltonian.py index e19e1723..9401f4a1 100644 --- a/emu_mps/hamiltonian.py +++ b/emu_mps/hamiltonian.py @@ -27,6 +27,12 @@ class Operators: class HamiltonianMPOFactors(ABC): + """Abstract generator for MPO factors of a two-body Hamiltonian. + + Subclasses implement the local MPO tensor at each position in the chain: + first site, left half, middle, right half, and last site. + """ + def __init__(self, interaction_matrix: torch.Tensor, dim: int = 2): assert interaction_matrix.ndim == 2, "interaction matrix is not a matrix" assert ( @@ -40,6 +46,7 @@ def __init__(self, interaction_matrix: torch.Tensor, dim: int = 2): self.identity = Operators.id if self.dim == 2 else Operators.id_3x3 def __iter__(self) -> Iterator[torch.Tensor]: + """Yield the full ordered list of MPO factors for the Hamiltonian.""" yield self.first_factor() for n in range(1, self.middle): @@ -55,22 +62,27 @@ def __iter__(self) -> Iterator[torch.Tensor]: @abstractmethod def first_factor(self) -> torch.Tensor: + """Return the MPO factor for the first site.""" pass @abstractmethod def left_factor(self, n: int) -> torch.Tensor: + """Return the MPO factor for site ``n`` in the left half of the chain.""" pass @abstractmethod def middle_factor(self) -> torch.Tensor: + """Return the MPO factor at the central site bridging both halves.""" pass @abstractmethod def right_factor(self, n: int) -> torch.Tensor: + """Return the MPO factor for site ``n`` in the right half of the chain.""" pass @abstractmethod def last_factor(self) -> torch.Tensor: + """Return the MPO factor for the last site.""" pass @@ -440,7 +452,7 @@ def update_H( Defaults to a zero tensor. """ - assert noise.shape == (2, 2) or (3, 3) + assert noise.shape in [(2, 2), (3, 3)] nqubits = omega.size(dim=0) a = torch.tensordot(omega * torch.cos(phi), Operators.sx, dims=0) From 13acac004ca567da43e9805b15278289cee3482f Mon Sep 17 00:00:00 2001 From: Kemal Bidzhiev Date: Thu, 26 Mar 2026 12:14:45 +0100 Subject: [PATCH 02/13] has right interaction --- emu_mps/hamiltonian.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/emu_mps/hamiltonian.py b/emu_mps/hamiltonian.py index 9401f4a1..2f5f98ea 100644 --- a/emu_mps/hamiltonian.py +++ b/emu_mps/hamiltonian.py @@ -85,13 +85,16 @@ def last_factor(self) -> torch.Tensor: """Return the MPO factor for the last site.""" pass + def _has_right_interaction(self, site: int) -> bool: + return bool(self.interaction_matrix[site, site + 1 :].any()) + class RydbergHamiltonianMPOFactors(HamiltonianMPOFactors): def first_factor(self) -> torch.Tensor: - has_right_interaction = self.interaction_matrix[0, 1:].any() - fac = torch.zeros( - 1, self.dim, self.dim, 3 if has_right_interaction else 2, dtype=dtype - ) + has_right_interaction = self._has_right_interaction(site=0) + bond_dim = 3 if has_right_interaction else 2 + fac = torch.zeros(1, self.dim, self.dim, bond_dim, dtype=dtype) + fac[0, :, :, 1] = self.identity if has_right_interaction: fac[0, :2, :2, 2] = Operators.n @@ -99,7 +102,7 @@ def first_factor(self) -> torch.Tensor: return fac def left_factor(self, n: int) -> torch.Tensor: - has_right_interaction = self.interaction_matrix[n, n + 1 :].any() + has_right_interaction = self._has_right_interaction(site=n) current_left_interactions = self.interaction_matrix[:n, n:].any(dim=1) left_interactions_to_keep = self.interaction_matrix[:n, n + 1 :].any(dim=1) @@ -212,7 +215,7 @@ def last_factor(self) -> torch.Tensor: class XYHamiltonianMPOFactors(HamiltonianMPOFactors): def first_factor(self) -> torch.Tensor: - has_right_interaction = self.interaction_matrix[0, 1:].any() + has_right_interaction = self._has_right_interaction(site=0) fac = torch.zeros( 1, self.dim, self.dim, 4 if has_right_interaction else 2, dtype=dtype ) @@ -224,7 +227,7 @@ def first_factor(self) -> torch.Tensor: return fac def left_factor(self, n: int) -> torch.Tensor: - has_right_interaction = self.interaction_matrix[n, n + 1 :].any() + has_right_interaction = self._has_right_interaction(site=n) current_left_interactions = self.interaction_matrix[:n, n:].any(dim=1) left_interactions_to_keep = self.interaction_matrix[:n, n + 1 :].any(dim=1) From 8e76ac660e05f93aeb0fd144e28d4161edfa12cd Mon Sep 17 00:00:00 2001 From: Kemal Bidzhiev Date: Thu, 26 Mar 2026 12:45:11 +0100 Subject: [PATCH 03/13] bond dims --- emu_mps/hamiltonian.py | 95 ++++++++++++++++++++++++++---------------- 1 file changed, 60 insertions(+), 35 deletions(-) diff --git a/emu_mps/hamiltonian.py b/emu_mps/hamiltonian.py index 2f5f98ea..f2eab6b7 100644 --- a/emu_mps/hamiltonian.py +++ b/emu_mps/hamiltonian.py @@ -88,12 +88,17 @@ def last_factor(self) -> torch.Tensor: def _has_right_interaction(self, site: int) -> bool: return bool(self.interaction_matrix[site, site + 1 :].any()) + def _has_left_interaction(self, site: int) -> bool: + return bool(self.interaction_matrix[site, :site].any()) + class RydbergHamiltonianMPOFactors(HamiltonianMPOFactors): def first_factor(self) -> torch.Tensor: has_right_interaction = self._has_right_interaction(site=0) - bond_dim = 3 if has_right_interaction else 2 - fac = torch.zeros(1, self.dim, self.dim, bond_dim, dtype=dtype) + + left_bond_dim = 1 + right_bond_dim = 3 if has_right_interaction else 2 + fac = torch.zeros(left_bond_dim, self.dim, self.dim, right_bond_dim, dtype=dtype) fac[0, :, :, 1] = self.identity if has_right_interaction: @@ -106,11 +111,15 @@ def left_factor(self, n: int) -> torch.Tensor: current_left_interactions = self.interaction_matrix[:n, n:].any(dim=1) left_interactions_to_keep = self.interaction_matrix[:n, n + 1 :].any(dim=1) + left_bond_dim = int(current_left_interactions.sum().item() + 2) + right_bond_dim = int( + left_interactions_to_keep.sum().item() + int(has_right_interaction) + 2 + ) fac = torch.zeros( - int(current_left_interactions.sum().item() + 2), + left_bond_dim, self.dim, self.dim, - int(left_interactions_to_keep.sum().item() + int(has_right_interaction) + 2), + right_bond_dim, dtype=dtype, ) @@ -138,11 +147,13 @@ def middle_factor(self) -> torch.Tensor: current_left_interactions = self.interaction_matrix[:n, n:].any(dim=1) current_right_interactions = self.interaction_matrix[n + 1 :, : n + 1].any(dim=1) + left_bond_dim = int(current_left_interactions.sum().item() + 2) + right_bond_dim = int(current_right_interactions.sum().item() + 2) fac = torch.zeros( - int(current_left_interactions.sum().item() + 2), + left_bond_dim, self.dim, self.dim, - int(current_right_interactions.sum().item() + 2), + right_bond_dim, dtype=dtype, ) @@ -168,15 +179,19 @@ def middle_factor(self) -> torch.Tensor: return fac def right_factor(self, n: int) -> torch.Tensor: - has_left_interaction = self.interaction_matrix[n, :n].any() + has_left_interaction = self._has_left_interaction(site=n) current_right_interactions = self.interaction_matrix[n + 1 :, : n + 1].any(dim=1) right_interactions_to_keep = self.interaction_matrix[n + 1 :, :n].any(dim=1) + left_bond_dim = int( + right_interactions_to_keep.sum().item() + int(has_left_interaction) + 2 + ) + right_bond_dim = int(current_right_interactions.sum().item() + 2) fac = torch.zeros( - int(right_interactions_to_keep.sum().item() + int(has_left_interaction) + 2), + left_bond_dim, self.dim, self.dim, - int(current_right_interactions.sum().item() + 2), + right_bond_dim, dtype=dtype, ) @@ -199,10 +214,11 @@ def right_factor(self, n: int) -> torch.Tensor: return fac def last_factor(self) -> torch.Tensor: - has_left_interaction = self.interaction_matrix[-1, :-1].any() - fac = torch.zeros( - 3 if has_left_interaction else 2, self.dim, self.dim, 1, dtype=dtype - ) + has_left_interaction = self._has_left_interaction(site=-1) + + left_bond_dim = 3 if has_left_interaction else 2 + right_bond_dim = 1 + fac = torch.zeros(left_bond_dim, self.dim, self.dim, right_bond_dim, dtype=dtype) fac[0, :, :, 0] = self.identity if has_left_interaction: if self.qubit_count >= 3: @@ -216,9 +232,10 @@ def last_factor(self) -> torch.Tensor: class XYHamiltonianMPOFactors(HamiltonianMPOFactors): def first_factor(self) -> torch.Tensor: has_right_interaction = self._has_right_interaction(site=0) - fac = torch.zeros( - 1, self.dim, self.dim, 4 if has_right_interaction else 2, dtype=dtype - ) + + left_bond_dim = 1 + right_bond_dim = 4 if has_right_interaction else 2 + fac = torch.zeros(left_bond_dim, self.dim, self.dim, right_bond_dim, dtype=dtype) fac[0, :, :, 1] = self.identity if has_right_interaction: fac[0, :2, :2, 2] = Operators.creation @@ -231,15 +248,17 @@ def left_factor(self, n: int) -> torch.Tensor: current_left_interactions = self.interaction_matrix[:n, n:].any(dim=1) left_interactions_to_keep = self.interaction_matrix[:n, n + 1 :].any(dim=1) + left_bond_dim = int(2 * current_left_interactions.sum().item() + 2) + right_bond_dim = int( + 2 * left_interactions_to_keep.sum().item() + + 2 * int(has_right_interaction) + + 2 + ) fac = torch.zeros( - int(2 * current_left_interactions.sum().item() + 2), + left_bond_dim, self.dim, self.dim, - int( - 2 * left_interactions_to_keep.sum().item() - + 2 * int(has_right_interaction) - + 2 - ), + right_bond_dim, dtype=dtype, ) @@ -273,11 +292,14 @@ def middle_factor(self) -> torch.Tensor: current_left_interactions = self.interaction_matrix[:n, n:].any(dim=1) current_right_interactions = self.interaction_matrix[n + 1 :, : n + 1].any(dim=1) + left_bond_dim = int(2 * current_left_interactions.sum().item() + 2) + right_bond_dim = int(2 * current_right_interactions.sum().item() + 2) + fac = torch.zeros( - int(2 * current_left_interactions.sum().item() + 2), + left_bond_dim, self.dim, self.dim, - int(2 * current_right_interactions.sum().item() + 2), + right_bond_dim, dtype=dtype, ) @@ -316,19 +338,21 @@ def middle_factor(self) -> torch.Tensor: return fac def right_factor(self, n: int) -> torch.Tensor: - has_left_interaction = self.interaction_matrix[n, :n].any() + has_left_interaction = self._has_left_interaction(site=n) current_right_interactions = self.interaction_matrix[n + 1 :, : n + 1].any(dim=1) right_interactions_to_keep = self.interaction_matrix[n + 1 :, :n].any(dim=1) + left_bond_dim = int( + 2 * right_interactions_to_keep.sum().item() + + 2 * int(has_left_interaction) + + 2 + ) + right_bond_dim = int(2 * current_right_interactions.sum().item() + 2) fac = torch.zeros( - int( - 2 * right_interactions_to_keep.sum().item() - + 2 * int(has_left_interaction) - + 2 - ), + left_bond_dim, self.dim, self.dim, - int(2 * current_right_interactions.sum().item() + 2), + right_bond_dim, dtype=dtype, ) @@ -356,10 +380,11 @@ def right_factor(self, n: int) -> torch.Tensor: return fac def last_factor(self) -> torch.Tensor: - has_left_interaction = self.interaction_matrix[-1, :-1].any() - fac = torch.zeros( - 4 if has_left_interaction else 2, self.dim, self.dim, 1, dtype=dtype - ) + has_left_interaction = self._has_left_interaction(site=-1) + + left_bond_dim = 4 if has_left_interaction else 2 + right_bond_dim = 1 + fac = torch.zeros(left_bond_dim, self.dim, self.dim, right_bond_dim, dtype=dtype) fac[0, :, :, 0] = self.identity if has_left_interaction: if self.qubit_count >= 3: From 1169b761273ca73498d6648c557d69940347d926 Mon Sep 17 00:00:00 2001 From: Kemal Bidzhiev Date: Thu, 26 Mar 2026 15:50:14 +0100 Subject: [PATCH 04/13] spsm to sxsy --- emu_mps/hamiltonian.py | 113 ++++++++++++++++++++++++++++++----------- 1 file changed, 84 insertions(+), 29 deletions(-) diff --git a/emu_mps/hamiltonian.py b/emu_mps/hamiltonian.py index f2eab6b7..d72940af 100644 --- a/emu_mps/hamiltonian.py +++ b/emu_mps/hamiltonian.py @@ -27,7 +27,7 @@ class Operators: class HamiltonianMPOFactors(ABC): - """Abstract generator for MPO factors of a two-body Hamiltonian. + """Abstract class for MPO factors of a two-body Hamiltonian. Subclasses implement the local MPO tensor at each position in the chain: first site, left half, middle, right half, and last site. @@ -238,8 +238,10 @@ def first_factor(self) -> torch.Tensor: fac = torch.zeros(left_bond_dim, self.dim, self.dim, right_bond_dim, dtype=dtype) fac[0, :, :, 1] = self.identity if has_right_interaction: - fac[0, :2, :2, 2] = Operators.creation - fac[0, :2, :2, 3] = Operators.creation.T + # fac[0, :2, :2, 2] = Operators.creation + # fac[0, :2, :2, 3] = Operators.creation.T + fac[0, :2, :2, 2] = Operators.sx + fac[0, :2, :2, 3] = Operators.sy return fac @@ -265,16 +267,28 @@ def left_factor(self, n: int) -> torch.Tensor: fac[0, :, :, 0] = self.identity fac[1, :, :, 1] = self.identity if has_right_interaction: - fac[1, :2, :2, -2] = Operators.creation - fac[1, :2, :2, -1] = Operators.creation.T - + # fac[1, :2, :2, -2] = Operators.creation + # fac[1, :2, :2, -1] = Operators.creation.T + fac[1, :2, :2, -2] = Operators.sx + fac[1, :2, :2, -1] = Operators.sy + + # fac[2::2, :2, :2, 0] = ( + # self.interaction_matrix[:n][current_left_interactions, n, None, None] + # * Operators.creation.T + # ) + # fac[3::2, :2, :2, 0] = ( + # self.interaction_matrix[:n][current_left_interactions, n, None, None] + # * Operators.creation + # ) fac[2::2, :2, :2, 0] = ( self.interaction_matrix[:n][current_left_interactions, n, None, None] - * Operators.creation.T + * 2 + * Operators.sx ) fac[3::2, :2, :2, 0] = ( self.interaction_matrix[:n][current_left_interactions, n, None, None] - * Operators.creation + * 2 + * Operators.sy ) i = 2 @@ -306,32 +320,55 @@ def middle_factor(self) -> torch.Tensor: fac[0, :, :, 0] = self.identity fac[1, :, :, 1] = self.identity + # fac[2::2, :2, :2, 0] = ( + # self.interaction_matrix[:n][current_left_interactions, n, None, None] + # * Operators.creation.T + # ) + # fac[3::2, :2, :2, 0] = ( + # self.interaction_matrix[:n][current_left_interactions, n, None, None] + # * Operators.creation + # ) + fac[2::2, :2, :2, 0] = ( self.interaction_matrix[:n][current_left_interactions, n, None, None] - * Operators.creation.T + * 2 + * Operators.sx ) fac[3::2, :2, :2, 0] = ( self.interaction_matrix[:n][current_left_interactions, n, None, None] - * Operators.creation + * 2 + * Operators.sy ) - fac[1, :2, :2, 2::2] = self.interaction_matrix[n + 1 :][ - None, None, current_right_interactions, n - ] * Operators.creation.unsqueeze(-1) - fac[1, :2, :2, 3::2] = self.interaction_matrix[n + 1 :][ - None, None, current_right_interactions, n - ] * Operators.creation.T.unsqueeze(-1) + # fac[1, :2, :2, 2::2] = self.interaction_matrix[n + 1 :][ + # None, None, current_right_interactions, n + # ] * Operators.creation.unsqueeze(-1) + # fac[1, :2, :2, 3::2] = self.interaction_matrix[n + 1 :][ + # None, None, current_right_interactions, n + # ] * Operators.creation.T.unsqueeze(-1) + fac[1, :2, :2, 2::2] = ( + self.interaction_matrix[n + 1 :][None, None, current_right_interactions, n] + * 2 + * Operators.sx.unsqueeze(-1) + ) + fac[1, :2, :2, 3::2] = ( + self.interaction_matrix[n + 1 :][None, None, current_right_interactions, n] + * 2 + * Operators.sy.unsqueeze(-1) + ) fac[2::2, :, :, 2::2] = ( self.interaction_matrix[:n, n + 1 :][current_left_interactions, :][ :, None, None, current_right_interactions ] + * 2 * self.identity[None, ..., None] ) fac[3::2, :, :, 3::2] = ( self.interaction_matrix[:n, n + 1 :][current_left_interactions, :][ :, None, None, current_right_interactions ] + * 2 * self.identity[None, ..., None] ) @@ -359,15 +396,28 @@ def right_factor(self, n: int) -> torch.Tensor: fac[0, :, :, 0] = self.identity fac[1, :, :, 1] = self.identity if has_left_interaction: - fac[2, :2, :2, 0] = Operators.creation.T - fac[3, :2, :2, 0] = Operators.creation - - fac[1, :2, :2, 2::2] = self.interaction_matrix[n + 1 :][ - None, None, current_right_interactions, n - ] * Operators.creation.unsqueeze(-1) - fac[1, :2, :2, 3::2] = self.interaction_matrix[n + 1 :][ - None, None, current_right_interactions, n - ] * Operators.creation.T.unsqueeze(-1) + # fac[2, :2, :2, 0] = Operators.creation.T + # fac[3, :2, :2, 0] = Operators.creation + fac[2, :2, :2, 0] = Operators.sx + fac[3, :2, :2, 0] = Operators.sy + + # fac[1, :2, :2, 2::2] = self.interaction_matrix[n + 1 :][ + # None, None, current_right_interactions, n + # ] * Operators.creation.unsqueeze(-1) + # fac[1, :2, :2, 3::2] = self.interaction_matrix[n + 1 :][ + # None, None, current_right_interactions, n + # ] * Operators.creation.T.unsqueeze(-1) + + fac[1, :2, :2, 2::2] = ( + self.interaction_matrix[n + 1 :][None, None, current_right_interactions, n] + * 2 + * Operators.sx.unsqueeze(-1) + ) + fac[1, :2, :2, 3::2] = ( + self.interaction_matrix[n + 1 :][None, None, current_right_interactions, n] + * 2 + * Operators.sy.unsqueeze(-1) + ) i = 4 if has_left_interaction else 2 j = 2 @@ -387,12 +437,17 @@ def last_factor(self) -> torch.Tensor: fac = torch.zeros(left_bond_dim, self.dim, self.dim, right_bond_dim, dtype=dtype) fac[0, :, :, 0] = self.identity if has_left_interaction: + # TODO improve if self.qubit_count >= 3: - fac[2, :2, :2, 0] = Operators.creation.T - fac[3, :2, :2, 0] = Operators.creation + # fac[2, :2, :2, 0] = Operators.creation.T + # fac[3, :2, :2, 0] = Operators.creation + fac[2, :2, :2, 0] = Operators.sx + fac[3, :2, :2, 0] = Operators.sy else: - fac[2, :2, :2, 0] = self.interaction_matrix[0, 1] * Operators.creation.T - fac[3, :2, :2, 0] = self.interaction_matrix[0, 1] * Operators.creation + # fac[2, :2, :2, 0] = self.interaction_matrix[0, 1] * Operators.creation.T + # fac[3, :2, :2, 0] = self.interaction_matrix[0, 1] * Operators.creation + fac[2, :2, :2, 0] = self.interaction_matrix[0, 1] * 2 * Operators.sx + fac[3, :2, :2, 0] = self.interaction_matrix[0, 1] * 2 * Operators.sy return fac From 9b4446421930a80c4187c8e602a31df9c8b374d3 Mon Sep 17 00:00:00 2001 From: Kemal Bidzhiev Date: Thu, 26 Mar 2026 16:02:10 +0100 Subject: [PATCH 05/13] shrinking the code --- emu_mps/hamiltonian.py | 120 +++++++++-------------------------------- 1 file changed, 24 insertions(+), 96 deletions(-) diff --git a/emu_mps/hamiltonian.py b/emu_mps/hamiltonian.py index d72940af..a4d56e63 100644 --- a/emu_mps/hamiltonian.py +++ b/emu_mps/hamiltonian.py @@ -238,8 +238,6 @@ def first_factor(self) -> torch.Tensor: fac = torch.zeros(left_bond_dim, self.dim, self.dim, right_bond_dim, dtype=dtype) fac[0, :, :, 1] = self.identity if has_right_interaction: - # fac[0, :2, :2, 2] = Operators.creation - # fac[0, :2, :2, 3] = Operators.creation.T fac[0, :2, :2, 2] = Operators.sx fac[0, :2, :2, 3] = Operators.sy @@ -267,29 +265,12 @@ def left_factor(self, n: int) -> torch.Tensor: fac[0, :, :, 0] = self.identity fac[1, :, :, 1] = self.identity if has_right_interaction: - # fac[1, :2, :2, -2] = Operators.creation - # fac[1, :2, :2, -1] = Operators.creation.T fac[1, :2, :2, -2] = Operators.sx fac[1, :2, :2, -1] = Operators.sy - # fac[2::2, :2, :2, 0] = ( - # self.interaction_matrix[:n][current_left_interactions, n, None, None] - # * Operators.creation.T - # ) - # fac[3::2, :2, :2, 0] = ( - # self.interaction_matrix[:n][current_left_interactions, n, None, None] - # * Operators.creation - # ) - fac[2::2, :2, :2, 0] = ( - self.interaction_matrix[:n][current_left_interactions, n, None, None] - * 2 - * Operators.sx - ) - fac[3::2, :2, :2, 0] = ( - self.interaction_matrix[:n][current_left_interactions, n, None, None] - * 2 - * Operators.sy - ) + coeff = 2 * self.interaction_matrix[:n][current_left_interactions, n, None, None] + fac[2::2, :2, :2, 0] = coeff * Operators.sx + fac[3::2, :2, :2, 0] = coeff * Operators.sy i = 2 j = 2 @@ -320,57 +301,25 @@ def middle_factor(self) -> torch.Tensor: fac[0, :, :, 0] = self.identity fac[1, :, :, 1] = self.identity - # fac[2::2, :2, :2, 0] = ( - # self.interaction_matrix[:n][current_left_interactions, n, None, None] - # * Operators.creation.T - # ) - # fac[3::2, :2, :2, 0] = ( - # self.interaction_matrix[:n][current_left_interactions, n, None, None] - # * Operators.creation - # ) - - fac[2::2, :2, :2, 0] = ( - self.interaction_matrix[:n][current_left_interactions, n, None, None] - * 2 - * Operators.sx - ) - fac[3::2, :2, :2, 0] = ( - self.interaction_matrix[:n][current_left_interactions, n, None, None] - * 2 - * Operators.sy - ) + coeff = 2 * self.interaction_matrix[:n][current_left_interactions, n, None, None] + fac[2::2, :2, :2, 0] = coeff * Operators.sx + fac[3::2, :2, :2, 0] = coeff * Operators.sy - # fac[1, :2, :2, 2::2] = self.interaction_matrix[n + 1 :][ - # None, None, current_right_interactions, n - # ] * Operators.creation.unsqueeze(-1) - # fac[1, :2, :2, 3::2] = self.interaction_matrix[n + 1 :][ - # None, None, current_right_interactions, n - # ] * Operators.creation.T.unsqueeze(-1) - fac[1, :2, :2, 2::2] = ( - self.interaction_matrix[n + 1 :][None, None, current_right_interactions, n] - * 2 - * Operators.sx.unsqueeze(-1) - ) - fac[1, :2, :2, 3::2] = ( - self.interaction_matrix[n + 1 :][None, None, current_right_interactions, n] - * 2 - * Operators.sy.unsqueeze(-1) + coeff = ( + 2 + * self.interaction_matrix[n + 1 :][None, None, current_right_interactions, n] ) + fac[1, :2, :2, 2::2] = coeff * Operators.sx.unsqueeze(-1) + fac[1, :2, :2, 3::2] = coeff * Operators.sy.unsqueeze(-1) - fac[2::2, :, :, 2::2] = ( - self.interaction_matrix[:n, n + 1 :][current_left_interactions, :][ - :, None, None, current_right_interactions - ] - * 2 - * self.identity[None, ..., None] - ) - fac[3::2, :, :, 3::2] = ( - self.interaction_matrix[:n, n + 1 :][current_left_interactions, :][ + coeff = ( + 2 + * self.interaction_matrix[:n, n + 1 :][current_left_interactions, :][ :, None, None, current_right_interactions ] - * 2 - * self.identity[None, ..., None] ) + fac[2::2, :, :, 2::2] = coeff * self.identity[None, ..., None] + fac[3::2, :, :, 3::2] = coeff * self.identity[None, ..., None] return fac @@ -396,28 +345,15 @@ def right_factor(self, n: int) -> torch.Tensor: fac[0, :, :, 0] = self.identity fac[1, :, :, 1] = self.identity if has_left_interaction: - # fac[2, :2, :2, 0] = Operators.creation.T - # fac[3, :2, :2, 0] = Operators.creation fac[2, :2, :2, 0] = Operators.sx fac[3, :2, :2, 0] = Operators.sy - # fac[1, :2, :2, 2::2] = self.interaction_matrix[n + 1 :][ - # None, None, current_right_interactions, n - # ] * Operators.creation.unsqueeze(-1) - # fac[1, :2, :2, 3::2] = self.interaction_matrix[n + 1 :][ - # None, None, current_right_interactions, n - # ] * Operators.creation.T.unsqueeze(-1) - - fac[1, :2, :2, 2::2] = ( - self.interaction_matrix[n + 1 :][None, None, current_right_interactions, n] - * 2 - * Operators.sx.unsqueeze(-1) - ) - fac[1, :2, :2, 3::2] = ( - self.interaction_matrix[n + 1 :][None, None, current_right_interactions, n] - * 2 - * Operators.sy.unsqueeze(-1) + coeff = ( + 2 + * self.interaction_matrix[n + 1 :][None, None, current_right_interactions, n] ) + fac[1, :2, :2, 2::2] = coeff * Operators.sx.unsqueeze(-1) + fac[1, :2, :2, 3::2] = coeff * Operators.sy.unsqueeze(-1) i = 4 if has_left_interaction else 2 j = 2 @@ -437,17 +373,9 @@ def last_factor(self) -> torch.Tensor: fac = torch.zeros(left_bond_dim, self.dim, self.dim, right_bond_dim, dtype=dtype) fac[0, :, :, 0] = self.identity if has_left_interaction: - # TODO improve - if self.qubit_count >= 3: - # fac[2, :2, :2, 0] = Operators.creation.T - # fac[3, :2, :2, 0] = Operators.creation - fac[2, :2, :2, 0] = Operators.sx - fac[3, :2, :2, 0] = Operators.sy - else: - # fac[2, :2, :2, 0] = self.interaction_matrix[0, 1] * Operators.creation.T - # fac[3, :2, :2, 0] = self.interaction_matrix[0, 1] * Operators.creation - fac[2, :2, :2, 0] = self.interaction_matrix[0, 1] * 2 * Operators.sx - fac[3, :2, :2, 0] = self.interaction_matrix[0, 1] * 2 * Operators.sy + coeff = 1 if self.qubit_count >= 3 else 2 * self.interaction_matrix[0, 1] + fac[2, :2, :2, 0] = coeff * Operators.sx + fac[3, :2, :2, 0] = coeff * Operators.sy return fac From dd66dad06f4e7fe443ef936c6a3c4d9de5a5506b Mon Sep 17 00:00:00 2001 From: Kemal Bidzhiev Date: Thu, 26 Mar 2026 16:30:56 +0100 Subject: [PATCH 06/13] shrinking the code in rydberg --- emu_mps/hamiltonian.py | 87 ++++++++++++++++++------------------------ 1 file changed, 38 insertions(+), 49 deletions(-) diff --git a/emu_mps/hamiltonian.py b/emu_mps/hamiltonian.py index a4d56e63..c9120a5b 100644 --- a/emu_mps/hamiltonian.py +++ b/emu_mps/hamiltonian.py @@ -128,10 +128,8 @@ def left_factor(self, n: int) -> torch.Tensor: if has_right_interaction: fac[1, :2, :2, -1] = Operators.n - fac[2:, :2, :2, 0] = ( - self.interaction_matrix[:n][current_left_interactions, n, None, None] - * Operators.n - ) + coeff = self.interaction_matrix[:n][current_left_interactions, n, None, None] + fac[2:, :2, :2, 0] = coeff * Operators.n i = 2 j = 2 @@ -160,21 +158,18 @@ def middle_factor(self) -> torch.Tensor: fac[0, :, :, 0] = self.identity fac[1, :, :, 1] = self.identity - fac[2:, :2, :2, 0] = ( - self.interaction_matrix[:n][current_left_interactions, n, None, None] - * Operators.n - ) + coeff = self.interaction_matrix[:n][current_left_interactions, n, None, None] + fac[2:, :2, :2, 0] = coeff * Operators.n - fac[1, :2, :2, 2:] = self.interaction_matrix[n + 1 :][ + coeff = self.interaction_matrix[n + 1 :][ None, None, current_right_interactions, n - ] * Operators.n.unsqueeze(-1) + ] + fac[1, :2, :2, 2:] = coeff * Operators.n.unsqueeze(-1) - fac[2:, :, :, 2:] = ( - self.interaction_matrix[:n, n + 1 :][current_left_interactions, :][ - :, None, None, current_right_interactions - ] - * self.identity[None, ..., None] - ) + coeff = self.interaction_matrix[:n, n + 1 :][current_left_interactions, :][ + :, None, None, current_right_interactions + ] + fac[2:, :, :, 2:] = coeff * self.identity[None, ..., None] return fac @@ -200,9 +195,10 @@ def right_factor(self, n: int) -> torch.Tensor: if has_left_interaction: fac[2, :2, :2, 0] = Operators.n - fac[1, :2, :2, 2:] = self.interaction_matrix[n + 1 :][ + coeff = self.interaction_matrix[n + 1 :][ None, None, current_right_interactions, n - ] * Operators.n.unsqueeze(-1) + ] + fac[1, :2, :2, 2:] = coeff * Operators.n.unsqueeze(-1) i = 3 if has_left_interaction else 2 j = 2 @@ -221,10 +217,8 @@ def last_factor(self) -> torch.Tensor: fac = torch.zeros(left_bond_dim, self.dim, self.dim, right_bond_dim, dtype=dtype) fac[0, :, :, 0] = self.identity if has_left_interaction: - if self.qubit_count >= 3: - fac[2, :2, :2, 0] = Operators.n - else: - fac[2, :2, :2, 0] = self.interaction_matrix[0, 1] * Operators.n + coeff = 1 if self.qubit_count >= 3 else self.interaction_matrix[0, 1] + fac[2, :2, :2, 0] = coeff * Operators.n return fac @@ -268,9 +262,9 @@ def left_factor(self, n: int) -> torch.Tensor: fac[1, :2, :2, -2] = Operators.sx fac[1, :2, :2, -1] = Operators.sy - coeff = 2 * self.interaction_matrix[:n][current_left_interactions, n, None, None] - fac[2::2, :2, :2, 0] = coeff * Operators.sx - fac[3::2, :2, :2, 0] = coeff * Operators.sy + coeff = self.interaction_matrix[:n][current_left_interactions, n, None, None] + fac[2::2, :2, :2, 0] = coeff * 2 * Operators.sx + fac[3::2, :2, :2, 0] = coeff * 2 * Operators.sy i = 2 j = 2 @@ -301,25 +295,21 @@ def middle_factor(self) -> torch.Tensor: fac[0, :, :, 0] = self.identity fac[1, :, :, 1] = self.identity - coeff = 2 * self.interaction_matrix[:n][current_left_interactions, n, None, None] - fac[2::2, :2, :2, 0] = coeff * Operators.sx - fac[3::2, :2, :2, 0] = coeff * Operators.sy + coeff = self.interaction_matrix[:n][current_left_interactions, n, None, None] + fac[2::2, :2, :2, 0] = coeff * 2 * Operators.sx + fac[3::2, :2, :2, 0] = coeff * 2 * Operators.sy - coeff = ( - 2 - * self.interaction_matrix[n + 1 :][None, None, current_right_interactions, n] - ) - fac[1, :2, :2, 2::2] = coeff * Operators.sx.unsqueeze(-1) - fac[1, :2, :2, 3::2] = coeff * Operators.sy.unsqueeze(-1) - - coeff = ( - 2 - * self.interaction_matrix[:n, n + 1 :][current_left_interactions, :][ - :, None, None, current_right_interactions - ] - ) - fac[2::2, :, :, 2::2] = coeff * self.identity[None, ..., None] - fac[3::2, :, :, 3::2] = coeff * self.identity[None, ..., None] + coeff = self.interaction_matrix[n + 1 :][ + None, None, current_right_interactions, n + ] + fac[1, :2, :2, 2::2] = coeff * 2 * Operators.sx.unsqueeze(-1) + fac[1, :2, :2, 3::2] = coeff * 2 * Operators.sy.unsqueeze(-1) + + coeff = self.interaction_matrix[:n, n + 1 :][current_left_interactions, :][ + :, None, None, current_right_interactions + ] + fac[2::2, :, :, 2::2] = coeff * 2 * self.identity[None, ..., None] + fac[3::2, :, :, 3::2] = coeff * 2 * self.identity[None, ..., None] return fac @@ -348,12 +338,11 @@ def right_factor(self, n: int) -> torch.Tensor: fac[2, :2, :2, 0] = Operators.sx fac[3, :2, :2, 0] = Operators.sy - coeff = ( - 2 - * self.interaction_matrix[n + 1 :][None, None, current_right_interactions, n] - ) - fac[1, :2, :2, 2::2] = coeff * Operators.sx.unsqueeze(-1) - fac[1, :2, :2, 3::2] = coeff * Operators.sy.unsqueeze(-1) + coeff = self.interaction_matrix[n + 1 :][ + None, None, current_right_interactions, n + ] + fac[1, :2, :2, 2::2] = coeff * 2 * Operators.sx.unsqueeze(-1) + fac[1, :2, :2, 3::2] = coeff * 2 * Operators.sy.unsqueeze(-1) i = 4 if has_left_interaction else 2 j = 2 From 86be09f200def2ecc151ad4860b2b3112d17fd1e Mon Sep 17 00:00:00 2001 From: Kemal Bidzhiev Date: Thu, 26 Mar 2026 16:42:42 +0100 Subject: [PATCH 07/13] empty factor --- emu_mps/hamiltonian.py | 69 +++++++++++++----------------------------- 1 file changed, 21 insertions(+), 48 deletions(-) diff --git a/emu_mps/hamiltonian.py b/emu_mps/hamiltonian.py index c9120a5b..1311ce15 100644 --- a/emu_mps/hamiltonian.py +++ b/emu_mps/hamiltonian.py @@ -91,6 +91,15 @@ def _has_right_interaction(self, site: int) -> bool: def _has_left_interaction(self, site: int) -> bool: return bool(self.interaction_matrix[site, :site].any()) + def _empty_factor(self, left_bond_dim: int, right_bond_dim: int) -> torch.Tensor: + return torch.zeros( + left_bond_dim, + self.dim, + self.dim, + right_bond_dim, + dtype=dtype, + ) + class RydbergHamiltonianMPOFactors(HamiltonianMPOFactors): def first_factor(self) -> torch.Tensor: @@ -98,7 +107,7 @@ def first_factor(self) -> torch.Tensor: left_bond_dim = 1 right_bond_dim = 3 if has_right_interaction else 2 - fac = torch.zeros(left_bond_dim, self.dim, self.dim, right_bond_dim, dtype=dtype) + fac = self._empty_factor(left_bond_dim, right_bond_dim) fac[0, :, :, 1] = self.identity if has_right_interaction: @@ -115,13 +124,7 @@ def left_factor(self, n: int) -> torch.Tensor: right_bond_dim = int( left_interactions_to_keep.sum().item() + int(has_right_interaction) + 2 ) - fac = torch.zeros( - left_bond_dim, - self.dim, - self.dim, - right_bond_dim, - dtype=dtype, - ) + fac = self._empty_factor(left_bond_dim, right_bond_dim) fac[0, :, :, 0] = self.identity fac[1, :, :, 1] = self.identity @@ -147,13 +150,7 @@ def middle_factor(self) -> torch.Tensor: left_bond_dim = int(current_left_interactions.sum().item() + 2) right_bond_dim = int(current_right_interactions.sum().item() + 2) - fac = torch.zeros( - left_bond_dim, - self.dim, - self.dim, - right_bond_dim, - dtype=dtype, - ) + fac = self._empty_factor(left_bond_dim, right_bond_dim) fac[0, :, :, 0] = self.identity fac[1, :, :, 1] = self.identity @@ -182,13 +179,7 @@ def right_factor(self, n: int) -> torch.Tensor: right_interactions_to_keep.sum().item() + int(has_left_interaction) + 2 ) right_bond_dim = int(current_right_interactions.sum().item() + 2) - fac = torch.zeros( - left_bond_dim, - self.dim, - self.dim, - right_bond_dim, - dtype=dtype, - ) + fac = self._empty_factor(left_bond_dim, right_bond_dim) fac[0, :, :, 0] = self.identity fac[1, :, :, 1] = self.identity @@ -214,10 +205,10 @@ def last_factor(self) -> torch.Tensor: left_bond_dim = 3 if has_left_interaction else 2 right_bond_dim = 1 - fac = torch.zeros(left_bond_dim, self.dim, self.dim, right_bond_dim, dtype=dtype) + fac = self._empty_factor(left_bond_dim, right_bond_dim) fac[0, :, :, 0] = self.identity if has_left_interaction: - coeff = 1 if self.qubit_count >= 3 else self.interaction_matrix[0, 1] + coeff = self.interaction_matrix[0, 1] if self.qubit_count == 2 else 1 fac[2, :2, :2, 0] = coeff * Operators.n return fac @@ -229,7 +220,7 @@ def first_factor(self) -> torch.Tensor: left_bond_dim = 1 right_bond_dim = 4 if has_right_interaction else 2 - fac = torch.zeros(left_bond_dim, self.dim, self.dim, right_bond_dim, dtype=dtype) + fac = self._empty_factor(left_bond_dim, right_bond_dim) fac[0, :, :, 1] = self.identity if has_right_interaction: fac[0, :2, :2, 2] = Operators.sx @@ -248,13 +239,7 @@ def left_factor(self, n: int) -> torch.Tensor: + 2 * int(has_right_interaction) + 2 ) - fac = torch.zeros( - left_bond_dim, - self.dim, - self.dim, - right_bond_dim, - dtype=dtype, - ) + fac = self._empty_factor(left_bond_dim, right_bond_dim) fac[0, :, :, 0] = self.identity fac[1, :, :, 1] = self.identity @@ -284,13 +269,7 @@ def middle_factor(self) -> torch.Tensor: left_bond_dim = int(2 * current_left_interactions.sum().item() + 2) right_bond_dim = int(2 * current_right_interactions.sum().item() + 2) - fac = torch.zeros( - left_bond_dim, - self.dim, - self.dim, - right_bond_dim, - dtype=dtype, - ) + fac = self._empty_factor(left_bond_dim, right_bond_dim) fac[0, :, :, 0] = self.identity fac[1, :, :, 1] = self.identity @@ -324,13 +303,7 @@ def right_factor(self, n: int) -> torch.Tensor: + 2 ) right_bond_dim = int(2 * current_right_interactions.sum().item() + 2) - fac = torch.zeros( - left_bond_dim, - self.dim, - self.dim, - right_bond_dim, - dtype=dtype, - ) + fac = self._empty_factor(left_bond_dim, right_bond_dim) fac[0, :, :, 0] = self.identity fac[1, :, :, 1] = self.identity @@ -359,10 +332,10 @@ def last_factor(self) -> torch.Tensor: left_bond_dim = 4 if has_left_interaction else 2 right_bond_dim = 1 - fac = torch.zeros(left_bond_dim, self.dim, self.dim, right_bond_dim, dtype=dtype) + fac = self._empty_factor(left_bond_dim, right_bond_dim) fac[0, :, :, 0] = self.identity if has_left_interaction: - coeff = 1 if self.qubit_count >= 3 else 2 * self.interaction_matrix[0, 1] + coeff = 2 * self.interaction_matrix[0, 1] if self.qubit_count == 2 else 1 fac[2, :2, :2, 0] = coeff * Operators.sx fac[3, :2, :2, 0] = coeff * Operators.sy From c2f5aa1171a6851ca0c0a10741644f6b40dce8cc Mon Sep 17 00:00:00 2001 From: Kemal Bidzhiev Date: Thu, 26 Mar 2026 17:08:25 +0100 Subject: [PATCH 08/13] helpers --- emu_mps/hamiltonian.py | 190 ++++++++++++++++++++++++----------------- 1 file changed, 110 insertions(+), 80 deletions(-) diff --git a/emu_mps/hamiltonian.py b/emu_mps/hamiltonian.py index 1311ce15..599d4ca7 100644 --- a/emu_mps/hamiltonian.py +++ b/emu_mps/hamiltonian.py @@ -42,20 +42,20 @@ def __init__(self, interaction_matrix: torch.Tensor, dim: int = 2): self.interaction_matrix = interaction_matrix.clone() self.interaction_matrix.fill_diagonal_(0.0) # or assert self.qubit_count = self.interaction_matrix.shape[0] - self.middle = self.qubit_count // 2 + self.middle_site = self.qubit_count // 2 self.identity = Operators.id if self.dim == 2 else Operators.id_3x3 def __iter__(self) -> Iterator[torch.Tensor]: """Yield the full ordered list of MPO factors for the Hamiltonian.""" yield self.first_factor() - for n in range(1, self.middle): + for n in range(1, self.middle_site): yield self.left_factor(n) if self.qubit_count >= 3: yield self.middle_factor() - for n in range(self.middle + 1, self.qubit_count - 1): + for n in range(self.middle_site + 1, self.qubit_count - 1): yield self.right_factor(n) yield self.last_factor() @@ -100,6 +100,36 @@ def _empty_factor(self, left_bond_dim: int, right_bond_dim: int) -> torch.Tensor dtype=dtype, ) + def _left_interaction_masks(self, n: int) -> tuple[torch.Tensor, torch.Tensor]: + current_left_interactions = self.interaction_matrix[:n, n:].any(dim=1) + left_interactions_to_keep = self.interaction_matrix[:n, n + 1 :].any(dim=1) + return current_left_interactions, left_interactions_to_keep + + def _right_interaction_masks(self, n: int) -> tuple[torch.Tensor, torch.Tensor]: + current_right_interactions = self.interaction_matrix[n + 1 :, : n + 1].any(dim=1) + right_interactions_to_keep = self.interaction_matrix[n + 1 :, :n].any(dim=1) + return current_right_interactions, right_interactions_to_keep + + def _left_interaction_coefficients( + self, n: int, current_left_interactions: torch.Tensor + ) -> torch.Tensor: + return self.interaction_matrix[:n][current_left_interactions, n, None, None] + + def _right_interaction_coefficients( + self, n: int, current_right_interactions: torch.Tensor + ) -> torch.Tensor: + return self.interaction_matrix[n + 1 :][None, None, current_right_interactions, n] + + def _middle_interaction_coefficients( + self, + n: int, + current_left_interactions: torch.Tensor, + current_right_interactions: torch.Tensor, + ) -> torch.Tensor: + return self.interaction_matrix[:n, n + 1 :][current_left_interactions, :][ + :, None, None, current_right_interactions + ] + class RydbergHamiltonianMPOFactors(HamiltonianMPOFactors): def first_factor(self) -> torch.Tensor: @@ -107,13 +137,13 @@ def first_factor(self) -> torch.Tensor: left_bond_dim = 1 right_bond_dim = 3 if has_right_interaction else 2 - fac = self._empty_factor(left_bond_dim, right_bond_dim) + factor = self._empty_factor(left_bond_dim, right_bond_dim) - fac[0, :, :, 1] = self.identity + factor[0, :, :, 1] = self.identity if has_right_interaction: - fac[0, :2, :2, 2] = Operators.n + factor[0, :2, :2, 2] = Operators.n - return fac + return factor def left_factor(self, n: int) -> torch.Tensor: has_right_interaction = self._has_right_interaction(site=n) @@ -124,51 +154,51 @@ def left_factor(self, n: int) -> torch.Tensor: right_bond_dim = int( left_interactions_to_keep.sum().item() + int(has_right_interaction) + 2 ) - fac = self._empty_factor(left_bond_dim, right_bond_dim) + factor = self._empty_factor(left_bond_dim, right_bond_dim) - fac[0, :, :, 0] = self.identity - fac[1, :, :, 1] = self.identity + factor[0, :, :, 0] = self.identity + factor[1, :, :, 1] = self.identity if has_right_interaction: - fac[1, :2, :2, -1] = Operators.n + factor[1, :2, :2, -1] = Operators.n coeff = self.interaction_matrix[:n][current_left_interactions, n, None, None] - fac[2:, :2, :2, 0] = coeff * Operators.n + factor[2:, :2, :2, 0] = coeff * Operators.n i = 2 j = 2 for current_left_interaction in current_left_interactions.nonzero().flatten(): if left_interactions_to_keep[current_left_interaction]: - fac[i, :, :, j] = self.identity + factor[i, :, :, j] = self.identity j += 1 i += 1 - return fac + return factor def middle_factor(self) -> torch.Tensor: - n = self.middle + n = self.middle_site current_left_interactions = self.interaction_matrix[:n, n:].any(dim=1) current_right_interactions = self.interaction_matrix[n + 1 :, : n + 1].any(dim=1) left_bond_dim = int(current_left_interactions.sum().item() + 2) right_bond_dim = int(current_right_interactions.sum().item() + 2) - fac = self._empty_factor(left_bond_dim, right_bond_dim) + factor = self._empty_factor(left_bond_dim, right_bond_dim) - fac[0, :, :, 0] = self.identity - fac[1, :, :, 1] = self.identity + factor[0, :, :, 0] = self.identity + factor[1, :, :, 1] = self.identity coeff = self.interaction_matrix[:n][current_left_interactions, n, None, None] - fac[2:, :2, :2, 0] = coeff * Operators.n + factor[2:, :2, :2, 0] = coeff * Operators.n coeff = self.interaction_matrix[n + 1 :][ None, None, current_right_interactions, n ] - fac[1, :2, :2, 2:] = coeff * Operators.n.unsqueeze(-1) + factor[1, :2, :2, 2:] = coeff * Operators.n.unsqueeze(-1) - coeff = self.interaction_matrix[:n, n + 1 :][current_left_interactions, :][ - :, None, None, current_right_interactions - ] - fac[2:, :, :, 2:] = coeff * self.identity[None, ..., None] + coeff = self._middle_interaction_coefficients( + n, current_left_interactions, current_right_interactions + ) + factor[2:, :, :, 2:] = coeff * self.identity[None, ..., None] - return fac + return factor def right_factor(self, n: int) -> torch.Tensor: has_left_interaction = self._has_left_interaction(site=n) @@ -179,39 +209,39 @@ def right_factor(self, n: int) -> torch.Tensor: right_interactions_to_keep.sum().item() + int(has_left_interaction) + 2 ) right_bond_dim = int(current_right_interactions.sum().item() + 2) - fac = self._empty_factor(left_bond_dim, right_bond_dim) + factor = self._empty_factor(left_bond_dim, right_bond_dim) - fac[0, :, :, 0] = self.identity - fac[1, :, :, 1] = self.identity + factor[0, :, :, 0] = self.identity + factor[1, :, :, 1] = self.identity if has_left_interaction: - fac[2, :2, :2, 0] = Operators.n + factor[2, :2, :2, 0] = Operators.n coeff = self.interaction_matrix[n + 1 :][ None, None, current_right_interactions, n ] - fac[1, :2, :2, 2:] = coeff * Operators.n.unsqueeze(-1) + factor[1, :2, :2, 2:] = coeff * Operators.n.unsqueeze(-1) i = 3 if has_left_interaction else 2 j = 2 for current_right_interaction in current_right_interactions.nonzero().flatten(): if right_interactions_to_keep[current_right_interaction]: - fac[i, :, :, j] = self.identity + factor[i, :, :, j] = self.identity i += 1 j += 1 - return fac + return factor def last_factor(self) -> torch.Tensor: has_left_interaction = self._has_left_interaction(site=-1) left_bond_dim = 3 if has_left_interaction else 2 right_bond_dim = 1 - fac = self._empty_factor(left_bond_dim, right_bond_dim) - fac[0, :, :, 0] = self.identity + factor = self._empty_factor(left_bond_dim, right_bond_dim) + factor[0, :, :, 0] = self.identity if has_left_interaction: coeff = self.interaction_matrix[0, 1] if self.qubit_count == 2 else 1 - fac[2, :2, :2, 0] = coeff * Operators.n + factor[2, :2, :2, 0] = coeff * Operators.n - return fac + return factor class XYHamiltonianMPOFactors(HamiltonianMPOFactors): @@ -220,13 +250,13 @@ def first_factor(self) -> torch.Tensor: left_bond_dim = 1 right_bond_dim = 4 if has_right_interaction else 2 - fac = self._empty_factor(left_bond_dim, right_bond_dim) - fac[0, :, :, 1] = self.identity + factor = self._empty_factor(left_bond_dim, right_bond_dim) + factor[0, :, :, 1] = self.identity if has_right_interaction: - fac[0, :2, :2, 2] = Operators.sx - fac[0, :2, :2, 3] = Operators.sy + factor[0, :2, :2, 2] = Operators.sx + factor[0, :2, :2, 3] = Operators.sy - return fac + return factor def left_factor(self, n: int) -> torch.Tensor: has_right_interaction = self._has_right_interaction(site=n) @@ -239,58 +269,58 @@ def left_factor(self, n: int) -> torch.Tensor: + 2 * int(has_right_interaction) + 2 ) - fac = self._empty_factor(left_bond_dim, right_bond_dim) + factor = self._empty_factor(left_bond_dim, right_bond_dim) - fac[0, :, :, 0] = self.identity - fac[1, :, :, 1] = self.identity + factor[0, :, :, 0] = self.identity + factor[1, :, :, 1] = self.identity if has_right_interaction: - fac[1, :2, :2, -2] = Operators.sx - fac[1, :2, :2, -1] = Operators.sy + factor[1, :2, :2, -2] = Operators.sx + factor[1, :2, :2, -1] = Operators.sy coeff = self.interaction_matrix[:n][current_left_interactions, n, None, None] - fac[2::2, :2, :2, 0] = coeff * 2 * Operators.sx - fac[3::2, :2, :2, 0] = coeff * 2 * Operators.sy + factor[2::2, :2, :2, 0] = coeff * 2 * Operators.sx + factor[3::2, :2, :2, 0] = coeff * 2 * Operators.sy i = 2 j = 2 for current_left_interaction in current_left_interactions.nonzero().flatten(): if left_interactions_to_keep[current_left_interaction]: - fac[i, :, :, j] = self.identity - fac[i + 1, :, :, j + 1] = self.identity + factor[i, :, :, j] = self.identity + factor[i + 1, :, :, j + 1] = self.identity j += 2 i += 2 - return fac + return factor def middle_factor(self) -> torch.Tensor: - n = self.middle + n = self.middle_site current_left_interactions = self.interaction_matrix[:n, n:].any(dim=1) current_right_interactions = self.interaction_matrix[n + 1 :, : n + 1].any(dim=1) left_bond_dim = int(2 * current_left_interactions.sum().item() + 2) right_bond_dim = int(2 * current_right_interactions.sum().item() + 2) - fac = self._empty_factor(left_bond_dim, right_bond_dim) + factor = self._empty_factor(left_bond_dim, right_bond_dim) - fac[0, :, :, 0] = self.identity - fac[1, :, :, 1] = self.identity + factor[0, :, :, 0] = self.identity + factor[1, :, :, 1] = self.identity coeff = self.interaction_matrix[:n][current_left_interactions, n, None, None] - fac[2::2, :2, :2, 0] = coeff * 2 * Operators.sx - fac[3::2, :2, :2, 0] = coeff * 2 * Operators.sy + factor[2::2, :2, :2, 0] = coeff * 2 * Operators.sx + factor[3::2, :2, :2, 0] = coeff * 2 * Operators.sy coeff = self.interaction_matrix[n + 1 :][ None, None, current_right_interactions, n ] - fac[1, :2, :2, 2::2] = coeff * 2 * Operators.sx.unsqueeze(-1) - fac[1, :2, :2, 3::2] = coeff * 2 * Operators.sy.unsqueeze(-1) + factor[1, :2, :2, 2::2] = coeff * 2 * Operators.sx.unsqueeze(-1) + factor[1, :2, :2, 3::2] = coeff * 2 * Operators.sy.unsqueeze(-1) - coeff = self.interaction_matrix[:n, n + 1 :][current_left_interactions, :][ - :, None, None, current_right_interactions - ] - fac[2::2, :, :, 2::2] = coeff * 2 * self.identity[None, ..., None] - fac[3::2, :, :, 3::2] = coeff * 2 * self.identity[None, ..., None] + coeff = self._middle_interaction_coefficients( + n, current_left_interactions, current_right_interactions + ) + factor[2::2, :, :, 2::2] = coeff * 2 * self.identity[None, ..., None] + factor[3::2, :, :, 3::2] = coeff * 2 * self.identity[None, ..., None] - return fac + return factor def right_factor(self, n: int) -> torch.Tensor: has_left_interaction = self._has_left_interaction(site=n) @@ -303,43 +333,43 @@ def right_factor(self, n: int) -> torch.Tensor: + 2 ) right_bond_dim = int(2 * current_right_interactions.sum().item() + 2) - fac = self._empty_factor(left_bond_dim, right_bond_dim) + factor = self._empty_factor(left_bond_dim, right_bond_dim) - fac[0, :, :, 0] = self.identity - fac[1, :, :, 1] = self.identity + factor[0, :, :, 0] = self.identity + factor[1, :, :, 1] = self.identity if has_left_interaction: - fac[2, :2, :2, 0] = Operators.sx - fac[3, :2, :2, 0] = Operators.sy + factor[2, :2, :2, 0] = Operators.sx + factor[3, :2, :2, 0] = Operators.sy coeff = self.interaction_matrix[n + 1 :][ None, None, current_right_interactions, n ] - fac[1, :2, :2, 2::2] = coeff * 2 * Operators.sx.unsqueeze(-1) - fac[1, :2, :2, 3::2] = coeff * 2 * Operators.sy.unsqueeze(-1) + factor[1, :2, :2, 2::2] = coeff * 2 * Operators.sx.unsqueeze(-1) + factor[1, :2, :2, 3::2] = coeff * 2 * Operators.sy.unsqueeze(-1) i = 4 if has_left_interaction else 2 j = 2 for current_right_interaction in current_right_interactions.nonzero().flatten(): if right_interactions_to_keep[current_right_interaction]: - fac[i, :, :, j] = self.identity - fac[i + 1, :, :, j + 1] = self.identity + factor[i, :, :, j] = self.identity + factor[i + 1, :, :, j + 1] = self.identity i += 2 j += 2 - return fac + return factor def last_factor(self) -> torch.Tensor: has_left_interaction = self._has_left_interaction(site=-1) left_bond_dim = 4 if has_left_interaction else 2 right_bond_dim = 1 - fac = self._empty_factor(left_bond_dim, right_bond_dim) - fac[0, :, :, 0] = self.identity + factor = self._empty_factor(left_bond_dim, right_bond_dim) + factor[0, :, :, 0] = self.identity if has_left_interaction: coeff = 2 * self.interaction_matrix[0, 1] if self.qubit_count == 2 else 1 - fac[2, :2, :2, 0] = coeff * Operators.sx - fac[3, :2, :2, 0] = coeff * Operators.sy + factor[2, :2, :2, 0] = coeff * Operators.sx + factor[3, :2, :2, 0] = coeff * Operators.sy - return fac + return factor def make_H( From ce97952abfce6e4221d2181551f45dca9dcacc21 Mon Sep 17 00:00:00 2001 From: Kemal Bidzhiev Date: Fri, 27 Mar 2026 11:00:47 +0100 Subject: [PATCH 09/13] clenaup --- emu_mps/hamiltonian.py | 123 +++++++++++++++++++++++++---------------- 1 file changed, 76 insertions(+), 47 deletions(-) diff --git a/emu_mps/hamiltonian.py b/emu_mps/hamiltonian.py index 599d4ca7..e69a269a 100644 --- a/emu_mps/hamiltonian.py +++ b/emu_mps/hamiltonian.py @@ -4,9 +4,8 @@ """ from abc import abstractmethod, ABC -from typing import Iterator, Literal, Union +from typing import Iterator -from pulser.channels.base_channel import States from emu_base import HamiltonianType import torch from emu_mps.mpo import MPO @@ -14,14 +13,10 @@ dtype = torch.complex128 -Eigenstate = Union[States, Literal["0", "1"]] - - class Operators: id = torch.eye(2, dtype=dtype) id_3x3 = torch.eye(3, dtype=dtype) n = torch.tensor([[0.0, 0.0], [0.0, 1.0]], dtype=dtype) - creation = torch.tensor([[0.0, 1.0], [0.0, 0.0]], dtype=dtype) sx = torch.tensor([[0.0, 0.5], [0.5, 0.0]], dtype=dtype) sy = torch.tensor([[0.0, -0.5j], [0.5j, 0.0]], dtype=dtype) @@ -34,17 +29,25 @@ class HamiltonianMPOFactors(ABC): """ def __init__(self, interaction_matrix: torch.Tensor, dim: int = 2): - assert interaction_matrix.ndim == 2, "interaction matrix is not a matrix" - assert ( - interaction_matrix.shape[0] == interaction_matrix.shape[1] - ), "interaction matrix is not square" + self._validate_interaction_matrix(interaction_matrix) + + if dim not in (2, 3): + raise ValueError(f"dim must be 2 or 3, got {dim}") self.dim = dim + self.interaction_matrix = interaction_matrix.clone() self.interaction_matrix.fill_diagonal_(0.0) # or assert - self.qubit_count = self.interaction_matrix.shape[0] - self.middle_site = self.qubit_count // 2 + self.num_sites = self.interaction_matrix.shape[0] + self.middle_site = self.num_sites // 2 self.identity = Operators.id if self.dim == 2 else Operators.id_3x3 + @staticmethod + def _validate_interaction_matrix(matrix: torch.Tensor) -> None: + if matrix.ndim != 2: + raise ValueError("interaction_matrix must be 2-dimensional.") + if matrix.shape[0] != matrix.shape[1]: + raise ValueError("interaction_matrix must be square.") + def __iter__(self) -> Iterator[torch.Tensor]: """Yield the full ordered list of MPO factors for the Hamiltonian.""" yield self.first_factor() @@ -52,10 +55,10 @@ def __iter__(self) -> Iterator[torch.Tensor]: for n in range(1, self.middle_site): yield self.left_factor(n) - if self.qubit_count >= 3: + if self.num_sites >= 3: yield self.middle_factor() - for n in range(self.middle_site + 1, self.qubit_count - 1): + for n in range(self.middle_site + 1, self.num_sites - 1): yield self.right_factor(n) yield self.last_factor() @@ -63,27 +66,22 @@ def __iter__(self) -> Iterator[torch.Tensor]: @abstractmethod def first_factor(self) -> torch.Tensor: """Return the MPO factor for the first site.""" - pass @abstractmethod def left_factor(self, n: int) -> torch.Tensor: """Return the MPO factor for site ``n`` in the left half of the chain.""" - pass @abstractmethod def middle_factor(self) -> torch.Tensor: """Return the MPO factor at the central site bridging both halves.""" - pass @abstractmethod def right_factor(self, n: int) -> torch.Tensor: """Return the MPO factor for site ``n`` in the right half of the chain.""" - pass @abstractmethod def last_factor(self) -> torch.Tensor: """Return the MPO factor for the last site.""" - pass def _has_right_interaction(self, site: int) -> bool: return bool(self.interaction_matrix[site, site + 1 :].any()) @@ -100,14 +98,30 @@ def _empty_factor(self, left_bond_dim: int, right_bond_dim: int) -> torch.Tensor dtype=dtype, ) - def _left_interaction_masks(self, n: int) -> tuple[torch.Tensor, torch.Tensor]: - current_left_interactions = self.interaction_matrix[:n, n:].any(dim=1) - left_interactions_to_keep = self.interaction_matrix[:n, n + 1 :].any(dim=1) + def _left_interaction_masks(self, site: int) -> tuple[torch.Tensor, torch.Tensor]: + """ + For a site in the left half: + - current_left_interactions[i] tells whether site i < site interacts + with current/future sites + - left_interactions_to_keep[i] tells whether that interaction channel + remains active after this site + """ + current_left_interactions = self.interaction_matrix[:site, site:].any(dim=1) + left_interactions_to_keep = self.interaction_matrix[:site, site + 1 :].any(dim=1) return current_left_interactions, left_interactions_to_keep - def _right_interaction_masks(self, n: int) -> tuple[torch.Tensor, torch.Tensor]: - current_right_interactions = self.interaction_matrix[n + 1 :, : n + 1].any(dim=1) - right_interactions_to_keep = self.interaction_matrix[n + 1 :, :n].any(dim=1) + def _right_interaction_masks(self, site: int) -> tuple[torch.Tensor, torch.Tensor]: + """ + For a site in the right half: + - current_right_interactions[j] tells whether site j > site interacts + with current/past sites + - right_interactions_to_keep[j] tells whether that interaction channel + remains active before this site + """ + current_right_interactions = self.interaction_matrix[site + 1 :, : site + 1].any( + dim=1 + ) + right_interactions_to_keep = self.interaction_matrix[site + 1 :, :site].any(dim=1) return current_right_interactions, right_interactions_to_keep def _left_interaction_coefficients( @@ -147,8 +161,9 @@ def first_factor(self) -> torch.Tensor: def left_factor(self, n: int) -> torch.Tensor: has_right_interaction = self._has_right_interaction(site=n) - current_left_interactions = self.interaction_matrix[:n, n:].any(dim=1) - left_interactions_to_keep = self.interaction_matrix[:n, n + 1 :].any(dim=1) + current_left_interactions, left_interactions_to_keep = ( + self._left_interaction_masks(n) + ) left_bond_dim = int(current_left_interactions.sum().item() + 2) right_bond_dim = int( @@ -161,7 +176,7 @@ def left_factor(self, n: int) -> torch.Tensor: if has_right_interaction: factor[1, :2, :2, -1] = Operators.n - coeff = self.interaction_matrix[:n][current_left_interactions, n, None, None] + coeff = self._left_interaction_coefficients(n, current_left_interactions) factor[2:, :2, :2, 0] = coeff * Operators.n i = 2 @@ -175,8 +190,8 @@ def left_factor(self, n: int) -> torch.Tensor: def middle_factor(self) -> torch.Tensor: n = self.middle_site - current_left_interactions = self.interaction_matrix[:n, n:].any(dim=1) - current_right_interactions = self.interaction_matrix[n + 1 :, : n + 1].any(dim=1) + current_left_interactions, _ = self._left_interaction_masks(n) + current_right_interactions, _ = self._right_interaction_masks(n) left_bond_dim = int(current_left_interactions.sum().item() + 2) right_bond_dim = int(current_right_interactions.sum().item() + 2) @@ -185,7 +200,7 @@ def middle_factor(self) -> torch.Tensor: factor[0, :, :, 0] = self.identity factor[1, :, :, 1] = self.identity - coeff = self.interaction_matrix[:n][current_left_interactions, n, None, None] + coeff = self._left_interaction_coefficients(n, current_left_interactions) factor[2:, :2, :2, 0] = coeff * Operators.n coeff = self.interaction_matrix[n + 1 :][ @@ -202,8 +217,9 @@ def middle_factor(self) -> torch.Tensor: def right_factor(self, n: int) -> torch.Tensor: has_left_interaction = self._has_left_interaction(site=n) - current_right_interactions = self.interaction_matrix[n + 1 :, : n + 1].any(dim=1) - right_interactions_to_keep = self.interaction_matrix[n + 1 :, :n].any(dim=1) + current_right_interactions, right_interactions_to_keep = ( + self._right_interaction_masks(site=n) + ) left_bond_dim = int( right_interactions_to_keep.sum().item() + int(has_left_interaction) + 2 @@ -238,7 +254,7 @@ def last_factor(self) -> torch.Tensor: factor = self._empty_factor(left_bond_dim, right_bond_dim) factor[0, :, :, 0] = self.identity if has_left_interaction: - coeff = self.interaction_matrix[0, 1] if self.qubit_count == 2 else 1 + coeff = self.interaction_matrix[0, 1] if self.num_sites == 2 else 1 factor[2, :2, :2, 0] = coeff * Operators.n return factor @@ -260,8 +276,9 @@ def first_factor(self) -> torch.Tensor: def left_factor(self, n: int) -> torch.Tensor: has_right_interaction = self._has_right_interaction(site=n) - current_left_interactions = self.interaction_matrix[:n, n:].any(dim=1) - left_interactions_to_keep = self.interaction_matrix[:n, n + 1 :].any(dim=1) + current_left_interactions, left_interactions_to_keep = ( + self._left_interaction_masks(n) + ) left_bond_dim = int(2 * current_left_interactions.sum().item() + 2) right_bond_dim = int( @@ -277,7 +294,7 @@ def left_factor(self, n: int) -> torch.Tensor: factor[1, :2, :2, -2] = Operators.sx factor[1, :2, :2, -1] = Operators.sy - coeff = self.interaction_matrix[:n][current_left_interactions, n, None, None] + coeff = self._left_interaction_coefficients(n, current_left_interactions) factor[2::2, :2, :2, 0] = coeff * 2 * Operators.sx factor[3::2, :2, :2, 0] = coeff * 2 * Operators.sy @@ -293,8 +310,8 @@ def left_factor(self, n: int) -> torch.Tensor: def middle_factor(self) -> torch.Tensor: n = self.middle_site - current_left_interactions = self.interaction_matrix[:n, n:].any(dim=1) - current_right_interactions = self.interaction_matrix[n + 1 :, : n + 1].any(dim=1) + current_left_interactions, _ = self._left_interaction_masks(n) + current_right_interactions, _ = self._right_interaction_masks(n) left_bond_dim = int(2 * current_left_interactions.sum().item() + 2) right_bond_dim = int(2 * current_right_interactions.sum().item() + 2) @@ -304,7 +321,7 @@ def middle_factor(self) -> torch.Tensor: factor[0, :, :, 0] = self.identity factor[1, :, :, 1] = self.identity - coeff = self.interaction_matrix[:n][current_left_interactions, n, None, None] + coeff = self._left_interaction_coefficients(n, current_left_interactions) factor[2::2, :2, :2, 0] = coeff * 2 * Operators.sx factor[3::2, :2, :2, 0] = coeff * 2 * Operators.sy @@ -324,8 +341,9 @@ def middle_factor(self) -> torch.Tensor: def right_factor(self, n: int) -> torch.Tensor: has_left_interaction = self._has_left_interaction(site=n) - current_right_interactions = self.interaction_matrix[n + 1 :, : n + 1].any(dim=1) - right_interactions_to_keep = self.interaction_matrix[n + 1 :, :n].any(dim=1) + current_right_interactions, right_interactions_to_keep = ( + self._right_interaction_masks(site=n) + ) left_bond_dim = int( 2 * right_interactions_to_keep.sum().item() @@ -365,7 +383,7 @@ def last_factor(self) -> torch.Tensor: factor = self._empty_factor(left_bond_dim, right_bond_dim) factor[0, :, :, 0] = self.identity if has_left_interaction: - coeff = 2 * self.interaction_matrix[0, 1] if self.qubit_count == 2 else 1 + coeff = 2 * self.interaction_matrix[0, 1] if self.num_sites == 2 else 1 factor[2, :2, :2, 0] = coeff * Operators.sx factor[3, :2, :2, 0] = coeff * Operators.sy @@ -383,12 +401,18 @@ def make_H( Constructs and returns a Matrix Product Operator (MPO) representing the neutral atoms Hamiltonian, parameterized by `omega`, `delta`, and `phi`. - The Hamiltonian H is given by: - H = ∑ⱼΩⱼ[cos(ϕⱼ)σˣⱼ + sin(ϕⱼ)σʸⱼ] - ∑ⱼΔⱼnⱼ + ∑ᵢ﹥ⱼC⁶/rᵢⱼ⁶ nᵢnⱼ + The linear terms of the Hamiltonian is + H_0 = ∑ⱼΩⱼ[cos(ϕⱼ)σˣⱼ + sin(ϕⱼ)σʸⱼ] - ∑ⱼΔⱼnⱼ + + The Rydberg Hamiltonian H is given by: + H_R = H_0 + ∑ᵢ﹥ⱼC⁶/rᵢⱼ⁶ nᵢnⱼ + + The XY Hamiltonian H is given by: + H_XY = H_0 + ∑ᵢ﹥ⱼC₃/rᵢⱼ³ 2(SˣᵢSˣⱼ + SʸᵢSʸⱼ) If noise is considered, the Hamiltonian includes an additional term to support the Monte Carlo WaveFunction algorithm: - H = ∑ⱼΩⱼ[cos(ϕⱼ)σˣⱼ + sin(ϕⱼ)σʸⱼ] - ∑ⱼΔⱼnⱼ + ∑ᵢ﹥ⱼC⁶/rᵢⱼ⁶ nᵢnⱼ - 0.5i∑ₘ ∑ᵤ Lₘᵘ⁺ Lₘᵘ + H = H_{R|XY} - 0.5i∑ₘ ∑ᵤ Lₘᵘ⁺ Lₘᵘ where Lₘᵘ are the Lindblad operators representing the noise, m for noise channel and u for the number of atoms @@ -425,6 +449,8 @@ def make_H( num_gpus_to_use=num_gpus_to_use, ) + raise ValueError(f"Unsupported hamiltonian_type: {hamiltonian_type}") + def update_H( hamiltonian: MPO, @@ -455,7 +481,10 @@ def update_H( Defaults to a zero tensor. """ - assert noise.shape in [(2, 2), (3, 3)] + if noise.shape not in {(2, 2), (3, 3)}: + raise ValueError( + f"noise must have shape (2, 2) or (3, 3), got {tuple(noise.shape)}" + ) nqubits = omega.size(dim=0) a = torch.tensordot(omega * torch.cos(phi), Operators.sx, dims=0) From f69facb1b947ae4bb569dc583b39f82876a35556 Mon Sep 17 00:00:00 2001 From: Kemal Bidzhiev Date: Fri, 27 Mar 2026 11:22:42 +0100 Subject: [PATCH 10/13] test --- test/emu_mps/test_hamiltonian.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/test/emu_mps/test_hamiltonian.py b/test/emu_mps/test_hamiltonian.py index 5c3879d3..fdd5c981 100644 --- a/test/emu_mps/test_hamiltonian.py +++ b/test/emu_mps/test_hamiltonian.py @@ -702,3 +702,22 @@ def test_truncation_nn(basis): sv, expected, ) + + +@pytest.mark.parametrize("h_type", (HamiltonianType.Rydberg, HamiltonianType.XY)) +@pytest.mark.parametrize("nqubits", range(3, 12)) +@pytest.mark.parametrize("phys_dim", (2, 3)) +def test_MPO_factor_hermitean(h_type, nqubits, phys_dim): + torch.manual_seed(0) + J = torch.randn(nqubits, nqubits) + ham = make_H( + interaction_matrix=(J + J.T) / 2, + num_gpus_to_use=0, + hamiltonian_type=h_type, + dim=phys_dim, + ) + # Strictly speaking MPO tensor nodes are not required to be Hermitian, + # but the resulting Hamiltonian must be Hermitian. + # Hermiticity of MPO nodes allows to compress baths in TDVP. + for factor in ham.factors: + assert torch.allclose(factor, factor.transpose(1, 2).conj()) From b13f01ce9088c016460423bcb245e3e28f044d02 Mon Sep 17 00:00:00 2001 From: Kemal Bidzhiev Date: Fri, 27 Mar 2026 16:45:25 +0100 Subject: [PATCH 11/13] Mauro comments --- emu_mps/hamiltonian.py | 28 +++++++++++----------------- test/emu_mps/test_hamiltonian.py | 2 +- 2 files changed, 12 insertions(+), 18 deletions(-) diff --git a/emu_mps/hamiltonian.py b/emu_mps/hamiltonian.py index e69a269a..543287e3 100644 --- a/emu_mps/hamiltonian.py +++ b/emu_mps/hamiltonian.py @@ -203,9 +203,7 @@ def middle_factor(self) -> torch.Tensor: coeff = self._left_interaction_coefficients(n, current_left_interactions) factor[2:, :2, :2, 0] = coeff * Operators.n - coeff = self.interaction_matrix[n + 1 :][ - None, None, current_right_interactions, n - ] + coeff = self._right_interaction_coefficients(n, current_right_interactions) factor[1, :2, :2, 2:] = coeff * Operators.n.unsqueeze(-1) coeff = self._middle_interaction_coefficients( @@ -232,9 +230,7 @@ def right_factor(self, n: int) -> torch.Tensor: if has_left_interaction: factor[2, :2, :2, 0] = Operators.n - coeff = self.interaction_matrix[n + 1 :][ - None, None, current_right_interactions, n - ] + coeff = self._right_interaction_coefficients(n, current_right_interactions) factor[1, :2, :2, 2:] = coeff * Operators.n.unsqueeze(-1) i = 3 if has_left_interaction else 2 @@ -325,9 +321,7 @@ def middle_factor(self) -> torch.Tensor: factor[2::2, :2, :2, 0] = coeff * 2 * Operators.sx factor[3::2, :2, :2, 0] = coeff * 2 * Operators.sy - coeff = self.interaction_matrix[n + 1 :][ - None, None, current_right_interactions, n - ] + coeff = self._right_interaction_coefficients(n, current_right_interactions) factor[1, :2, :2, 2::2] = coeff * 2 * Operators.sx.unsqueeze(-1) factor[1, :2, :2, 3::2] = coeff * 2 * Operators.sy.unsqueeze(-1) @@ -359,9 +353,7 @@ def right_factor(self, n: int) -> torch.Tensor: factor[2, :2, :2, 0] = Operators.sx factor[3, :2, :2, 0] = Operators.sy - coeff = self.interaction_matrix[n + 1 :][ - None, None, current_right_interactions, n - ] + coeff = self._right_interaction_coefficients(n, current_right_interactions) factor[1, :2, :2, 2::2] = coeff * 2 * Operators.sx.unsqueeze(-1) factor[1, :2, :2, 3::2] = coeff * 2 * Operators.sy.unsqueeze(-1) @@ -401,18 +393,20 @@ def make_H( Constructs and returns a Matrix Product Operator (MPO) representing the neutral atoms Hamiltonian, parameterized by `omega`, `delta`, and `phi`. - The linear terms of the Hamiltonian is + The linear term of the Hamiltonian is H_0 = ∑ⱼΩⱼ[cos(ϕⱼ)σˣⱼ + sin(ϕⱼ)σʸⱼ] - ∑ⱼΔⱼnⱼ - The Rydberg Hamiltonian H is given by: - H_R = H_0 + ∑ᵢ﹥ⱼC⁶/rᵢⱼ⁶ nᵢnⱼ + The Rydberg Hamiltonian is given by: + H_Ryd = H_0 + ∑ᵢ﹥ⱼC⁶/rᵢⱼ⁶ nᵢnⱼ - The XY Hamiltonian H is given by: + The XY Hamiltonian is given by: H_XY = H_0 + ∑ᵢ﹥ⱼC₃/rᵢⱼ³ 2(SˣᵢSˣⱼ + SʸᵢSʸⱼ) If noise is considered, the Hamiltonian includes an additional term to support the Monte Carlo WaveFunction algorithm: - H = H_{R|XY} - 0.5i∑ₘ ∑ᵤ Lₘᵘ⁺ Lₘᵘ + H = H_Ryd - 0.5i∑ₘ ∑ᵤ Lₘᵘ⁺ Lₘᵘ + or + H = H_XY - 0.5i∑ₘ ∑ᵤ Lₘᵘ⁺ Lₘᵘ where Lₘᵘ are the Lindblad operators representing the noise, m for noise channel and u for the number of atoms diff --git a/test/emu_mps/test_hamiltonian.py b/test/emu_mps/test_hamiltonian.py index fdd5c981..3cb0180c 100644 --- a/test/emu_mps/test_hamiltonian.py +++ b/test/emu_mps/test_hamiltonian.py @@ -705,7 +705,7 @@ def test_truncation_nn(basis): @pytest.mark.parametrize("h_type", (HamiltonianType.Rydberg, HamiltonianType.XY)) -@pytest.mark.parametrize("nqubits", range(3, 12)) +@pytest.mark.parametrize("nqubits", range(2, 12)) @pytest.mark.parametrize("phys_dim", (2, 3)) def test_MPO_factor_hermitean(h_type, nqubits, phys_dim): torch.manual_seed(0) From 1742fa98eafeb472cfb89653732e0c5c5917e238 Mon Sep 17 00:00:00 2001 From: Kemal Bidzhiev Date: Mon, 30 Mar 2026 13:49:38 +0200 Subject: [PATCH 12/13] xy docstring --- emu_mps/hamiltonian.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/emu_mps/hamiltonian.py b/emu_mps/hamiltonian.py index 543287e3..3822542d 100644 --- a/emu_mps/hamiltonian.py +++ b/emu_mps/hamiltonian.py @@ -257,6 +257,17 @@ def last_factor(self) -> torch.Tensor: class XYHamiltonianMPOFactors(HamiltonianMPOFactors): + """ + Note: + The XY Hamiltonian is implemented using the Pauli matrices X and Y, + which ensures that the MPO nodes are Hermitian. This representation + is later used to compress energy baths in TDVP and DMRG. + + An alternative approach is to use the sigma^+ and sigma^- operators. + In that case, restoring Hermiticity of the MPO nodes requires adding + the Hermitian conjugate via a direct sum. + """ + def first_factor(self) -> torch.Tensor: has_right_interaction = self._has_right_interaction(site=0) From 20712dea7c4596d359142227fd6a2e5d4bf66ac4 Mon Sep 17 00:00:00 2001 From: Kemal Bidzhiev Date: Mon, 30 Mar 2026 15:35:25 +0200 Subject: [PATCH 13/13] next comment --- emu_mps/hamiltonian.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/emu_mps/hamiltonian.py b/emu_mps/hamiltonian.py index 3822542d..99d267ee 100644 --- a/emu_mps/hamiltonian.py +++ b/emu_mps/hamiltonian.py @@ -69,7 +69,8 @@ def first_factor(self) -> torch.Tensor: @abstractmethod def left_factor(self, n: int) -> torch.Tensor: - """Return the MPO factor for site ``n`` in the left half of the chain.""" + """Return the MPO factor for site ``n`` in the left half of the chain + except the first factor.""" @abstractmethod def middle_factor(self) -> torch.Tensor: @@ -77,7 +78,8 @@ def middle_factor(self) -> torch.Tensor: @abstractmethod def right_factor(self, n: int) -> torch.Tensor: - """Return the MPO factor for site ``n`` in the right half of the chain.""" + """Return the MPO factor for site ``n`` in the right half of the chain + except the last factor.""" @abstractmethod def last_factor(self) -> torch.Tensor: @@ -102,7 +104,7 @@ def _left_interaction_masks(self, site: int) -> tuple[torch.Tensor, torch.Tensor """ For a site in the left half: - current_left_interactions[i] tells whether site i < site interacts - with current/future sites + with current/next sites - left_interactions_to_keep[i] tells whether that interaction channel remains active after this site """ @@ -114,7 +116,7 @@ def _right_interaction_masks(self, site: int) -> tuple[torch.Tensor, torch.Tenso """ For a site in the right half: - current_right_interactions[j] tells whether site j > site interacts - with current/past sites + with current/previous sites - right_interactions_to_keep[j] tells whether that interaction channel remains active before this site """ @@ -269,6 +271,7 @@ class XYHamiltonianMPOFactors(HamiltonianMPOFactors): """ def first_factor(self) -> torch.Tensor: + print("SxSy Hamiltonian") has_right_interaction = self._has_right_interaction(site=0) left_bond_dim = 1