Source code for oqupy.mps_mpo

# Copyright 2022 The TEMPO Collaboration
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# See the License for the specific language governing permissions and
# limitations under the License.
Module for MPSs and MPOs.

from copy import deepcopy
from typing import Any, List, Optional, Text, Tuple

import numpy as np
from numpy import ndarray
from scipy import linalg

from oqupy.base_api import BaseAPIClass
from oqupy.config import NpDtype
SystemChain = Any # from oqupy.system import SystemChain

[docs]class Gate: """ Class representing an n-site gate in MPO form. The axes of the MPO tensors are (L = left bond leg, I = input leg, O = output leg, R = right bond leg): * for n=1: (I, O) * for n=2: (I, O, R), (L, I, O) * for n>2: (I, O, R), (L, I, O, R), ..., (L, I, O, R), (L, I, O) Parameters ---------- sites: List[int] The site numbers onto which the MPO Gate acts. tensors: List[ndarray] The MPO tensors of the gate. """ def __init__( self, sites: List[int], tensors: List[ndarray]) -> None: """Create a n-site gate in MPO form. """ for site in sites: assert isinstance(site, int) for site_l, site_r in zip(sites[:-1], sites[1:]): assert site_l < site_r self._sites = sites assert len(sites) == len(tensors) self._length = len(sites) self._span = sites[-1] - sites[0] tmp_tensors = [] for tensor in tensors: tmp_tensor = np.array(tensor, dtype=NpDtype) tmp_tensors.append(tmp_tensor) for i, tensor in enumerate(tmp_tensors): number_of_legs = 4 number_of_legs -= 1 if i == 0 else 0 number_of_legs -= 1 if i == self._length-1 else 0 assert len(tensor.shape) == number_of_legs for tensor_l, tensor_r in zip(tmp_tensors[:-1], tmp_tensors[1:]): assert tensor_l.shape[-1] == tensor_r.shape[0] self._tensors = tmp_tensors def __len__(self): """Number of sites onto which the gate acts. """ return self._length @property def span(self) -> int: """The span of sites onto which the gate acts. """ return self._span @property def sites(self) -> List[int]: """The sites onto which the gate acts. """ return self._sites @property def tensors(self) -> List[ndarray]: """The tensors of the MPO gate. """ return self._tensors
[docs]class SiteGate(Gate): """ An MPO gate acting on a single site. Parameters ---------- site: int The site onto which the MPO gate acts. tensor: ndarray The single site MPO (which is a matrix). """ def __init__( self, site: int, tensor: ndarray) -> None: """Create a single site MPO gate. """ super().__init__([site], [tensor])
[docs]class NnGate(Gate): """ An MPO gate acting on a pair of neighboring sites. Parameters ---------- site: int The index of the left site. tensor: ndarray The two MPO tensors of shape. """ def __init__( self, site: int, tensors: Tuple[(ndarray, ndarray)]) -> None: """Create a nearest neighbor gate. """ super().__init__([site, site+1], [tensors[0], tensors[1]])
[docs]class GateLayer: """ A layer of gates. Parameters ---------- parallel: bool Whether of not the gates are suitable for a parallel application. gates: List[Gates] List of gates. """ def __init__( self, parallel: bool, gates: List[Gate]) -> None: """Create a GateLayer object. """ assert isinstance(parallel, bool) assert isinstance(gates, list) for gate in gates: assert isinstance(gate, Gate) self._parallel = parallel self._gates = gates @property def parallel(self) -> bool: """Whether of not the gates are suitable for a parallel application. """ return self._parallel @property def gates(self) -> List[Gate]: """List of gates. """ return self._gates
[docs]def compute_nn_gate( liouvillian: ndarray, site: int, hs_dim_l: int, hs_dim_r: int, dt: float, epsrel: float) -> NnGate: """ Compute nearest neighbor gate from Liouvillian. Parameters ---------- liouvillian: ndarray The two site Liouvillian. site: int The index of the left site. hs_dim_l: int The Hilbert space dimension of the left site. hs_dim_r: int The Hilbert space dimension of the right site. dt: float The time step length. epsrel: float The relative singular value truncation tolerance. Returns ------- nn_gate: NnGate Nearest neighbor gate. """ # exponentiate the liouvillian to become a propagator propagator = linalg.expm(dt*liouvillian) # split leg 0 and leg 1 each into left and right. propagator.shape = [hs_dim_l**2, # left output hs_dim_r**2, # right output hs_dim_l**2, # left input hs_dim_r**2] # right input temp = np.swapaxes(propagator, 1, 2) temp = temp.reshape([hs_dim_l**2 * hs_dim_l**2, hs_dim_r**2 * hs_dim_r**2]) u_full, s_full, vh_full = linalg.svd(temp, full_matrices=False) chi = _truncation_index(s_full, epsrel) s = s_full[:chi] u=u_full[:,:chi] vh=vh_full[:chi,:] sqrt_s = np.sqrt(s) u_sqrt_s = u * sqrt_s sqrt_s_vh =(sqrt_s * vh.T).T # left tensor with legs: left output, left input, bond tensor_l = u_sqrt_s.reshape(hs_dim_l**2, hs_dim_l**2, chi) # right tensor with legs: bond, right output, right input tensor_r = sqrt_s_vh.reshape(chi, hs_dim_r**2, hs_dim_r**2) return NnGate(site=site, tensors=(tensor_l, tensor_r))
[docs]def compute_trotter_layers( nn_full_liouvillians: List[ndarray], hs_dims: List[int], dt: float, epsrel: float) -> Tuple[GateLayer, GateLayer]: """ Compute even and odd Trotter layers. Parameters ---------- nn_full_liouvillians: List[ndarrays] Full list of nearest neighbor Liouvillians. hs_dims: List[int] Hilbert space dimensions of the chain sites. dt: float The time step length. epsrel: float The relative singular value truncation tolerance. Returns ------- gate_layer_even: GateLayer Gate layer with nearest neighbor gates with left sites having even indices. gate_layer_odd: GateLayer Gate layer with nearest neighbor gates with left sites having odd indices. """ all_gates = [] for i, liouv in enumerate(nn_full_liouvillians): gate = compute_nn_gate(liouvillian = liouv, site=i, hs_dim_l=hs_dims[i], hs_dim_r=hs_dims[i+1], dt=dt, epsrel=epsrel) all_gates.append(gate) gates_even = all_gates[0::2] gates_odd = all_gates[1::2] gate_layer_even = GateLayer(parallel=True, gates=gates_even) gate_layer_odd = GateLayer(parallel=True, gates=gates_odd) return [gate_layer_even, gate_layer_odd]
[docs]class TebdPropagator: """ TEBD (Time Evolving Block Decimation) Propagators consist of a list of GateLayers. Parameters ---------- gate_layers: List[GateLayer] The gate layers that make up a full time step propagation in a TEBD tensor network. """ def __init__( self, gate_layers: List[GateLayer]) -> None: """Create a TebdPropagators object. """ self._gate_layers = gate_layers @property def gate_layers(self) -> List[GateLayer]: """ The gate layers that make up a full time step propagation in a TEBD tensor network. """ return self._gate_layers
[docs]def compute_tebd_propagator( system_chain: SystemChain, time_step: float, epsrel: float, order: int) -> TebdPropagator: """ Compute a TebdPropagator object for a given SystemChain. Parameters ---------- system_chain: SystemChain A SystemChain object that encodes the nearest neighbor Liouvillians. time_step: float The time step length of the full TEBD propagator. epsrel: float The relative singular value truncation tolerance. order: int The expansion order. Returns ------- tebd_propagator: TebdPropagator The TEBD Propagator. """ nn_full_liouvillians = system_chain.get_nn_full_liouvillians() hs_dims = system_chain.hs_dims if order == 1: layers = compute_trotter_layers( nn_full_liouvillians=nn_full_liouvillians, hs_dims=hs_dims, dt=time_step, epsrel=epsrel) propagator = TebdPropagator(gate_layers=[layers[0], layers[1]]) elif order == 2: layers = compute_trotter_layers( nn_full_liouvillians=nn_full_liouvillians, hs_dims=hs_dims, dt=time_step/2.0, epsrel=epsrel) propagator = TebdPropagator(gate_layers=[layers[0], layers[1], layers[1], layers[0]]) else: raise NotImplementedError(f"Trotter layers of order {order} are " \ + "not implemented.") return propagator
[docs]class AugmentedMPS(BaseAPIClass): """ An augmented matrix product state (as introduced in the supplemental material of [Fux2022]). The full gamma tensors (one for each site) have the following axis: (L = left bond leg, P = physical leg, R = right bond leg, A = augmented leg). Depending on the rank of the input gamma tensor, it is completed with legs of dimension 1 according to the following interpretation of the input: * rank = 1: Product state vectorized density matrix, i.e. (1, P, 1, 1) * rank = 2: Product state density matrix, i.e. (1, p*p, 1, 1) * rank = 3: Canonical MPS gamma tensor, i.e. (L, P, R, 1) * rank = 4: Augmented MPS gamma tensor, i.e. (L, P, R, A) If no lambdas are given, they are assumed to be identities. A single lambda can be `None` (identity), a vector (giving the diagonal values) or a matrix (which must be diagonal). Parameters ---------- gammas: List[ndarray] The input gamma tensors. lambdas: List[ndarray] The input lambda diagonal matrices. name: str An optional name for the augmented MPS. description: str An optional description of the augmented MPS. """ def __init__( self, gammas: List[ndarray], lambdas: Optional[List[ndarray]] = None, name: Optional[Text] = None, description: Optional[Text] = None) -> None: """Create an AugmentedMPS object. """ # input parsing self._n = len(gammas) assert self._n > 0 if lambdas is not None: assert len(lambdas) == len(gammas)-1 tmp_gammas = [] for g in gammas: tmp_gamma = np.array(g, dtype=NpDtype) shape = deepcopy(tmp_gamma.shape) rank = len(shape) if rank == 4: tmp_gamma = np.swapaxes(tmp_gamma,2,3) elif rank == 3: tmp_gamma.shape = (shape[0], shape[1], 1, shape[2]) elif rank == 2: tmp_gamma.shape = (1, shape[0]*shape[1], 1, 1) elif rank == 1: tmp_gamma.shape = (1, shape[0], 1, 1) else: raise ValueError() tmp_gammas.append(tmp_gamma) bond_dims = [] for g1, g2 in zip(tmp_gammas[:-1], tmp_gammas[1:]): assert g1.shape[3] == g2.shape[0] bond_dims.append(g1.shape[3]) if lambdas is None: lambdas = [None] * (self._n - 1) tmp_lambdas = [] for bond_dim, l in zip(bond_dims, lambdas): if l is None: tmp_lambda = np.ones(bond_dim, dtype=NpDtype) else: tmp_lambda = np.array(l, dtype=NpDtype) shape = tmp_lambda.shape rank = len(shape) if rank == 2: assert np.all( tmp_lambda == np.diag(np.diagonal(tmp_lambda))) tmp_lambda = np.diagonal(tmp_lambda) elif rank == 1: pass else: raise ValueError() assert np.all(tmp_lambda > 0.0), \ "All lambda matrix diagonal values must be positive. " assert len(tmp_lambda) == bond_dim tmp_lambdas.append(tmp_lambda) self._gammas = tmp_gammas self._lambdas = tmp_lambdas super().__init__(name, description) @property def gammas(self) -> ndarray: """"The gamma tensors.""" return self._gammas @property def lambdas(self) -> ndarray: """"The values of the lambda matrices diagonals.""" return self._lambdas
def _truncation_index(s: ndarray, epsrel: float) -> int: """Helper function to figure out the right singular value cutoff. """ absrel = s[0] * epsrel cummulative_square = np.cumsum(np.flip(s)**2) chi = np.count_nonzero(cummulative_square > absrel**2) return chi