Charge conservation with np_conserved

The basic idea is quickly summarized: By inspecting the Hamiltonian, you can identify symmetries, which correspond to conserved quantities, called charges. These charges divide the tensors into different sectors. This can be used to infer for example a block-diagonal structure of certain matrices, which in turn speeds up SVD or diagonalization a lot. Even for more general (non-square-matrix) tensors, charge conservation imposes restrictions which blocks of a tensor can be non-zero. Only those blocks need to be saved, which ultimately (= for large enough arrays) leads to a speedup of many routines, e.g., tensordot.

This introduction covers our implementation of charges; explaining mathematical details of the underlying symmetry is beyond its scope. We refer you to the corresponding chapter in our [TeNPyNotes] for a more general introduction of the idea (also stating the “charge rule” introduced below). [singh2010] explains why it works form a mathematical point of view, [singh2011] has the focus on a \(U(1)\) symmetry and might be easier to read.

What you really need to know about np_conserved

The good news is: It is not necessary to understand all the details explained in the following sections if you just want to use TeNPy for “standard” simulations like TEBD and DMRG. In praxis, you will likely not have to define the charges by yourself. For most simulations using TeNPy, the charges are initially defined in the Site; and there are many pre-defined sites like the :class:SpinHalfSite, which you can just use. The sites in turn are initialized by the Model class you are using (see also Models). From there, all the necessary charge information is automatically propagated along with the tensors.

However, you should definitely know a few basic facts about the usage of charge conservation in TeNPy:

Instead of using numpy arrays, tensors are represented by the Array class. This class is defined in np_conserved (the name standing for “numpy with charge conservation”). Internally, it stores only non-zero blocks of the tensor, which are “compatible” with the charges of the indices. It has to have a well defined overall charge qtotal. This excludes certain operators (like \(S^x\) for Sz conservation) and MPS which are a superpositions of states in different charge sectors.
There is a class ChargeInfo holding the general information what kind of charges we have, and a LegCharge for the charge data on a given leg. The leg holds a flag qconj which is +1 or -1, depending on whether the leg goes into the tensor (representing a vector space) or out of the tensor (representing the corresponding dual vector space).
Besides the array class methods, there are a bunch of functions like tensordot(), svd() or eigh() to manipulate tensors. These function have a very similar call structure as the corresponding numpy functions, but they act on our tensor Array class, and preserve the block structure (and exploit it for speed, wherever possible).
The only allowed “reshaping” operations for those tensors are to combine legs and to split previously combined legs. See the corresponding section below.
It is convenient to use string labels instead of numbers to refer to the various legs of a tensor. The rules how these labels change during the various operations are also described a section below.

Introduction to combine_legs, split_legs and LegPipes

Often, it is necessary to “combine” multiple legs into one: for example to perform a SVD, a tensor needs to be viewed as a matrix. For a flat array, this can be done with np.reshape, e.g., if A has shape (10, 3, 7) then B = np.reshape(A, (30, 7)) will result in a (view of the) array with one less dimension, but a “larger” first leg. By default (order='C'), this results in

B[i*3 + j , k] == A[i, j, k] for i in range(10) for j in range(3) for k in range(7)

While for a np.array, also a reshaping (10, 3, 7) -> (2, 21, 5) would be allowed, it does not make sense physically. The only sensible “reshape” operation on an Array are

to combine multiple legs into one leg pipe (LegPipe) with combine_legs(), or
to split a pipe of previously combined legs with split_legs().

Each leg has a Hilbert space, and a representation of the symmetry on that Hilbert space. Combining legs corresponds to the tensor product operation, and for abelian groups, the corresponding “fusion” of the representation is the simple addition of charge.

Fusion is not a lossless process, so if we ever want to split the combined leg, we need some additional data to tell us how to reverse the tensor product. This data is saved in the class LegPipe, derived from the LegCharge and used as new leg. Details of the information contained in a LegPipe are given in the class doc string.

The rough usage idea is as follows:

You can call combine_legs() without supplying any LegPipes, combine_legs will then make them for you.

Nevertheless, if you plan to perform the combination over and over again on sets of legs you know to be identical [with same charges etc, up to an overall -1 in qconj on all incoming and outgoing Legs] you might make a LegPipe anyway to save on the overhead of computing it each time.
In any way, the resulting Array will have a LegPipe as a LegCharge on the combined leg. Thus, it – and all tensors inheriting the leg (e.g. the results of svd, tensordot etc.) – will have the information how to split the LegPipe back to the original legs.
Once you performed the necessary operations, you can call split_legs(). This uses the information saved in the LegPipe to split the legs, recovering the original legs.

For a LegPipe, conj() changes qconj for the outgoing pipe and the incoming legs. If you need a LegPipe with the same incoming qconj, use outer_conj().

Leg labeling

It’s convenient to name the legs of a tensor: for instance, we can name legs 0, 1, 2 to be 'a', 'b', 'c': \(T_{i_a,i_b,i_c}\). That way we don’t have to remember the ordering! Under tensordot, we can then call

U = npc.tensordot(S, T, axes = [ [...],  ['b'] ] )

without having to remember where exactly 'b' is. Obviously U should then inherit the name of its legs from the uncontracted legs of S and T. So here is how it works:

Labels can only be strings. The labels should not include the characters . or ?. Internally, the labels are stored as dict a.labels = {label: leg_position, ...}. Not all legs need a label.
To set the labels, call
```
A.set_labels(['a', 'b', None, 'c', ... ])
```
which will set up the labeling {'a': 0, 'b': 1, 'c': 3 ...}.
(Where implemented) the specification of axes can use either the labels or the index positions. For instance, the call tensordot(A, B, [['a', 2, 'c'], [...]]) will interpret 'a' and 'c' as labels (calling get_leg_indices() to find their positions using the dict) and 2 as ‘the 2nd leg’. That’s why we require labels to be strings!
Labels will be intelligently inherited through the various operations of np_conserved.
- Under transpose, labels are permuted.
- Under tensordot, labels are inherited from uncontracted legs. If there is a collision, both labels are dropped.
- Under combine_legs, labels get concatenated with a . delimiter and surrounded by brackets. Example: let a.labels = {'a': 1, 'b': 2, 'c': 3}. Then if b = a.combine_legs([[0, 1], [2]]), it will have b.labels = {'(a.b)': 0, '(c)': 1}. If some sub-leg of a combined leg isn’t named, then a '?#' label is inserted (with # the leg index), e.g., 'a.?0.c'.
- Under split_legs, the labels are split using the delimiters (and the '?#' are dropped).
- Under conj, iconj: take 'a' -> 'a*', 'a*' -> 'a', and '(a.(b*.c))' -> '(a*.(b.c*))'
- Under svd, the outer labels are inherited, and inner labels can be optionally passed.
- Under pinv, the labels are transposed.

Indexing of an Array

Although it is usually not necessary to access single entries of an Array, you can of course do that. In the simplest case, this is something like A[0, 2, 1] for a rank-3 Array A. However, accessing single entries is quite slow and usually not recommended. For small Arrays, it may be convenient to convert them back to flat numpy arrays with to_ndarray().

On top of that very basic indexing, Array supports slicing and some kind of advanced indexing, which is however different from the one of numpy arrays (described here). Unlike numpy arrays, our Array class does not broadcast existing index arrays – this would be terribly slow. Also, np.newaxis is not supported, since inserting new axes requires additional information for the charges.

Instead, we allow just indexing of the legs independent of each other, of the form A[i0, i1, ...]. If all indices i0, i1, ... are integers, the single corresponding entry (of type dtype) is returned.

However, the individual ‘indices’ i0 for the individual legs can also be one of what is described in the following list. In that case, a new Array with less data (specified by the indices) is returned.

The ‘indices’ can be:

an int: fix the index of that axis, return array with one less dimension. See also take_slice().
a slice(None) or :: keep the complete axis
an Ellipsis or ...: shorthand for slice(None) for missing axes to fix the len
an 1D bool ndarray mask: apply a mask to that axis, see iproject().
a slice(start, stop, step) or start:stop:step: keep only the indices specified by the slice. This is also implemented with iproject.
an 1D int ndarray mask: keep only the indices specified by the array. This is also implemented with iproject.

For slices and 1D arrays, additional permutations may be performed with the help of permute().

If the number of indices is less than rank, the remaining axes remain free, so for a rank 4 Array A, A[i0, i1] == A[i0, i1, ...] == A[i0, i1, :, :].

Note that indexing always copies the data – even if int contains just slices, in which case numpy would return a view. However, assigning with A[:, [3, 5], 3] = B should work as you would expect.

Warning

Due to numpy’s advanced indexing, for 1D integer arrays a0 and a1 the following holds

A[a0, a1].to_ndarray() == A.to_ndarray()[np.ix_(a0, a1)] != A.to_ndarray()[a0, a1]

For a combination of slices and arrays, things get more complicated with numpys advanced indexing. In that case, a simple np.ix_(...) doesn’t help any more to emulate our version of indexing.

Details of the np_conserved implementation

Notations

Lets fix the notation of certain terms for this introduction and the doc-strings in np_conserved. This might be helpful if you know the basics from a different context. If you’re new to the subject, keep reading even if you don’t understand each detail, and come back to this section when you encounter the corresponding terms again.

A Array is a multi-dimensional array representing a tensor with the entries:

\[T_{a_0, a_1, ... a_{rank-1}} \quad \text{ with } \quad a_i \in \lbrace 0, ..., n_i-1 \rbrace\]

Each leg \(a_i\) corresponds the a vector space of dimension n_i.

An index of a leg is a particular value \(a_i \in \lbrace 0, ... ,n_i-1\rbrace\).

The rank is the number of legs, the shape is \((n_0, ..., n_{rank-1})\).

We restrict ourselves to abelian charges with entries in \(\mathbb{Z}\) or in \(\mathbb{Z}_m\). The nature of a charge is specified by \(m\); we set \(m=1\) for charges corresponding to \(\mathbb{Z}\). The number of charges is referred to as qnumber as a short hand, and the collection of \(m\) for each charge is called qmod. The qnumber, qmod and possibly descriptive names of the charges are saved in an instance of ChargeInfo.

To each index of each leg, a value of the charge(s) is associated. A charge block is a contiguous slice corresponding to the same charge(s) of the leg. A qindex is an index in the list of charge blocks for a certain leg. A charge sector is for given charge(s) is the set of all qindices of that charge(s). A leg is blocked if all charge sectors map one-to-one to qindices. Finally, a leg is sorted, if the charges are sorted lexicographically. Note that a sorted leg is always blocked. We can also speak of the complete array to be blocked by charges or legcharge-sorted, which means that all of its legs are blocked or sorted, respectively. The charge data for a single leg is collected in the class LegCharge. A LegCharge has also a flag qconj, which tells whether the charges point inward (+1) or outward (-1). What that means, is explained later in Which entries of the npc Array can be non-zero?.

For completeness, let us also summarize also the internal structure of an Array here: The array saves only non-zero blocks, collected as a list of np.array in self._data. The qindices necessary to map these blocks to the original leg indices are collected in self._qdata An array is said to be qdata-sorted if its self._qdata is lexicographically sorted. More details on this follow later. However, note that you usually shouldn’t access _qdata and _data directly - this is only necessary from within tensordot, svd, etc. Also, an array has a total charge, defining which entries can be non-zero - details in Which entries of the npc Array can be non-zero?.

Finally, a leg pipe (implemented in LegPipe) is used to formally combine multiple legs into one leg. Again, more details follow later.

Physical Example

For concreteness, you can think of the Hamiltonian \(H = -t \sum_{<i,j>} (c^\dagger_i c_j + H.c.) + U n_i n_j\) with \(n_i = c^\dagger_i c_i\). This Hamiltonian has the global \(U(1)\) gauge symmetry \(c_i \rightarrow c_i e^{i\phi}\). The corresponding charge is the total number of particles \(N = \sum_i n_i\). You would then introduce one charge with \(m=1\).

Note that the total charge is a sum of local terms, living on single sites. Thus, you can infer the charge of a single physical site: it’s just the value \(q_i = n_i \in \mathbb{N}\) for each of the states.

Note that you can only assign integer charges. Consider for example the spin 1/2 Heisenberg chain. Here, you can naturally identify the magnetization \(S^z = \sum_i S^z_i\) as the conserved quantity, with values \(S^z_i = \pm \frac{1}{2}\). Obviously, if \(S^z\) is conserved, then so is \(2 S^z\), so you can use the charges \(q_i = 2 S^z_i \in \lbrace-1, +1 \rbrace\) for the down and up states, respectively. Alternatively, you can also use a shift and define \(q_i = S^z_i + \frac{1}{2} \in \lbrace 0, 1 \rbrace\).

As another example, consider BCS like terms \(\sum_k (c^\dagger_k c^\dagger_{-k} + H.c.)\). These terms break the total particle conservation, but they preserve the total parity, i.e., \(N \mod 2\) is conserved. Thus, you would introduce a charge with \(m = 2\) in this case.

In the above examples, we had only a single charge conserved at a time, but you might be lucky and have multiple conserved quantities, e.g. if you have two chains coupled only by interactions. TeNPy is designed to handle the general case of multiple charges. When giving examples, we will restrict to one charge, but everything generalizes to multiple charges.

The different formats for LegCharge

As mentioned above, we assign charges to each index of each leg of a tensor. This can be done in three formats: qflat, as qind and as qdict. Let me explain them with examples, for simplicity considering only a single charge (the most inner array has one entry for each charge).

qflat form: simply a list of charges for each index.

An example:

qflat = [[-2], [-1], [-1], [0], [0], [0], [0], [3], [3]]

This tells you that the leg has size 9, the charges for are [-2], [-1], [-1], ..., [3] for the indices 0, 1, 2, 3,..., 8. You can identify four charge blocks slice(0, 1), slice(1, 3), slice(3, 7), slice(7, 9) in this example, which have charges [-2], [-1], [0], [3]. In other words, the indices 1, 2 (which are in slice(1, 3)) have the same charge value [-1]. A qindex would just enumerate these blocks as 0, 1, 2, 3.

qind form: a 1D array slices and a 2D array charges.

This is a more compact version than the qflat form: the slices give a partition of the indices and the charges give the charge values. The same example as above would simply be:

slices = [0, 1, 3, 7, 9]
charges = [[-2], [-1], [0], [3]]

Note that slices includes 0 as first entry and the number of indices (here 9) as last entries. Thus it has len block_number + 1, where block_number (given by block_number) is the number of charge blocks in the leg, i.e. a qindex runs from 0 to block_number-1. On the other hand, the 2D array charges has shape (block_number, qnumber), where qnumber is the number of charges (given by qnumber).

In that way, the qind form maps an qindex, say qi, to the indices slice(slices[qi], slices[qi+1]) and the charge(s) charges[qi].

qdict form: a dictionary in the other direction than qind, taking charge tuples to slices.

Again for the same example:

{(-2,): slice(0, 1),
 (-1,): slice(1, 3),
 (0,) : slice(3, 7),
 (3,) : slice(7, 9)}

Since the keys of a dictionary are unique, this form is only possible if the leg is completely blocked.

The LegCharge saves the charge data of a leg internally in qind form, directly in the attribute slices and charges. However, it also provides convenient functions for conversion between from and to the qflat and qdict form.

The above example was nice since all charges were sorted and the charge blocks were ‘as large as possible’. This is however not required.

The following example is also a valid qind form:

slices = [0, 1, 3, 5, 7, 9]
charges = [[-2], [-1], [0], [0], [3]]

This leads to the same qflat form as the above examples, thus representing the same charges on the leg indices. However, regarding our Arrays, this is quite different, since it divides the leg into 5 (instead of previously 4) charge blocks. We say the latter example is not bunched, while the former one is bunched.

To make the different notions of sorted and bunched clearer, consider the following (valid) examples:

charges	bunched	sorted	blocked
`[[-2], [-1], [0], [1], [3]]`	`True`	`True`	`True`
`[[-2], [-1], [0], [0], [3]]`	`False`	`True`	`False`
`[[-2], [0], [-1], [1], [3]]`	`True`	`False`	`True`
`[[-2], [0], [-1], [0], [3]]`	`True`	`False`	`False`

If a leg is bunched and sorted, it is automatically blocked (but not vice versa). See also below for further comments on that.

Which entries of the npc Array can be non-zero?

The reason for the speedup with np_conserved lies in the fact that it saves only the blocks ‘compatible’ with the charges. But how is this ‘compatible’ defined?

Assume you have a tensor, call it \(T\), and the LegCharge for all of its legs, say \(a, b, c, ...\).

Remember that the LegCharge associates to each index of the leg a charge value (for each of the charges, if qnumber > 1). Let a.to_qflat()[ia] denote the charge(s) of index ia for leg a, and similar for other legs.

In addition, the LegCharge has a flag qconj. This flag qconj is only a sign, saved as +1 or -1, specifying whether the charges point ‘inward’ (+1, default) or ‘outward’ (-1) of the tensor.

Then, the total charge of an entry T[ia, ib, ic, ...] of the tensor is defined as:

qtotal[ia, ib, ic, ...] = a.to_qflat()[ia] * a.qconj + b.to_qflat()[ib] * b.qconj + c.to_qflat()[ic] * c.qconj + ...  modulo qmod

The rule which entries of the a Array can be non-zero (i.e., are ‘compatible’ with the charges), is then very simple:

In other words, there is a single total charge .qtotal attribute of a Array. All indices ia, ib, ic, ... for which the above defined qtotal[ia, ib, ic, ...] matches this total charge, are said to be compatible with the charges and can be non-zero. All other indices are incompatible with the charges and must be zero.

In case of multiple charges, qnumber > 1, is a straight-forward generalization: an entry can only be non-zero if it is compatible with each of the defined charges.

The pesky qconj - contraction as an example

Why did we introduce the qconj flag? Remember it’s just a sign telling whether the charge points inward or outward. So whats the reasoning?

The short answer is, that LegCharges actually live on bonds (i.e., legs which are to be contracted) rather than individual tensors. Thus, it is convenient to share the LegCharges between different legs and even tensors, and just adjust the sign of the charges with qconj.

As an example, consider the contraction of two tensors, \(C_{ia,ic} = \sum_{ib} A_{ia,ib} B_{ib,ic}\). For simplicity, say that the total charge of all three tensors is zero. What are the implications of the above rule for non-zero entries? Or rather, how can we ensure that C complies with the above rule? An entry C[ia,ic] will only be non-zero, if there is an ib such that both A[ia,ib] and B[ib,ic] are non-zero, i.e., both of the following equations are fulfilled:

A.qtotal == A.legs[0].to_qflat()[ia] * A.legs[0].qconj + A.legs[1].to_qflat()[ib] * A.legs[1].qconj  modulo qmod
B.qtotal == B.legs[0].to_qflat()[ib] * B.legs[0].qconj + B.legs[1].to_qflat()[ic] * B.legs[1].qconj  modulo qmod

(A.legs[0] is the LegCharge saving the charges of the first leg (with index ia) of A.)

For the uncontracted legs, we just keep the charges as they are:

C.legs = [A.legs[0], B.legs[1]]

It is then straight-forward to check, that the rule is fulfilled for \(C\), if the following condition is met:

A.qtotal + B.qtotal - C.qtotal == A.legs[1].to_qflat()[ib] A.b.qconj + B.legs[0].to_qflat()[ib] B.b.qconj  modulo qmod

The easiest way to meet this condition is (1) to require that A.b and B.b share the same charges b.to_qflat(), but have opposite qconj, and (2) to define C.qtotal = A.qtotal + B.qtotal. This justifies the introduction of qconj: when you define the tensors, you have to define the LegCharge for the b only once, say for A.legs[1]. For B.legs[0] you simply use A.legs[1].conj() which creates a copy of the LegCharge with shared slices and charges, but opposite qconj. As a more impressive example, all ‘physical’ legs of an MPS can usually share the same LegCharge (up to different qconj if the local Hilbert space is the same). This leads to the following convention:

Assigning charges to non-physical legs

From the above physical examples, it should be clear how you assign charges to physical legs. But what about other legs, e.g, the virtual bond of an MPS (or an MPO)?

The charge of these bonds must be derived by using the ‘rule for non-zero entries’, as far as they are not arbitrary. As a concrete example, consider an MPS on just two spin 1/2 sites:

|        _____         _____
|   x->- | A | ->-y->- | B | ->-z
|        -----         -----
|          ^             ^
|          |p            |p

The two legs p are the physical legs and share the same charge, as they both describe the same local Hilbert space. For better distinction, let me label the indices of them by \(\uparrow=0\) and \(\downarrow=1\). As noted above, we can associate the charges 1 (\(p=\uparrow\)) and -1 (\(p=\downarrow\)), respectively, so we define:

chinfo = npc.ChargeInfo([1], ['2*Sz'])
p  = npc.LegCharge.from_qflat(chinfo, [1, -1], qconj=+1)

For the qconj signs, we stick to the convention used in our MPS code and indicated by the arrows in above ‘picture’: physical legs are incoming (qconj=+1), and from left to right on the virtual bonds. This is achieved by using [p, x, y.conj()] as legs for A, and [p, y, z.conj()] for B, with the default qconj=+1 for all p, x, y, z: y.conj() has the same charges as y, but opposite qconj=-1.

The legs x and z of an L=2 MPS, are ‘dummy’ legs with just one index 0. The charge on one of them, as well as the total charge of both A and B is arbitrary (i.e., a gauge freedom), so we make a simple choice: total charge 0 on both arrays, as well as for \(x=0\), x = npc.LegCharge.from_qflat(chinfo, [0], qconj=+1).

The charges on the bonds y and z then depend on the state the MPS represents. Here, we consider a singlet \(\psi = (|\uparrow \downarrow\rangle - |\downarrow \uparrow\rangle)/\sqrt{2}\) as a simple example. A possible MPS representation is given by:

A[up, :, :]   = [[1/2.**0.5, 0]]     B[up, :, :]   = [[0], [-1]]
A[down, :, :] = [[0, 1/2.**0.5]]     B[down, :, :] = [[1], [0]]

There are two non-zero entries in A, for the indices \((a, x, y) = (\uparrow, 0, 0)\) and \((\downarrow, 0, 1)\). For \((a, x, y) = (\uparrow, 0, 0)\), we want:

A.qtotal = 0 = p.to_qflat()[up] * p.qconj + x.to_qflat()[0] * x.qconj + y.conj().to_qflat()[0] * y.conj().qconj
             = 1                * (+1)    + 0               * (+1)    + y.conj().to_qflat()[0] * (-1)

This fixes the charge of y=0 to 1. A similar calculation for \((a, x, y) = (\downarrow, 0, 1)\) yields the charge -1 for y=1. We have thus all the charges of the leg y and can define y = npc.LegCharge.from_qflat(chinfo, [1, -1], qconj=+1).

Now take a look at the entries of B. For the non-zero entry \((b, y, z) = (\uparrow, 1, 0)\), we want:

B.qtotal = 0 = p.to_qflat()[up] * p.qconj + y.to_qflat()[1] * y.qconj + z.conj().to_qflat()[0] * z.conj().qconj
             = 1                * (+1)    + (-1)            * (+1)    + z.conj().to_qflat()[0] * (-1)

This implies the charge 0 for z = 0, thus z = npc.LegCharge.form_qflat(chinfo, [0], qconj=+1). Finally, note that the rule for \((b, y, z) = (\downarrow, 0, 0)\) is automatically fulfilled! This is an implication of the fact that the singlet has a well defined value for \(S^z_a + S^z_b\). For other states without fixed magnetization (e.g., \(|\uparrow \uparrow\rangle + |\downarrow \downarrow\rangle\)) this would not be the case, and we could not use charge conservation.

As an exercise, you can calculate the charge of z in the case that A.qtotal=5, B.qtotal = -1 and charge 2 for x=0. The result is -2.

Note

This section is meant be an pedagogical introduction. In you program, you can use the functions detect_legcharge() (which does exactly what’s described above) or detect_qtotal() (if you know all LegCharges, but not qtotal).

Array creation

Direct creation

Making an new Array requires both the tensor entries (data) and charge data.

The default initialization a = Array(...) creates an empty Array, where all entries are zero (equivalent to zeros()). (Non-zero) data can be provided either as a dense np.array to from_ndarray(), or by providing a numpy function such as np.random, np.ones etc. to from_func().

In both cases, the charge data is provided by one ChargeInfo, and a LegCharge instance for each of the legs.

Note

The charge data instances are not copied, in order to allow it to be shared between different Arrays. Consequently, you must make copies of the charge data, if you manipulate it directly. (However, methods like sort() do that for you.)

Indirect creation by manipulating existing arrays

Of course, a new Array can also created using the charge data from existing Arrays, for example with zeros_like() or creating a (deep or shallow) copy(). Further, there are many higher level functions like tensordot() or svd(), which also return new Arrays.

Complete blocking of Charges

While the code was designed in such a way that each charge sector has a different charge, the code should still run correctly if multiple charge sectors (for different qindex) correspond to the same charge. In this sense Array can act like a sparse array class to selectively store subblocks. Algorithms which need a full blocking should state that explicitly in their doc-strings. (Some functions (like svd and eigh) require complete blocking internally, but if necessary they just work on a temporary copy returned by as_completely_blocked()).

If you expect the tensor to be dense subject to charge constraints (as for MPS), it will be most efficient to fully block by charge, so that work is done on large chunks.

However, if you expect the tensor to be sparser than required by charge (as for an MPO), it may be convenient not to completely block, which forces smaller matrices to be stored, and hence many zeroes to be dropped. Nevertheless, the algorithms were not designed with this in mind, so it is not recommended in general. (If you want to use it, run a benchmark to check whether it is really faster!)

If you haven’t created the array yet, you can call sort() (with bunch=True) on each LegCharge which you want to block. This sorts by charges and thus induces a permutation of the indices, which is also returned as an 1D array perm. For consistency, you have to apply this permutation to your flat data as well.

Alternatively, you can simply call sort_legcharge() on an existing Array. It calls sort() internally on the specified legs and performs the necessary permutations directly to (a copy of) self. Yet, you should keep in mind, that the axes are permuted afterwards.

Internal Storage schema of npc Arrays

The actual data of the tensor is stored in _data. Rather than keeping a single np.array (which would have many zeros in it), we store only the non-zero sub blocks. So _data is a python list of np.array’s. The order in which they are stored in the list is not physically meaningful, and so not guaranteed (more on this later). So to figure out where the sub block sits in the tensor, we need the _qdata structure (on top of the LegCharges in legs).

Consider a rank 3 tensor T, with the first leg like:

legs[0].slices = np.array([0, 1, 4, ...])
legs[0].charges = np.array([[-2], [1], ...])

Each row of charges gives the charges for a charge block of the leg, with the actual indices of the total tensor determined by the slices. The qindex simply enumerates the charge blocks of a lex. Picking a qindex (and thus a charge block) from each leg, we have a subblock of the tensor.

For each (non-zero) subblock of the tensor, we put a (numpy) ndarray entry in the _data list. Since each subblock of the tensor is specified by rank qindices, we put a corresponding entry in _qdata, which is a 2D array of shape (#stored_blocks, rank). Each row corresponds to a non-zero subblock, and there are rank columns giving the corresponding qindex for each leg.

Example: for a rank 3 tensor we might have:

T._data = [t1, t2, t3, t4, ...]
T._qdata = np.array([[3, 2, 1],
                     [1, 1, 1],
                     [4, 2, 2],
                     [2, 1, 2],
                     ...       ])

The third subblock has an ndarray t3, and qindices [4 2 2] for the three legs.

To find the position of t3 in the actual tensor you can use get_slice():
```
T.legs[0].get_slice(4), T.legs[1].get_slice(2), T.legs[2].get_slice(2)
```
The function leg.get_charges(qi) simply returns slice(leg.slices[qi], leg.slices[qi+1])
To find the charges of t3, we an use get_charge():
```
T.legs[0].get_charge(2), T.legs[1].get_charge(2), T.legs[2].get_charge(2)
```
The function leg.get_charge(qi) simply returns leg.charges[qi]*leg.qconj.

Note

Outside of np_conserved, you should use the API to access the entries. If you really need to iterate over all blocks of an Array T, try for (block, blockslices, charges, qindices) in T: do_something().

The order in which the blocks stored in _data/_qdata is arbitrary (although of course _data and _qdata must be in correspondence). However, for many purposes it is useful to sort them according to some convention. So we include a flag ._qdata_sorted to the array. So, if sorted (with isort_qdata(), the _qdata example above goes to

_qdata = np.array([[1, 1, 1],
                   [3, 2, 1],
                   [2, 1, 2],
                   [4, 2, 2],
                   ...       ])

Note that np.lexsort chooses the right-most column to be the dominant key, a convention we follow throughout.

If _qdata_sorted == True, _qdata and _data are guaranteed to be lexsorted. If _qdata_sorted == False, there is no guarantee. If an algorithm modifies _qdata, it must set _qdata_sorted = False (unless it guarantees it is still sorted). The routine sort_qdata() brings the data to sorted form.

A full example code for spin-1/2

Below follows a full example demonstrating the creation and contraction of Arrays. (It’s the file a_np_conserved.py in the examples folder of the tenpy source.)

"""An example code to demonstrate the usage of :class:`~tenpy.linalg.np_conserved.Array`.

This example includes the following steps:
1) create Arrays for an Neel MPS
2) create an MPO representing the nearest-neighbour AFM Heisenberg Hamiltonian
3) define 'environments' left and right
4) contract MPS and MPO to calculate the energy
5) extract two-site hamiltonian ``H2`` from the MPO
6) calculate ``exp(-1.j*dt*H2)`` by diagonalization of H2
7) apply ``exp(H2)`` to two sites of the MPS and truncate with svd

Note that this example uses only np_conserved, but no other modules.
Compare it to the example `b_mps.py`,
which does the same steps using a few predefined classes like MPS and MPO.
"""
# Copyright (C) TeNPy Developers, Apache license

import tenpy.linalg.np_conserved as npc
import numpy as np

# model parameters
Jxx, Jz = 1., 1.
L = 20
dt = 0.1
cutoff = 1.e-10
print("Jxx={Jxx}, Jz={Jz}, L={L:d}".format(Jxx=Jxx, Jz=Jz, L=L))

print("1) create Arrays for an Neel MPS")

#  vL ->--B-->- vR
#         |
#         ^
#         |
#         p

# create a ChargeInfo to specify the nature of the charge
chinfo = npc.ChargeInfo([1], ['2*Sz'])  # the second argument is just a descriptive name

# create LegCharges on physical leg and even/odd bonds
p_leg = npc.LegCharge.from_qflat(chinfo, [[1], [-1]])  # charges for up, down
v_leg_even = npc.LegCharge.from_qflat(chinfo, [[0]])
v_leg_odd = npc.LegCharge.from_qflat(chinfo, [[1]])

B_even = npc.zeros([v_leg_even, v_leg_odd.conj(), p_leg],
                   labels=['vL', 'vR', 'p'])  # virtual left/right, physical
B_odd = npc.zeros([v_leg_odd, v_leg_even.conj(), p_leg], labels=['vL', 'vR', 'p'])
B_even[0, 0, 0] = 1.  # up
B_odd[0, 0, 1] = 1.  # down

Bs = [B_even, B_odd] * (L // 2) + [B_even] * (L % 2)  # (right-canonical)
Ss = [np.ones(1)] * L  # Ss[i] are singular values between Bs[i-1] and Bs[i]

# Side remark:
# An MPS is expected to have non-zero entries everywhere compatible with the charges.
# In general, we recommend to use `sort_legcharge` (or `as_completely_blocked`)
# to ensure complete blocking. (But the code will also work, if you don't do it.)
# The drawback is that this might introduce permutations in the indices of single legs,
# which you have to keep in mind when converting dense numpy arrays to and from npc.Arrays.

print("2) create an MPO representing the AFM Heisenberg Hamiltonian")

#         p*
#         |
#         ^
#         |
#  wL ->--W-->- wR
#         |
#         ^
#         |
#         p

# create physical spin-1/2 operators Sz, S+, S-
Sz = npc.Array.from_ndarray([[0.5, 0.], [0., -0.5]], [p_leg, p_leg.conj()], labels=['p', 'p*'])
Sp = npc.Array.from_ndarray([[0., 1.], [0., 0.]], [p_leg, p_leg.conj()], labels=['p', 'p*'])
Sm = npc.Array.from_ndarray([[0., 0.], [1., 0.]], [p_leg, p_leg.conj()], labels=['p', 'p*'])
Id = npc.eye_like(Sz, labels=Sz.get_leg_labels())  # identity

mpo_leg = npc.LegCharge.from_qflat(chinfo, [[0], [2], [-2], [0], [0]])

W_grid = [[Id,   Sp,   Sm,   Sz,   None          ],
          [None, None, None, None, 0.5 * Jxx * Sm],
          [None, None, None, None, 0.5 * Jxx * Sp],
          [None, None, None, None, Jz * Sz       ],
          [None, None, None, None, Id            ]]  # yapf:disable

W = npc.grid_outer(W_grid, [mpo_leg, mpo_leg.conj()], grid_labels=['wL', 'wR'])
# wL/wR = virtual left/right of the MPO
Ws = [W] * L

print("3) define 'environments' left and right")

#  .---->- vR     vL ->----.
#  |                       |
#  envL->- wR     wL ->-envR
#  |                       |
#  .---->- vR*    vL*->----.

envL = npc.zeros([W.get_leg('wL').conj(), Bs[0].get_leg('vL').conj(), Bs[0].get_leg('vL')],
                 labels=['wR', 'vR', 'vR*'])
envL[0, :, :] = npc.diag(1., envL.legs[1])
envR = npc.zeros([W.get_leg('wR').conj(), Bs[-1].get_leg('vR').conj(), Bs[-1].get_leg('vR')],
                 labels=['wL', 'vL', 'vL*'])
envR[-1, :, :] = npc.diag(1., envR.legs[1])

print("4) contract MPS and MPO to calculate the energy <psi|H|psi>")
contr = envL
for i in range(L):
    # contr labels: wR, vR, vR*
    contr = npc.tensordot(contr, Bs[i], axes=('vR', 'vL'))
    # wR, vR*, vR, p
    contr = npc.tensordot(contr, Ws[i], axes=(['p', 'wR'], ['p*', 'wL']))
    # vR*, vR, wR, p
    contr = npc.tensordot(contr, Bs[i].conj(), axes=(['p', 'vR*'], ['p*', 'vL*']))
    # vR, wR, vR*
    # note that the order of the legs changed, but that's no problem with labels:
    # the arrays are automatically transposed as necessary
E = npc.inner(contr, envR, axes=(['vR', 'wR', 'vR*'], ['vL', 'wL', 'vL*']))
print("E =", E)

print("5) calculate two-site hamiltonian ``H2`` from the MPO")
# label left, right physical legs with p, q
W0 = W.replace_labels(['p', 'p*'], ['p0', 'p0*'])
W1 = W.replace_labels(['p', 'p*'], ['p1', 'p1*'])
H2 = npc.tensordot(W0, W1, axes=('wR', 'wL')).itranspose(['wL', 'wR', 'p0', 'p1', 'p0*', 'p1*'])
H2 = H2[0, -1]  # (If H has single-site terms, it's not that simple anymore)
print("H2 labels:", H2.get_leg_labels())

print("6) calculate exp(H2) by diagonalization of H2")
# diagonalization requires to view H2 as a matrix
H2 = H2.combine_legs([('p0', 'p1'), ('p0*', 'p1*')], qconj=[+1, -1])
print("labels after combine_legs:", H2.get_leg_labels())
E2, U2 = npc.eigh(H2)
print("Eigenvalues of H2:", E2)
U_expE2 = U2.scale_axis(np.exp(-1.j * dt * E2), axis=1)  # scale_axis ~= apply an diagonal matrix
exp_H2 = npc.tensordot(U_expE2, U2.conj(), axes=(1, 1))
exp_H2.iset_leg_labels(H2.get_leg_labels())
exp_H2 = exp_H2.split_legs()  # by default split all legs which are `LegPipe`
# (this restores the originial labels ['p0', 'p1', 'p0*', 'p1*'] of `H2` in `exp_H2`)

print("7) apply exp(H2) to even/odd bonds of the MPS and truncate with svd")
# (this implements one time step of first order TEBD)
for even_odd in [0, 1]:
    for i in range(even_odd, L - 1, 2):
        B_L = Bs[i].scale_axis(Ss[i], 'vL').ireplace_label('p', 'p0')
        B_R = Bs[i + 1].replace_label('p', 'p1')
        theta = npc.tensordot(B_L, B_R, axes=('vR', 'vL'))
        theta = npc.tensordot(exp_H2, theta, axes=(['p0*', 'p1*'], ['p0', 'p1']))
        # view as matrix for SVD
        theta = theta.combine_legs([('vL', 'p0'), ('p1', 'vR')], new_axes=[0, 1], qconj=[+1, -1])
        # now theta has labels '(vL.p0)', '(p1.vR)'
        U, S, V = npc.svd(theta, inner_labels=['vR', 'vL'])
        # truncate
        keep = S > cutoff
        S = S[keep]
        invsq = np.linalg.norm(S)
        Ss[i + 1] = S / invsq
        U = U.iscale_axis(S / invsq, 'vR')
        Bs[i] = U.split_legs('(vL.p0)').iscale_axis(Ss[i]**(-1), 'vL').ireplace_label('p0', 'p')
        Bs[i + 1] = V.split_legs('(p1.vR)').ireplace_label('p1', 'p')
print("finished")

Charge conservation with np_conserved

What you really need to know about np_conserved

Introduction to combine_legs, split_legs and LegPipes

Leg labeling

Indexing of an Array

Details of the np_conserved implementation

Notations

Physical Example

The different formats for LegCharge

Which entries of the npc Array can be non-zero?

The pesky qconj - contraction as an example

Assigning charges to non-physical legs

Array creation

Direct creation

Indirect creation by manipulating existing arrays

Complete blocking of Charges

Internal Storage schema of npc Arrays

See also

A full example code for spin-1/2