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 blockdiagonal structure of certain matrices, which in turn speeds up SVD or diagonalization a lot. Even for more general (nonsquarematrix) tensors, charge conservation imposes restrictions which blocks of a tensor can be nonzero. 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 predefined 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 innp_conserved
(the name standing for “numpy with charge conservation”). Internally, it stores only nonzero blocks of the tensor, which are “compatible” with the charges of the indices. It has to have a well defined overall chargeqtotal
. 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 aLegCharge
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()
oreigh()
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
) withcombine_legs()
, orto 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 dicta.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 (callingget_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: leta.labels = {'a': 1, 'b': 2, 'c': 3}
. Then ifb = a.combine_legs([[0, 1], [2]])
, it will haveb.labels = {'(a.b)': 0, '(c)': 1}
. If some subleg 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 rank3 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 axisan
Ellipsis
or...
: shorthand forslice(None)
for missing axes to fix the lenan 1D bool ndarray
mask
: apply a mask to that axis, seeiproject()
.a
slice(start, stop, step)
orstart: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 docstrings 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 multidimensional array representing a tensor with the entries:
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_i1\rbrace\).
The rank is the number of legs, the shape is \((n_0, ..., n_{rank1})\).
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 onetoone 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 legchargesorted, 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 nonzero?.
For completeness, let us also summarize also the internal structure of an Array
here:
The array saves only nonzero 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 qdatasorted 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 nonzero  details in Which entries of the npc Array can be nonzero?.
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 \lbrace1, +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 indices0, 1, 2, 3,..., 8
. You can identify four charge blocksslice(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 indices1, 2
(which are inslice(1, 3)
) have the same charge value[1]
. A qindex would just enumerate these blocks as0, 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 (here9
) as last entries. Thus it has lenblock_number + 1
, whereblock_number
(given byblock_number
) is the number of charge blocks in the leg, i.e. a qindex runs from 0 toblock_number1
. On the other hand, the 2D array charges has shape(block_number, qnumber)
, whereqnumber
is the number of charges (given byqnumber
).In that way, the qind form maps an qindex, say
qi
, to the indicesslice(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 

















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 nonzero?
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 nonzero
(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 nonzero.
All other indices are incompatible with the charges and must be zero.
In case of multiple charges, qnumber > 1, is a straightforward generalization: an entry can only be nonzero 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 nonzero entries?
Or rather, how can we ensure that C
complies with the above rule?
An entry C[ia,ic]
will only be nonzero,
if there is an ib
such that both A[ia,ib]
and B[ib,ic]
are nonzero, 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 straightforward 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 nonphysical 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 nonzero 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 nonzero 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 nonzero 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()
).
(Nonzero) 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 docstrings.
(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 nonzero 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 (nonzero) 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 nonzero 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 useget_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 returnsslice(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 returnsleg.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 rightmost 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.
See also
The module
tenpy.linalg.np_conserved
should contain all the API needed from the point of view of the algorithms. It contains the fundamentalArray
class and functions for working with them (creating and manipulating).The module
tenpy.linalg.charges
contains implementations for the charge structure, for example the classesChargeInfo
,LegCharge
, andLegPipe
. As noted above, the ‘public’ API is imported to (and accessible from)np_conserved
.
A full example code for spin1/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 nearestneighbour AFM Heisenberg Hamiltonian
3) define 'environments' left and right
4) contract MPS and MPO to calculate the energy
5) extract twosite 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, GNU GPLv3
import tenpy.linalg.np_conserved as npc
import numpy as np
# model parameters
Jxx, Jz = 1., 1.
L = 20
dt = 0.1
cutoff = 1.e10
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) # (rightcanonical)
Ss = [np.ones(1)] * L # Ss[i] are singular values between Bs[i1] and Bs[i]
# Side remark:
# An MPS is expected to have nonzero 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 spin1/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 <psiHpsi>")
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 twosite 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 singlesite 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")