Operator module

The operator module provides functionality to evaluate the action of operators on computational basis configurations. In terms of an operator acting on the pure state Hilbert space this corresponds to ‘’on-the-fly’’ generation of matrix elements. Generally, this amounts to a mapping

\[s \rightarrow \{s_j'\}, \{O_{s,s_j'}\}\]

Here, \(s\) denotes the input basis configuration and the \(s'\) are the connected basis configurations produced by the operator. The \(O_{s,s'}\) are the coefficients associated with the map \(s\to s'\).

Abstract operator class

The abstract operator class Operator defines the general interface of operators acting on computational basis configurations. Any specific implementation should be a child of Operator and implement a compile function.

class jVMC.operator.Operator(ElocBatchSize=-1)

This class defines an interface and provides functionality to compute operator matrix elements

This is an abstract class. It provides the general interface to compute operator matrix elements, but it lacks the implementation of operator definitions. Arbitrary operators can be constructed as classes that inherit from Operator and implement the compile() method.

The compile() method has to return a jit-able pure function with the following interface:

Arguments:

s: A single basis configuration.

*args: Further positional arguments. E.g. time in the case of time-dependent operators.

Returns:
A tuple sp, matEls, where sp is the list of connected basis configurations (as jax.numpy.array) and matEls the corresponding matrix elements.

Alternatively, compile() can return a tuple of two functions, the first as described above and and the second a preprocessor for the additional positional arguments *args. Assuming that compile() returns the tuple (f1, f2), f1 will be called as f1(s, f2(*args)) .

Important: Any child class inheriting from Operator has to call super().__init__() in its constructor.

Example:

Here we define a \(\hat\sigma_1^x\) operator acting on lattice site \(1\) of a spin-1/2 chain.

    def __init__(self, siteIdx):

        self.siteIdx = siteIdx

        super().__init__()  # Constructor of base class Operator has to be called!

    def compile(self):

        def get_s_primes(s, *args, idx):

            sp = s.copy()

            configShape = sp.shape
            sp = sp.ravel()

            # Define matrix element
            matEl = jnp.array([1., ], dtype=global_defs.tCpx)
            # Define mapping of Sx: 0->1, 1->0
            sMap = jnp.array([1, 0])
            # Perform mapping
            sp = jax.ops.index_update(sp, jax.ops.index[idx], sMap[s[idx]])

            return sp.reshape(configShape), matEl

        # Create a pure function that takes only a basis configuration as argument
        map_function = functools.partial(get_s_primes, idx=self.siteIdx)

        return map_function


mySx = SxOperator(siteIdx=1)


# Get operator mapping function
testFunction = mySx.compile()

# Simple test configuration s=|0000>
testConfig = jnp.array([0, 0, 0, 0], dtype=np.int32)

# Get connected s' and matrix elements <s'|Sx|s>
sp, matEl = testFunction(testConfig)

Then we indeed obtain the correct sp and matEl:

>>> print(sp)
[0 1 0 0]
>>> print(matEl)
[1.+0.j]

abstract compile(): An implementation of this method should return a jit-able function that returns for a given basis configuration s the connected configurations s' and the corresponding matrix elements.

get_O_loc(samples, psi, logPsiS=None, *args)

Compute \(O_{loc}(s)\).

If the instance parameter ElocBatchSize is larger than 0 \(O_{loc}(s)\) is computed in a batch-wise manner to avoid out-of-memory issues.

Arguments:

samples: Sample of computational basis configurations \(s\).
psi: Neural quantum state.
logPsiS: Logarithmic amplitudes \(\ln(\psi(s))\)
*args: Further positional arguments for the operator.

Returns:

\(O_{loc}(s)\) for each configuration \(s\).

get_O_loc_batched(samples, psi, logPsiS, batchSize, *args)

Compute \(O_{loc}(s)\) in batches.

Computes \(O_{loc}(s)=\sum_{s'} O_{s,s'}\frac{\psi(s')}{\psi(s)}\) in a batch-wise manner to avoid out-of-memory issues.

Arguments:

samples: Sample of computational basis configurations \(s\).
psi: Neural quantum state.
logPsiS: Logarithmic amplitudes \(\ln(\psi(s))\)
batchSize: Batch size.
*args: Further positional arguments for the operator.

Returns:

\(O_{loc}(s)\) for each configuration \(s\).

get_O_loc_unbatched(logPsiS, logPsiSP)

Compute \(O_{loc}(s)\).

This member function assumes that get_s_primes(s) has been called before, as internally stored matrix elements \(O_{s,s'}\) are used.

Computes \(O_{loc}(s)=\sum_{s'} O_{s,s'}\frac{\psi(s')}{\psi(s)}\), given the logarithmic wave function amplitudes of the involved configurations \(\ln(\psi(s))\) and \(\ln\psi(s')\)

Arguments:

logPsiS: Logarithmic amplitudes \(\ln(\psi(s))\)
logPsiSP: Logarithmic amplitudes \(\ln(\psi(s'))\)

Returns:

\(O_{loc}(s)\) for each configuration \(s\).

get_estimator_function(psi, *args)

Get a function that computes \(O_{loc}(\theta, s)\).

Returns a function that computes \(O_{loc}(\theta, s)=\sum_{s'} O_{s,s'}\frac{\psi_\theta(s')}{\psi_\theta(s)}\) for a given configuration \(s\) and parameters \(\theta\) of a given ansatz \(\psi_\theta(s)\).

Arguments:

psi: Neural quantum state.
*args: Further positional arguments for the operator.

Returns:

A function \(O_{loc}(\theta, s)\).

get_s_primes(s, *args)

Compute matrix elements

For a list of computational basis states \(s\) this member function computes the corresponding matrix elements \(O_{s,s'}=\langle s|\hat O|s'\rangle\) and their respective configurations \(s'\).

Arguments:

s: Array of computational basis states.
*args: Further positional arguments that are passed on to the specific operator implementation.

Returns:

An array holding all configurations \(s'\) and the corresponding matrix elements \(O_{s,s'}\).

Branch-free operator class

These operators are intended to act on the computational basis of a many-body Hilbert space. The many-body Hilbert space is a product of local Hilbert spaces, \(\mathscr H=\bigotimes_l \mathscr h_l\). Operators generally take the form \(\hat O=\sum_k \hat O_k\), where \(\hat O_k=\prod_l \hat o_l^k\) are operator strings made up of elementary operators \(\hat o_l\) acting on the factors \(\mathscr h_l\).

Note

The key assumption of the BranchFreeOperator is that the elementary operators \(\hat o_l\) are ‘’branch-free’’, meaning that the correspondig matrices have only one non-zero entry per row.

In variational Monte Carlo a key quantity is

\(O_{loc}(s) = \sum_{s'}O_{s,s'}\psi(s')/\psi(s)\)

where \(s,s'\) denote computational basis states.

This module provides functionality to assemble general operators from elementary operators, to compute matrix elements \(O_{s,s'}\), and finally to compute \(O_{loc}(s)\).

Elementary operators

At the core the operator class works with a general description of elementary operators \(\hat o_l\). Elementary operators are defined by dictionaries that hold four items:

'idx': Index of the corresponding local Hilbert space.

'map': An array (jax.numpy.array) defining the mapping between local basis indices. The array entry \(j\) holds the image of the basis index \(j\) under the map.

'matEls': An array (jax.numpy.array) defining the corresponding matrix elements.

'diag': Boolean stating whether this operator is a diagonal operator (exploiting this information enhances efficiency).

For concreteness, consider the Pauli operator

\(\hat\sigma^x=\begin{pmatrix}0&1\\1&0\end{pmatrix}\)

acting on lattice site \(l=1\). The corresponding dictionary is:

Sx = {
        'idx': 1,
        'map': jax.numpy.array([1,0],dtype=np.int32),
        'matEl': jax.numpy.array([1.,1.],dtype=jVMC.global_defs.tReal),
        'diag': False
     }

A number of frequently used operators is pre-defined in this module, see below.

Operator strings

Operator strings are treated as tuples of elementary operators. For example, using the pre-defined Pauli operators Sx and Sz, the operator string \(\hat\sigma_1^x\hat\sigma_2^z\) is obtained as:

X1Z2 = ( Sx(1), Sz(2) )

Prefactors can be added to operator strings using the scal_opstr() function. For example, to obtain \(\frac{1}{2}\hat\sigma_1^x\hat\sigma_2^z\):

X1Z2_with_prefactor = scal_opstr(0.5, X1Z2)

Alternatively, a function that depends on some arguments, f(*args), can be given as prefactor:

X1Z2_with_prefactor = scal_opstr(f, X1Z2)

In that case, the respective arguments have to be passed as additional arguments, whenever the get_s_primes function is called.

Assembling operators

Finally, arbitrary operators can be assembled from operator strings. Consider, e.g., the Hamiltonian of the spin-1/2 quantum Ising chain of length L,

\(\hat H=-\sum_{l=0}^{L-2}\hat\sigma_l^z\hat\sigma_{l+1}^z-g\sum_{l=0}^{L-1}\hat\sigma_l^x\)

Again, using the pre-defined Pauli operators Sx and Sz, an Operator object for this Hamiltonian can be obtained as:

hamiltonian = Operator()
for l in range(L-1):
    hamiltonian.add( scal_opstr( -1., ( Sz(l), Sz(l+1) ) ) )
    hamiltonian.add( scal_opstr( -g, ( Sx(l), ) ) )
hamiltonian.add( scal_opstr( -g, ( Sx(L-1), ) ) )

Detailed documentation

class jVMC.operator.BranchFreeOperator(lDim=2, **kwargs)

This class provides functionality to compute operator matrix elements

Initializer arguments:

lDim: Dimension of local Hilbert space.

add(opDescr)

Add another operator to the operator

Arguments:

opDescr: Operator string to be added to the operator.

compile(): Compiles a operator mapping function from the given operator strings.

jVMC.operator.Id(idx=0, lDim=2)

Returns an identity operator

Args:

idx: Index of the local Hilbert space.
lDim: Dimension of local Hilbert space.

Returns:: Dictionary defining an identity operator

jVMC.operator.Sx(idx)

Returns a \(\hat\sigma^x\) Pauli operator

Args:

idx: Index of the local Hilbert space.

Returns:: Dictionary defining \(\hat\sigma^x\) Pauli operator

jVMC.operator.Sy(idx)

Returns a \(\hat\sigma^x\) Pauli operator

Args:

idx: Index of the local Hilbert space.

Returns:: Dictionary defining \(\hat\sigma^x\) Pauli operator

jVMC.operator.Sz(idx)

Returns a \(\hat\sigma^z\) Pauli operator

Args:

idx: Index of the local Hilbert space.

Returns:: Dictionary defining \(\hat\sigma^z\) Pauli operator

jVMC.operator.Sp(idx)

Returns a \(S^+\) spin-1/2 ladder operator

Args:

idx: Index of the local Hilbert space.

Returns:: Dictionary defining \(S^+\) ladder operator

jVMC.operator.Sm(idx)

Returns a \(S^-\) spin-1/2 ladder operator

Args:

idx: Index of the local Hilbert space.

Returns:: Dictionary defining \(S^-\) ladder operator

jVMC.operator.scal_opstr(a, op)

Add prefactor to operator string

Arguments:

a: Scalar prefactor or function.
op: Operator string.

Returns:

Operator string rescaled by a.

POVM operator class

The POVM-formalism builds on the observation that a density matrix \(\rho\) may equivalently (without information loss) be represented as a probability distribution \(P\) with \(4^N\) entries. Elements of this distribution are obtained as expectation values of the density matrix \(P^\textbf{a}=\mathrm{tr}(\rho M^\textbf{a})\) for the product POVM-operators \(M^\textbf{a}=M^{a_1}\otimes M^{a_2}\otimes..\otimes M^{a_N}\). In order to simulate dissipative dynamics in Open Quantum Systems (OQS) the Lindbladian evolution is expressed as \(\dot{P}^\textbf{a}=\mathcal{L}^\textbf{ab}P^\textbf{b}\), where \(\mathcal{L}\) is the operator that generates the evolution. The POVM operator class comprises an implementation of \(\mathcal{L}^\textbf{ab}\), whose main purpose it is to output all off-diagonal configurations \(\textbf{b}=b_1b_2..b_N\) associated with an input configuration \(\textbf{a}=a_1a_2..a_N\). This is achieved by building on the sparsity that is typically present in many-body systems, as this sparsity crucially is preserved in the POVM-formalism allowing the described operations to scale linearly in the system size \(N\).

Elementary operators

Elementary operators in the POVM-formalism are build by reexpressing the equation of motion of the density operator \(\dot{\rho}=\mathcal{L}[\rho]\) in terms of \(P\), resulting in

\(\dot{P}^\textbf{a} = \mathrm{tr}(\mathcal{L}[M^\textbf{c}]M^\textbf{a}) T^{-1\textbf{cb}} P^\textbf{b} = \mathcal{L}^\textbf{ab}P^\textbf{b}\).

Typically \(\mathcal{L}[\rho]\) consists of 2-body (1-body) operators which translate to (real) \(16\times 16\) (\(4\times 4\)) matrices in the POVM-picture. Importantly, only objects of this size need to be stored, but \(n\)-body operators for \(n>2\) are also supported. Frequently encountered unitary and dissipative operators are pre-defined and can be constructed as explained below.

Assembling operators

Consider again the above example of a spin-1/2 quantum Ising chain of length \(L\) (here with periodic boundary conditions),

\(\hat H=-\sum_{l=1}^{L}\hat\sigma_l^z\hat\sigma_{(l+1) \% L}^z-g\sum_{l=1}^{L}\hat\sigma_l^x\).

However, let us now additionally assume a single-qubit decay on each spin in the chain, such that the time-derivative of the density matrix is

\(\dot{\rho} = -i[H, \rho] + \gamma \sum_i \sigma^-_i\rho\sigma^+_i - \frac{1}{2}\lbrace \sigma^+_i\sigma^-_i,\rho\rbrace\)

Using the POVM-operator class, the expression for the operator that corresponds to these dynamics reads:

Lindbladian = jVMC.operator.POVMOperator()
for l in range(L):
    Lindbladian.add({"name": "ZZ", "strength": -1.0, "sites": (l, (l + 1) % L)})
    Lindbladian.add({"name": "X", "strength": -g, "sites": (l,)})
    Lindbladian.add({"name": "decaydown", "strength": gamma, "sites": (l,)})

Adding terms to the Lindbladian is done using dictionaries, which have three entries: The name of the operator to be added, its prefactor and the site-ids that are involved. Valid names that are recognized by default are the unitary 2-body (1-body) operators ["XX", "YY", "ZZ"] (["X", "Y", "Z"]) corresponding to couplings and external magnetic fields aswell as the single-particle dissipation terms ["dephasing", "decaydown", "decayup"] corresponding to the dissipation operators [\(\sigma^z\), \(\sigma^-\), \(\sigma^+\)].

Adding new operators

To add custom operators it is best to use the matrix_to_povm function and the appropiate adding method of the POVM object. For example adding the operator \(A = \sigma_i^x \sigma_j^z\) as part of the unitary term can be done by:

povm = POVM({"dim": "1D", "L": L})
M, T_inv = povm.M, povm.T_inv
M_2body = jax.numpy.array([[jax.numpy.kron(M[i], M[j]) for j in range(4)]
                            for i in range(4)]).reshape(4**2, 2**2, 2**2)
T_inv_2body = jax.numpy.kron(T_inv, T_inv)

omega = matrix_to_povm(A, M_2body, T_inv_2body, mode='unitary')
povm.add_unitary("XZ", omega)

Valid values for the mode parameter are unitary (uni), dissipative (dis), observable (obs) and imaginary (imag). The name for a operator added to the POVM object can not be used twice even for different modes, trying to add an operator with a already used name will raise a ValueError.

Detailed documentation

class jVMC.operator.POVM(system_data, theta=0, phi=0, name='SIC')

This class provides POVM - operators and related matrices.

Initializer arguments:

system_data: dictionary with lattice dimension "dim" (either "1D" or "2D") and length "L" specifying the computation of correlators.
theta: angle theta on the Bloch - sphere
phi: angle phi on the Bloch - sphere
name: specifier of the POVM

evaluate_observable(operator, states): Obtain X, Y and Z magnetizations and their correlators up to the specified length in system data["L"]. Note that this function assumes translational invariance and averages the results obtained for single spins.

set_standard_povm_operators(): Obtain matrices required for dynamics and observables.

class jVMC.operator.POVMOperator(povm, ldim=4, **kwargs)

This class provides functionality to compute operator matrix elements.

Initializer arguments:

povm: An instance of the POVM-class.

lDim: Dimension of local Hilbert space.

add(opDescr)

Add another operator to the operator.

Args:

opDescr: Operator dictionary to be added to the operator.

compile(): Compiles an operator mapping function from the previously added dictionaries.

jVMC.operator.measure_povm(povm, sampler, sampleConfigs=None, probs=None, observables=None)

For a set of sampled configurations, compute the associated expectation values for a given set of observables. If none is provided, magnetizations and correlations for X, Y and Z are computed.

Args:

povm: the povm that holds the jitted evaluation function
sampler: the sampler used for the computation of expectation values
sampleConfigs: optional, if configurations have already been generated
probs: optional, the associated probabilities
observables: optional, observables for which expectation values are computed

jVMC.operator.get_M(theta, phi, name)

Returns 4 POVM measurement operators.

Args:

theta: angle theta on the Bloch - sphere
phi: angle phi on the Bloch - sphere
name: specifier of the POVM

Returns:

jnp.array with the leading axis giving the different POVM-Measurement operators.

jVMC.operator.get_dissipators(M, T_inv)

Get the dissipation operators in the POVM-formalism.

Args:

M: POVM-Measurement operators
T_inv: Inverse POVM-Overlap matrix

Returns:

Dictionary with common (unscaled) 1-body dissipation channels.

jVMC.operator.get_unitaries(M, T_inv)

Get common 1- and 2-body unitary operators in the POVM formalism.

Args:

M: POVM-Measurement operators
T_inv: Inverse POVM-Overlap matrix

Returns:

Dictionary with common 1- and 2-body unitary interactions.

jVMC.operator.get_observables(M, T_inv)

Get X, Y and Z observables in the POVM-formalism.

Args:

M: POVM-Measurement operators
T_inv: Inverse POVM-Overlap matrix

Returns:

Dictionary giving the X, Y and Z magnetization.

jVMC.operator.get_1_particle_distributions(state, povm)

compute 1 particle POVM-representations, mainly used for defining initial states.

Args:

state: the desired state on the bloch-sphere, set together from ["x", "y", "z"] + "_" + ["up", "down"]
povm: the povm for which the distribution is desired

jVMC.operator.get_paulis(): Returns the Pauli matrices.

jVMC.operator.matrix_to_povm(A, M, T_inv, mode='unitary', dtype=<class 'numpy.float64'>)

Get operator from a matrix representation in POVM-formalism for the Lindblad equation or an observable.

In unitary mode this function implements \(\Omega^{ab} = -i T^{-1 bc} \mathrm{Tr}(A [M^c, M^a])\), in dissipative mode \(\Omega^{ab} = T^{-1 bc} \mathrm{Tr}(A M^c A^\dagger M^a- 1/2 A^\dagger A \{M^c, M^a\})\), in observable mode \(O^a = T^{-1 ab} \mathrm{Tr}(M^b A)\) and in imaginary time mode \(\Omega^{ab} = - T^{-1 bc} \mathrm{Tr}(A \{M^c, M^a\})\) where the Einstein sum convention is assumed.

Args:

A: Matrix representation of desired operator
M: POVM-Measurement operators
T_inv: Inverse POVM-Overlap matrix
mode: String specifying the conversion mode. Possible values are ‘unitary’ (‘uni’), ‘dissipative’ (‘dis’), ‘observable’ (‘obs’) and ‘imaginary’ (‘imag’)

Returns:

jax.numpy.ndarray