Symmray

<img src="https://github.com/jcmgray/symmray/blob/HEAD/docs/_static/logo-title.png?raw=true" alt="symmray logo" width="640"/>

symmray is minimal library for block sparse, abelian symmetric and fermionic arrays, designed to look as much as possible like standard ndarrays, whose blocks can be backed by numpy, torch or any other autoray compatible library.

Installation

Installing the latest version directly from github:

If you want to checkout the latest version of features and fixes, you can install directly from the github repository:

pip install -U git+https://github.com/jcmgray/symmray.git

Installing a local, editable development version:

If you want to make changes to the source code and test them out, you can install a local editable version of the package:

git clone https://github.com/jcmgray/symmray.git
pip install --no-deps -U -e symmray/

Usage

symmray objects are designed so that, as much as possible, one can interact with them in the same way as standard arrays. You can use the functions from the symmray namespace directly:

import symmray as sr

# create some random arrays:
a, b, c, d, e, f = [sr.utils.rand_index("Z2", d) for d in [2, 3, 4, 5, 6, 7]]
x = sr.utils.get_rand("Z2", shape=[a, b, c.conj(), d], fermionic=True)
y = sr.utils.get_rand("Z2", shape=[e, a.conj(), f, c], fermionic=True)
x, y
# (Z2FermionicArray(shape~(2, 3, 4, 5):[++-+], charge=0, num_blocks=8),
#  Z2FermionicArray(shape~(6, 2, 7, 4):[+--+], charge=0, num_blocks=8))

# contract them:
z = sr.tensordot(x, y, axes=[(0, 2), (1, 3)])
z
# Z2FermionicArray(shape~(3, 5, 6, 7):[+++-], charge=0, num_blocks=8)

print(z)
# Z2FermionicArray(ndim=4, charge=0, indices=[
#     (3 = 2+1 : +[0,1])
#     (5 = 3+2 : +[0,1])
#     (6 = 2+4 : +[0,1])
#     (7 = 6+1 : -[0,1])
# ], num_blocks=8, backend=numpy, dtype=float64)

# fuse and decompose:
sr.linalg.svd(z.fuse((0, 3), (1, 2)))
# (Z2FermionicArray(shape~(21, 21):[++], charge=0, num_blocks=2),
#  BlockVector(total_size=21, num_blocks=2),
#  Z2FermionicArray(shape~(21, 30):[-+], charge=0, num_b locks=2))

or you can use the automatic dispatch library autoray to support multiple backends including symmray:

import autoray as ar

z = ar.do("tensordot", x, y, axes=[(0, 2), (1, 3)])

or you can use the python array api:

xp = x.__array_namespace__()
z = xp.tensordot(x, y, axes=[(0, 2), (1, 3)])

symmray also uses autoray internally to handle manipulating blocks within an array, meaning that these can be numpy, torch, jax or any other autoray compatible library.

Whilst block sparse arrays do not have such a well defined notion of shape as dense arrays, for ease and compatibility with other libraries, symmray arrays do have a .shape attribute which is the shape of the dense array that would be returned by calling to_dense on the array, and a similarly defined .size. Likewise, symmray supports fusing and unfusing of indices via reshape (as long as it is clear what is meant by the new shape).

Quick-start tensor networks and Hamiltonians

symmray provides constructors for various quimb.tensor.TensorNetwork networks:

• TN_abelian_from_edges_rand
• TN_fermionic_from_edges_rand
• PEPS_abelian_rand (2D specific)
• PEPS_fermionic_rand (2D specific)

Along with constructors for common hamiltonians:

• ham_tfim_from_edges
• ham_heisenberg_from_edges
• ham_fermi_hubbard_from_edges
• ham_fermi_hubbard_spinless_from_edges (AKA 't-V' model)

These constructors automatically chooose various defaults. See the examples folder for usage.

Block sparse abelian symmetric arrays

The core AbelianArray object consists of 4 main components:

1. .indices: a sequence of BlockIndex instances describing the charge distribution and 'dualness' of each dimension.
2. .charge: an overall charge for the array, which sets which combinations of index charges ('sectors') are allowed.
3. .blocks: a dict mapping each non-zero sector to a 'raw' array.
4. .symmetry: an object defining allowed charges and how they combine.

Specific subclasses of AbelianArray have a static .symmetry class attribute.

The BlockIndex object consists of 2 main components:

1. .chargemap: a dict mapping each charge to its size. The total size of the index is the sum of the sizes of all charges.

2. .dual: a boolean indicating whether the index is 'dual' or not. By convention:

   • dual=False means index flows 'outwards' / (+ve) / ket-like
   • dual=True means index flows 'inwards' / (-ve) / bra-like

One convenient way to create AbelianArray instances is via the from_fill_fn method, which takes a function with signature fn(shape) -> array_like and uses it to fill each valid sector of the array.

import symmray as sr
import numpy as np

indices = (
    sr.BlockIndex(chargemap={0: 3, 1: 4}, dual=False),
    sr.BlockIndex(chargemap={0: 5, 1: 6}, dual=True),
)

x = sr.Z2Array.from_fill_fn(
    fill_fn=np.ones,
    indices=indices,
    charge=1,
)
x
# Z2Array(shape~(7, 11):[+-], charge=1, num_blocks=2)

x.blocks
# {(0,
#   1): array([[1., 1., 1., 1., 1., 1.],
#         [1., 1., 1., 1., 1., 1.],
#         [1., 1., 1., 1., 1., 1.]]),
#  (1,
#   0): array([[1., 1., 1., 1., 1.],
#         [1., 1., 1., 1., 1.],
#         [1., 1., 1., 1., 1.],
#         [1., 1., 1., 1., 1.]])}

We can pictorially represent this like so:

You can also create AbelianArray instances using the methods:

• AbelianArray(indices, charge, blocks)
• AbelianArray.from_blocks(blocks, duals, charges) which calculates the index chargemaps from the sectors and blocks themselves.
• AbelianArray.random(indices, charge) which uses a random fill_fn
• AbelianArray.from_dense(array, index_maps, duals, charge) which converts a dense array to a block sparse array given a mapping for each axis, which specifies the charge of each linear index in the dense array.

Key functions which match the numpy versions are:

• conj
• reshape
• tensordot
• trace
• transpose

With additional symmray specific key functions:

• fuse
• multiply_diagonal

Fuse in particular is a crucial function for A. performing efficient contractions, B. performing linear algebra decompositions, as well as various other tensor network manipulations. You can either fuse and unfuse directly, or by using reshape. Note that if a symmray array is quite sparse (e.g. with U1 symmetry), then the resulting fused/reshaped shape will not necessarily match the dense specification.

The key function tensordot can use one of two methods.

• tensordot(x, y, axes, mode="fused"): fuse the two arrays into block diagonal matrices and then unfuse the result. This can be much faster, though possibly requires explicitly filling missing blocks with zeros.
• tensordot(x, y, axes, mode="blockwise"): compute the contraction by directly looping over the blocks of x and y and contracting them. This has quite high overhead for large numbers of blocks.

Fermionic arrays

The approach to fermionic arrays symmray takes is equivalent to the 'Grassmann' or graded algebra approach. This associates a fermionic parity to each charge, combined with the directionality specified by dual, allows all fermionic swaps and the relevant sector phase changes to be handled essentially locally.

The core FermionicArray is a subclass of AbelianArray and instantiated in the same way:

indices = (
    sr.BlockIndex(chargemap={-1: 2, 0: 2, 1: 3}, dual=False),
    sr.BlockIndex(chargemap={0: 2, 2: 3, 3: 4}, dual=True),
)

x = sr.U1FermionicArray.random(
    indices=indices,
    charge=-2,
)
print(x)
# U1FermionicArray(ndim=2, charge=-2, indices=[
#     (7 = 2+2+3 : +[-1,0,1])
#     (9 = 2+3+4 : -[0,2,3])
# ], num_blocks=2, backend=numpy, dtype=float64)

Phases are lazily tracked into the attribute .phases when:

• transposing
• fusing
• conjugating
• contracting via tensordot or __matmul__ / @
• tracing
• linear algebra decompositions

via the methods:

• FermionicArray.phase_flip: virtually insert 'parity' tensors on some axes
• FermionicArray.phase_transpose: compute the phase of a 'virtual' transpose

And inserted when needed using:

• FermionicArray.phase_sync: actually multiply the blocks by the phases.

Multiple odd-parity fermionic arrays - label and dummy_modes

If you want to work with networks involving multiple odd-parity tensors (and care about the global phase) then you must supply any sortable object as label to the FermionicArray constructor (typically lattice coordinate for example). This triggers the default creation of a single dummy mode, stored in the attribute dummy_modes : tuple[FermionicOperator, ...]. As fermionic arrays are contracted, these dummy modes are accumulated as if they were size 1 indices prepended to the array. The modes are also sorted according to their label and dualness, yielding a global phase with respect to the sorted ordering. Adjacent conjugate pairs can also be contracted away.

For example, if we create two fermionic arrays with different labels:

# in