diff --git a/docs/make.jl b/docs/make.jl index 34b025580..d72b851a6 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -9,7 +9,7 @@ end using Documenter using Random using TensorKit -using TensorKit: FusionTreePair, FusionTreeBlock, Index2Tuple +using TensorKit: FusionTreePair, FusionTreeBlock, Index2Tuple, IndexTuple using TensorKit.TensorKitSectors using TensorKit.MatrixAlgebraKit using DocumenterInterLinks @@ -26,8 +26,14 @@ pages = [ "man/intro.md", "man/tutorial.md", "man/spaces.md", "man/symmetries.md", "man/sectors.md", "man/gradedspaces.md", - "man/fusiontrees.md", "man/tensors.md", - "man/tensormanipulations.md", + "man/fusiontrees.md", + "Tensors" => [ + "man/tensors.md", + "man/linearalgebra.md", + "man/indexmanipulations.md", + "man/factorizations.md", + "man/contractions.md", + ], ], "Library" => [ "lib/sectors.md", "lib/fusiontrees.md", diff --git a/docs/src/appendix/categories.md b/docs/src/appendix/categories.md index 5a111b4fe..0bf48da85 100644 --- a/docs/src/appendix/categories.md +++ b/docs/src/appendix/categories.md @@ -18,7 +18,7 @@ To start, a **category** ``C`` consists of * composition of morphisms ``f:W→V`` and ``g:X→W`` into ``(f ∘ g):X→V`` that is associative, such that for ``h:Y→X`` we have ``f ∘ (g ∘ h) = (f ∘ g) ∘ h`` * for each object ``V``, an identity morphism ``\mathrm{id}_V:V→V`` such that ``f ∘ \mathrm{id}_W = f = \mathrm{id}_V ∘ f``. -The morphisms in ``\mathrm{Hom}_C(V,V)`` are known as endomorphism and this set is also denoted as ``End_C(V)``. +The morphisms in ``\mathrm{Hom}_C(V,V)`` are known as endomorphisms and this set is also denoted as ``End_C(V)``. When the category ``C`` is clear, we can drop the subscript in ``\mathrm{Hom}(W,V)``. A morphism ``f:W→V`` is an isomorphism if there exists a morphism ``f^{-1}:V→W`` called its inverse, such that ``f^{-1} ∘ f = \mathrm{id}_W`` and ``f ∘ f^{-1} = \mathrm{id}_V``. @@ -182,7 +182,7 @@ We come back to this point below. A left (or right) duality in a (monoidal) category is now defined as an association of a left (or right) dual with every object of the category, with corresponding exact pairings, and a category admitting such a duality is a left (or right) **rigid category** (or left or right autonomous category). Given that left (or right) morphism transposition satisfies ``{}^{∨}(f ∘ g)= {}^{∨}g ∘ {}^{∨}f= {}^{∨}f ∘^{\mathrm{op}} {}^{∨}g`` and recalling ``{}^{∨}(V ⊗ W) = {}^{∨}W ⊗ {}^{∨}V`` (and similar for right duality), we can define duality in a functorial way. A (left or right) rigid category ``\mathcal{C}`` is a category which admits a (left or right) duality functor, i.e. a functor from ``\mathcal{C}`` to ``\mathcal{C}^{\mathrm{rev}}`` that maps objects to its (left or right) dual, and morphisms to its (left or right) transpose. -In particular, the snake rules can now be read as the functioral requirement that ``{}^{∨}(\mathrm{id}_V) = \mathrm{id}_{{}^{∨}V}``. +In particular, the snake rules can now be read as the functorial requirement that ``{}^{∨}(\mathrm{id}_V) = \mathrm{id}_{{}^{∨}V}``. In all of this, left and right duality can be completely distinct. Equivalently, the left dual of the left dual of an object ``V``, i.e. ``{}^{∨∨}V`` is not necessarily ``V`` itself, nor do the exact pairings enable us to construct an isomorphism between ``{}^{∨∨}V`` and ``V``. @@ -257,7 +257,7 @@ In order for there to be a unique map from ``V_1 ⊗ V_2 ⁠⊗ … V_N`` to any The resulting category is then referred to as a **symmetric tensor category**. In a graphical representation, it means that there is no distinction between over- and under- crossings and, as such, lines can just cross, where the crossing represents the action of ``τ_{V,W} = τ_{W,V}^{-1}``. -In the case of the category ``\mathbf{Vect}`` a valid braiding consists of just flipping the the objects/morphisms involved, e.g. for a simple cartesian tensor, permuting the tensor indices is equivalent to applying Julia's function `permutedims` on the underlying data. +In the case of the category ``\mathbf{Vect}`` a valid braiding consists of just flipping the objects/morphisms involved, e.g. for a simple cartesian tensor, permuting the tensor indices is equivalent to applying Julia's function `permutedims` on the underlying data. Less trivial braiding implementations arise in the context of tensors with symmetries (where the fusion tree needs to be reordered, as discussed in [Sectors, representation spaces and fusion trees](@ref s_sectorsrepfusion)) or in the case of ``\mathbf{SVect}``, which will again be studied in detail at the end of this section. The braiding of a space and a dual space also follows naturally, it is given by ``τ_{V^*,W} = λ_{W ⊗ V^*} ∘ (ϵ_V ⊗ \mathrm{id}_{W ⊗ V^*}) ∘ (\mathrm{id}_{V^*} ⊗ τ_{V,W}^{-1} ⊗ \mathrm{id}_{V^*}) ∘ (\mathrm{id}_{V^*⊗ W} ⊗ η_V) ∘ ρ_{V^* ⊗ W}^{-1}``, i.e. @@ -321,7 +321,7 @@ This construction of the pivotal structure can than be used to define the trace, ``` Note finally, that a ribbon category where the braiding is symmetric, is known as a **compact closed category**. -For a symmetric braiding, the trivial twist ``θ_V = \mathrm{id}_V`` is always a valid choice, but it might not be the choice that one necessarily want to use. +For a symmetric braiding, the trivial twist ``θ_V = \mathrm{id}_V`` is always a valid choice, but it might not be the choice that one necessarily wants to use. Let us study the case of ``\mathbf{SVect}`` again. Reinvoking our basis ``|m⟩ ∈ V`` and ``|n⟩ ∈ W``, the braiding ``τ_{V,W}`` is given by the Koszul sign rule, i.e. ``τ_{V,W}:|m⟩ ⊗_\mathrm{g} |n⟩ ↦ (-1)^{|m| |n|} |n⟩ ⊗_\mathrm{g} |m⟩``. Hence, braiding amounts to flipping the two spaces, but picks up an additional minus sign if both ``|m⟩ ∈ V_1`` and ``|n⟩ ∈ W_1``. @@ -356,7 +356,7 @@ In the graphical representation, the dagger of a morphism can be represented by where for completeness we have also depicted the graphical representation of the transpose, which is a very different operation. In particular, the dagger does not reverse the order of the tensor product. -Note that, for readibility, we have not mirrored or rotated the label in the box, but this implies that we need to use a type of box for which the action of mirroring or rotating can be observed. +Note that, for readability, we have not mirrored or rotated the label in the box, but this implies that we need to use a type of box for which the action of mirroring or rotating can be observed. A dagger monoidal category is one in which the associator and left and right unitor are unitary morphisms. Similarly, a dagger braided category also has a unitary braiding, and a dagger balanced category in addition has a unitary twist. @@ -391,7 +391,7 @@ The bosonic version is obtained by restricting to the subcategory ``\mathbf{Vect ## [Direct sums, simple objects and fusion categories](@id ss_fusion) -These last two section on fusion categories is also applicable, in a straightforward manner, to ``\mathbf{Vect}`` and ``\mathbf{SVect}``, but is rather meant to provide the background of working with symmetries. +These last two sections on fusion categories are also applicable, in a straightforward manner, to ``\mathbf{Vect}`` and ``\mathbf{SVect}``, but are rather meant to provide the background of working with symmetries. We first need two new concepts: * An object ``W ∈ \mathrm{Ob}(C)`` is a **direct sum** of objects ``V_1, V_2, …, V_k ∈ \mathrm{Ob}(C)`` if there exists a family morphisms ``x_α ∈ \mathrm{Hom}(V_α,W)`` and ``y^α ∈ \mathrm{Hom}(W,V_α)`` such that ``\mathrm{id}_W = ∑_{α=1}^{k} x_α ∘ y^α`` and ``y^α ∘ x_β = δ^α_β \mathrm{id}_{V_α}``. @@ -415,7 +415,7 @@ This property is always satisfied for a **pre-fusion category** ``C``, i.e. a mo * every object ``V ∈ \mathrm{Ob}(C)`` can be written as a direct sum of a finite family of elements from ``\mathcal{S}``. Note that in the direct sum decomposition of an object ``V``, a particular simple object ``V_i`` might appear multiple times. -This number is known as the multiplicity index ``N^V_i``, and equal to the rank of ``\mathrm{Hom}(V,V_i)`` or, equivalently, of ``\mathrm{Hom}(V_i,V)``. +This number is known as the multiplicity index ``N^V_i``, and is equal to the rank of ``\mathrm{Hom}(V,V_i)`` or, equivalently, of ``\mathrm{Hom}(V_i,V)``. Hence, we can choose inclusion and projection maps ``x_{i,μ}:V_i→V`` and ``y^{i,μ}:V→V_i`` for ``μ = 1,\ldots, N^V_i``, such that ``\mathrm{id}_V = \sum_{i}\sum_{μ=1}^{N_V^i} x_{i,μ} ∘ y^{i,μ}`` and ``y^{i,μ} ∘ x_{j,ν} = δ^i_j δ^μ_ν``. In particular, for a simple object ``V``, it either appears in ``\mathcal{S}`` or is isomorphic to an object ``S``. We thus have ``N^V_i = 1`` for one particular object ``V_i`` and ``N^V_j= 0`` for all other ``j``, with ``x_{i}`` and ``y^i = (x_i)^{-1}`` representing the isomorphism between ``V`` and ``V_i``. @@ -436,8 +436,7 @@ Before continuing, let us use some examples to sketch the relevance of the conce As mentioned, the categories ``\mathbf{Vect}_𝕜`` and ``\mathbf{SVect}_𝕜`` have ``I ≂ 𝕜`` as simple object. For ``\mathbf{Vect}``, this is the only simple object, i.e. any other vector space ``V`` over ``𝕜`` can be thought of as a direct sum over ``N^V_I = \mathrm{dim}(V)`` multiple copies of ``𝕜``. In ``\mathbf{SVect}``, the object ``J = 0 ⊕ 𝕜`` with ``J_0=0`` the zero dimensional space and ``J_1 ≂ 𝕜`` is another simple object. -Clearly, there are no non-zero grading preserving morphisms between ``I`` and ``J``, i. -. ``\mathrm{Hom}(I,J) = 0``, whereas ``\mathrm{Hom}(J,J) ≂ 𝕜``. Any other super vector space ``V=V_0 ⊕ V_1`` can be written as a direct sum over ``N^V_I = \mathrm{dim}(V_0)`` copies of ``I`` and ``N^V_J = \mathrm{dim}(V_1)`` copies of ``J``. +Clearly, there are no non-zero grading preserving morphisms between ``I`` and ``J``, i.e. ``\mathrm{Hom}(I,J) = 0``, whereas ``\mathrm{Hom}(J,J) ≂ 𝕜``. Any other super vector space ``V=V_0 ⊕ V_1`` can be written as a direct sum over ``N^V_I = \mathrm{dim}(V_0)`` copies of ``I`` and ``N^V_J = \mathrm{dim}(V_1)`` copies of ``J``. A more representative example is that of the category ``C = \mathbf{Rep}_{\mathsf{G}}``, the category of representations of a group ``\mathsf{G}``. Colloquially, this could be thought of as a subcategory of ``\mathbf{Vect}`` containing as objects vector spaces ``V`` on which a representation of ``\mathsf{G}`` is defined, denoted as ``u_V(g)`` for ``g ∈ \mathsf{G}``, and as morphisms the equivariant transformations, i.e. intertwiners between the representations on the source and target: @@ -460,7 +459,7 @@ A braided pivotal fusion category is spherical if and only if it is a ribbon cat ## [Topological data of a unitary pivotal fusion category](@id ss_topologicalfusion) -More explicitly, the different structures (monoidal structure, duals and pivotal structure, braiding and twists) in a fusion category can be characterized in terms of the simple objects, which we will henceforth denoted with just ``a`` instead of ``V_a``. +More explicitly, the different structures (monoidal structure, duals and pivotal structure, braiding and twists) in a fusion category can be characterized in terms of the simple objects, which we will henceforth denote with just ``a`` instead of ``V_a``. This gives rise to what is known as the *topological data* of a unitary pivotal fusion category, most importantly the ``N``, ``F`` and ``R`` symbols, which are introduced in this final section. ### Monoidal structure @@ -531,7 +530,7 @@ In a strict category, or in the graphical notation, the associator ``α`` is omi The matrix ``F^{abc}_d`` is thus a unitary matrix. The pentagon coherence equation can also be rewritten in terms of these matrix elements, and as such yields the celebrated pentagon equation for the F-symbols. In a similar fashion, the unitors result in ``N^{a1}_{b} = N^{1a}_b = δ^{a}_b`` (where we have now written ``1`` instead of ``I`` for the unit object) and the triangle equation leads to additional relations between the F- symbols involving the unit object. -In particular, if we identify ``X^{1a}_{a,1}:a→(1⊗a)`` with ``λ_a^†`` and ``X^{a1}_{a,1}:a→(a⊗1)`` with ``ρ_a^†``, the triangle equation and its collaries imply that ``[F^{1ab}_{c}]_{(11μ)}^{(cν1)} = δ^{ν}_{μ}``, and similar relations for ``F^{a1b}_c`` and ``F^{ab1}_c``, which are graphically represented as +In particular, if we identify ``X^{1a}_{a,1}:a→(1⊗a)`` with ``λ_a^†`` and ``X^{a1}_{a,1}:a→(a⊗1)`` with ``ρ_a^†``, the triangle equation and its corollaries imply that ``[F^{1ab}_{c}]_{(11μ)}^{(cν1)} = δ^{ν}_{μ}``, and similar relations for ``F^{a1b}_c`` and ``F^{ab1}_c``, which are graphically represented as ```@raw html Fmove1 @@ -556,7 +555,7 @@ In the context of tensors in quantum physics, we would like to be able to repres This means that ``\mathrm{Hom}(a^*,\bar{a})`` is isomorphic to ``𝕜`` and contains a single linearly independent element, ``Z_a``, which is a unitary isomorphism such that ``Z_a^\dagger ∘ Z_a = \mathrm{id}_{a^*}`` and ``Z_a ∘ Z_a^\dagger = \mathrm{id}_{\bar{a}}``. Using the transpose, we obtain ``Z_a^* ∈ \mathrm{Hom}(\bar{a}^*,a)``, and thus it is proportional to ``Z_{\bar{a}}``, i.e. ``Z_a^* = χ_a Z_{\bar{a}}`` with ``χ_a`` a complex phase (assuming ``𝕜 = ℂ``). Another transpose results in ``Z_{\bar{a}}^* = χ_{\bar{a}} Z_a`` with ``χ_{\bar{a}} = \overline{χ_{a}}``, where bar of a scalar quantity denotes its complex conjugate to avoid confusion with the transpose functor. -If ``a``and ``\bar{a}`` are distinct, we can essentially choose ``Z_{\bar{a}}`` such that ``χ_a`` is ``1``. +If ``a`` and ``\bar{a}`` are distinct, we can essentially choose ``Z_{\bar{a}}`` such that ``χ_a`` is ``1``. However, for ``a=\bar{a}``, the value of ``χ_a`` cannot be changed, but must satisfy ``χ_a^2 = 1``, or thus ``χ_a = ±1``. This value is a topological invariant known as the *Frobenius-Schur indicator*. Graphically, we represent this isomorphism and its relations as @@ -576,7 +575,7 @@ Here, we have denoted ``d_a = \mathrm{dim}(a) = \mathrm{tr}(\mathrm{id}_a)`` for With this information, we can then compute ``[F^{a\bar{a}a}_a]``, which has a single element (it's a ``1 × 1`` matrix), and find ``[F^{a\bar{a}a}_a] = \frac{χ_a}{d_a}``, where we've used ``\tilde{η}_a = ϵ_a^†`` and the snake rules. Hence, both the quantum dimensions and the Frobenius-Schur indicator are encoded in the F-symbol. Hence, they do not represent new independent data. -Again, the graphical representation is more enlightning: +Again, the graphical representation is more enlightening: ```@raw html ZtoF @@ -673,7 +672,7 @@ If ``a = \bar{a}``, we can furthermore relate the twist, the braiding and the Fr For the recurring example of ``\mathbf{Rep}_{\mathsf{G}}``, the braiding acts simply as the swap of the two vector spaces on which the representations are acting and is thus symmetric, i.e. ``τ_{b,a} ∘ τ_{a,b} = \mathrm{id}_{a⊗b}``. All the twists are simply ``θ_a = 1``. For an irrep that is self-dual, i.e. ``\bar{a}=a``, the final expression simplifies to ``R^{aa}_1 = χ_a`` and thus states that the fusion from ``a ⊗ a`` to the trivial sector is either symmetric under swaps if ``χ_a=1`` or antisymmetric if ``χ_a=-1``. -For the case of ``\mathsf{SU}_2``, the coupling of two spin ``j`` states to a singlet it symmetric for integer ``j`` and odd for half-integer ``j``. +For the case of ``\mathsf{SU}_2``, the coupling of two spin ``j`` states to a singlet is symmetric for integer ``j`` and odd for half-integer ``j``. With this, we conclude our exposition of unitary fusion categories. There are many fusion categories that do not originate from the representation theory of groups, but are related to quantum groups and the representation theory of quasi-triangular Hopf algebras. diff --git a/docs/src/appendix/symmetric_tutorial.md b/docs/src/appendix/symmetric_tutorial.md index de14dc28f..5aa4cc923 100644 --- a/docs/src/appendix/symmetric_tutorial.md +++ b/docs/src/appendix/symmetric_tutorial.md @@ -95,7 +95,7 @@ Representing these operators as `TensorMap`s, the invariance of ``H`` under a gl
ZZX_symm
``` -These identitities precisely mean that these local tensors transform trivially under a tensor product representation of ``\mathbb{Z}_2``. +These identities precisely mean that these local tensors transform trivially under a tensor product representation of ``\mathbb{Z}_2``. This implies that, recalling [the introduction on symmetries](@ref ss_symmetries), in an appropriate basis for the local physical vector space, our local tensors would become block-diagonal where each so-called *matrix block* is labeled by a ``\mathbb{Z}_2`` irrep. The appropriate local basis transformation is precisely the one that brings the local representation ``X`` into block-diagonal form. Clearly, this transformation is nothing more than the Hadamard transformation which maps the computational basis of ``Z`` eigenstates ``\{\ket{\uparrow}, \ket{\downarrow}\}`` to that of the ``X`` eigenstates ``\{\ket{+}, \ket{-}\}`` defined as ``\ket{+} = \frac{\ket{\uparrow} + \ket{\downarrow}}{\sqrt{2}}`` and ``\ket{-} = \frac{\ket{\uparrow} - \ket{\downarrow}}{\sqrt{2}}``. @@ -130,7 +130,7 @@ We will return to the implications of irreps with *non-Abelian* fusion rules [la !!! note Within TensorKit.jl, the nature of the fusion rules for charges of a given symmetry are represented by the [`FusionStyle`](@ref) of the corresponding `Sector` subtype. What we refer to as "Abelian" fusion rules in this tutorial corresponds to `UniqueFusion <: FusionStyle`. - We will also consider [examples](@ref ss_non_abelian) of two different kinds of non-Abelian" fusion rules, corresponding to `MultipleFusion <: FusionStyle` styles. + We will also consider [examples](@ref ss_non_abelian) of two different kinds of non-Abelian fusion rules, corresponding to `MultipleFusion <: FusionStyle` styles. For the case of the ``\mathbb{Z}_2`` irreps, the fusion rules are Abelian, and are given by addition modulo 2, ```math @@ -167,7 +167,7 @@ In TensorKit.jl, these reduced tensor elements corresponding to the fusion trees This view of the underlying symmetry structure in terms of fusion trees of irreps and corresponding reduced tensor elements is a very convenient way of working with the `TensorMap` type. In fact, this symmetry structure is inherently ingrained in a `TensorMap`, and goes beyond the group-like symmetries we have considered until now. -In this more general setting, we will refer to the labels that appear on this fusion trees as *charges* or *sectors*. +In this more general setting, we will refer to the labels that appear on these fusion trees as *charges* or *sectors*. These can be thought of as generalization of group irreps, and appear in the context of TensorKit.jl as instances of the [`Sector`](@ref) type. Consider a generic fusion tree of the form @@ -190,7 +190,7 @@ You will find such a `FusionTree` has the following properties encoded into its - `vertices::NTuple{L,T}`: list of fusion vertex labels of type `T` and length `L = N - 1` For our current application only `uncoupled` and `coupled` are relevant, since ``\mathbb{Z}_2`` irreps are self-dual and have Abelian fusion rules, so that irreps on the inner lines of a fusion tree are completely determined by the uncoupled irreps. -We will come back to these other properties when discussion more involved applications. +We will come back to these other properties when discussing more involved applications. Given some `TensorMap`, the method [`fusiontrees`](@ref) returns an iterator over all pairs of splitting and fusion trees that label the subblocks of `t`. ### Constructing a ``\mathbb{Z}_2``-symmetric `TensorMap` @@ -247,7 +247,7 @@ To assess the underlying structure of a symmetric tensor, it is often useful to subblocks(ZZ) ``` While all entries are zero, we see that all eight valid fusion trees with two incoming irreps and two outgoing irreps [of the type above](fusiontree) are listed with their corresponding subblock data. -Each of these subblocks is an array of shape ``(1, 1, 1, 1)`` since each irrep occuring in the space ``V`` has degeneracy 1. +Each of these subblocks is an array of shape ``(1, 1, 1, 1)`` since each irrep occurring in the space ``V`` has degeneracy 1. Using the [`fusiontrees`](@ref) method and the fact that we can index a `TensorMap` using a splitting/fusion tree pair, we can now fill in the nonzero subblocks of the operator by observing that the ``ZZ`` operator flips the irreps of the uncoupled charges in the domain with respect to the codomain, as shown in the diagrams above. Flipping a given `Z2Irrep` in the codomain can be implemented by fusing them with the sign irrep `Z2Irrep(1)`, giving: @@ -322,7 +322,7 @@ This Hamiltonian is invariant under conjugation by the global particle number op U = \sum_i N_i ``` This invariance corresponds to a ``\mathsf{U}_1`` particle number symmetry, which can again be manifestly imposed when constructing the Hamiltonian terms as `TensorMap`s. -From the representation theory of ``\mathsf{U}_1``, we know that its irreps are all one-dimensional and can be labeled by integers ``n`` where the tensor product of two irreps is corresponds to addition of these labels, giving the Abelian fusion rules +From the representation theory of ``\mathsf{U}_1``, we know that its irreps are all one-dimensional and can be labeled by integers ``n`` where the tensor product of two irreps corresponds to addition of these labels, giving the Abelian fusion rules ```math n_1 \otimes n_2 \cong (n_1 + n_2). @@ -401,8 +401,8 @@ However, it is possible to construct them as `TensorMap`s using an *auxiliary ve The creation operator ``a^+`` violates particle number conservation by mapping the occupation number ``n`` to ``n + 1``. From the point of view of representation theory, this process can be thought of as the *fusion* of an `U1Irrep(n)` with an `U1Irrep(1)`, naturally giving the fusion product `U1Irrep(n + 1)`. This means we can represent ``a^+`` as a `TensorMap(..., V ← V ⊗ A)`, where the auxiliary vector space `A` contains the ``+1`` irrep with degeneracy 1, `A = U1Space(1 => 1)`. -Similarly, the decrease in occupation number when acting with ``a^-`` can be thought of as the *splitting* of an `U1Irrep(n)` into an `U1Irrep(n - 1)` and an `U1Irrep(1)`, leading to a representation in terms of a `TensorMap(. -., A ⊗ V ← V)`. Based on these observations, we can represent the matrix elements \eqref{eq:bosonopmatel} as subblocks labeled by the ``\mathsf{U}_1`` fusion trees +Similarly, the decrease in occupation number when acting with ``a^-`` can be thought of as the *splitting* of an `U1Irrep(n)` into an `U1Irrep(n - 1)` and an `U1Irrep(1)`, leading to a representation in terms of a `TensorMap(..., A ⊗ V ← V)`. +Based on these observations, we can represent the matrix elements \eqref{eq:bosonopmatel} as subblocks labeled by the ``\mathsf{U}_1`` fusion trees. ```@raw html
bosonops
@@ -457,7 +457,7 @@ It is then simple to check that this is indeed what we expect. !!! note From the construction of the Hamiltonian operators [in terms of creation and annihilation operators](bosonham) we clearly see that they are invariant under a transformation ``a^\pm \to e^{\pm i\theta} a^\pm``. More generally, for a two-site operator that is defined as the contraction of two one-site operators across an auxiliary space, modifying the one-site operators by applying transformations ``Q`` and ``Q^{-1}`` on their respective auxiliary spaces for any invertible ``Q`` leaves the resulting contraction unchanged. - This ambiguity in the definition clearly shows that one should really always think in terms of the fully symmetric procucts of ``a^+`` and ``a^-`` rather than in terms of these operators themselves. + This ambiguity in the definition clearly shows that one should really always think in terms of the fully symmetric products of ``a^+`` and ``a^-`` rather than in terms of these operators themselves. In particular, one can always decompose such a symmetric product into the [form above](bosonham) by means of an SVD. @@ -502,7 +502,7 @@ The fusion rules of these irreps are the same as for ``\mathbb{Z}_2``. Similar to the previous case, the local symmetry operator ``Q_i`` is already diagonal, so the occupation number basis coincides with the irrep basis and we don't need an additional basis transform. The important difference with a regular ``\mathbb{Z}_2`` symmetry is that the irreps of ``f\mathbb{Z}_2`` have fermionic braiding statistics, in the sense that exchanging two odd irreps gives rise to a minus sign. -In TensorKit.jl, an ``f\mathbb{Z}_2``-graded vector spaces is represented as a `Vect[FermionParity]` space, where a given ``f\mathbb{Z}_2`` irrep can be represented as a [`FermionParity`](@ref FermionParity) sector instance. +In TensorKit.jl, an ``f\mathbb{Z}_2``-graded vector space is represented as a `Vect[FermionParity]` space, where a given ``f\mathbb{Z}_2`` irrep can be represented as a [`FermionParity`](@ref FermionParity) sector instance. Using the simplest instance of a vector space containing a single even and odd irrep, we can already demonstrate the corresponding fermionic braiding behavior by [performing a permutation](@ref TensorKit.permute) on a simple `TensorMap`. ```@example symmetric_tutorial @@ -523,7 +523,7 @@ In other words, when exchanging the two domain vector spaces, the reduced tensor We can directly construct the Hamiltonian terms as symmetric `TensorMap`s using the same procedure as before starting from their matrix elements in the occupation number basis. However, in this case we should be a bit more careful about the precise definition of the basis states in composite systems. Indeed, the tensor product structure of fermionic systems is inherently tricky to deal with, and should ideally be treated in the context of [*super vector spaces*](https://en.wikipedia.org/wiki/Super_vector_space). -For two sites, we can define the following basis states on top of the fermionic vacuuum ``\ket{00}``: +For two sites, we can define the following basis states on top of the fermionic vacuum ``\ket{00}``: ```math \begin{align*} \ket{01} &= c_2^+ \ket{00}, \\ @@ -611,7 +611,7 @@ end subblocks(N) ``` -We can easily all the reduced tensor elements are indeed correct. +We can easily check that all the reduced tensor elements are indeed correct. !!! note Working with fermionic systems is inherently tricky, as can already be seen from something as simple as computing matrix elements of fermionic operators. @@ -634,7 +634,7 @@ In particular, this means that the conversion of an operator, given its matrix e Therefore, we require some more discussion before we can actually move on to an example. We'll start by discussing the general structure of a `TensorMap` which is symmetric under a non-Abelian group symmetry. -We then given an example based on ``\mathsf{SU}_2``, where we construct the Heisenberg Hamiltonian using two different approaches. +We then give an example based on ``\mathsf{SU}_2``, where we construct the Heisenberg Hamiltonian using two different approaches. Finally, we show how the more intuitive approach can be used to obtain an elegant generalization to the ``\mathsf{SU}_N``-symmetric case. @@ -736,8 +736,7 @@ This procedure works for any group symmetry, and all we need are matrix elements In the following, we demonstrate this explicit procedure for the particular example of ``G = \mathsf{SU}_2``. However, it should be noted that, for other non-Abelian groups, the Clebsch-Gordan coefficients may not be as easy to compute (generically, no closed formulas exist). In addition, the procedure for manually projecting out the reduced tensor elements requires being particularly careful about the correspondence between the basis states used to define the original matrix elements and those implied by the Clebsch-Gordan coefficients. -Finally, for some symmetries supported in TensorKit. -l, there are simply no Clebsch-Gordan coefficients. +Finally, for some symmetries supported in TensorKit.jl, there are simply no Clebsch-Gordan coefficients. Therefore, it is often easier and sometimes simply necessary to directly construct the symmetric tensor and then fill in its reduced tensor elements based on some representation theory. We will cover some examples of this below. @@ -752,8 +751,8 @@ f = fusiontensor(SU2Irrep(1//2), SU2Irrep(1//2), SU2Irrep(1)) ``` We see that this fusion tensor has a size `2×2×3×1`, which contains an additional trailing `1` to what we might expect. -In the general case, `fusiontensor` returns a four-dimensional array, where the size of the first three dimensions corresponds to the dimensions of the irrep spaces under consideration, and the last index lables the different fusion channels, where its dimension corresponds to the number of distinct ways the irreps ``l_1`` and ``l_2`` can fuse to irrep ``k``. -This is precicely the extra label of the Clebsch-Gordan coefficients that is required in the the presence of fusion multiplicities. +In the general case, `fusiontensor` returns a four-dimensional array, where the size of the first three dimensions corresponds to the dimensions of the irrep spaces under consideration, and the last index labels the different fusion channels, where its dimension corresponds to the number of distinct ways the irreps ``l_1`` and ``l_2`` can fuse to irrep ``k``. +This is precisely the extra label of the Clebsch-Gordan coefficients that is required in the presence of fusion multiplicities. Since ``\mathsf{SU}_2`` is multiplicity-free, we can just discard this last index here. We can now explicitly verify that this `fusiontensor` indeed does what we expect it to do: @@ -873,7 +872,7 @@ First, we rewrite the exchange interaction in the following way: \vec{S}_i \cdot \vec{S}_j = \frac{1}{2} \left( \left( \vec{S}_i + \vec{S}_j \right)^2 - \vec{S}_i^2 - \vec{S}_j^2 \right) \end{equation} ``` -Here, ``\vec{S}_i`` and ``\vec{S}_j`` are spin operators on the physical irrep, while total spin operator ``\vec{S}_i + \vec{S}_j`` can be decomposed onto the different coupled irreps ``k``. +Here, ``\vec{S}_i`` and ``\vec{S}_j`` are spin operators on the physical irrep, while the total spin operator ``\vec{S}_i + \vec{S}_j`` can be decomposed onto the different coupled irreps ``k``. It is a well known fact that the quadratic sum of the generators of ``\mathsf{SU}_2``, often refered to as the [*quadratic Casimir*](https://en.wikipedia.org/wiki/Representation_theory_of_SU(2)#The_Casimir_element), commutes with all generators. By [Schur's lemma](https://en.wikipedia.org/wiki/Schur%27s_lemma), it must then act proportionally to the identity on every irrep, where the corresponding eigenvalue is determined by the spin irrep label. In particular, we have for each irrep ``l`` @@ -909,7 +908,7 @@ We end this subsection with some comments on the generalization of the above dis As foreshadowed above, the irreps of ``\mathsf{SU}_N`` in general have an even more complicated structure. In particular, they can admit so-called *fusion multiplicities*, where the fusion of two irreps can have not only multiple distinct outcomes, but they can even fuse to a given irrep in multiple inequivalent ways. We can demonstrate this behavior for the adjoint representation of ``\mathsf{SU}_3``. -For this we can use the the [SUNRepresentations.jl](https://github.com/QuantumKitHub/SUNRepresentations.jl) package which provides an interface for working with irreps of ``\mathsf{SU}_N`` and their Clebsch-Gordan coefficients. +For this we can use the [SUNRepresentations.jl](https://github.com/QuantumKitHub/SUNRepresentations.jl) package which provides an interface for working with irreps of ``\mathsf{SU}_N`` and their Clebsch-Gordan coefficients. A particular representation is represented by an `SUNIrrep{N}` which can be used with TensorKit.jl. The eight-dimensional adjoint representation of ``\mathsf{SU}_3`` is given by ```@setup symmetric_tutorial @@ -947,7 +946,7 @@ For any ``N``, the [quadratic Casimir](https://en.wikipedia.org/wiki/Casimir_ele \Omega = \sum_k T^k T^k ``` commutes with all ``\mathsf{SU}_N`` generators, meaning it has a well defined eigenvalue in each irrep. -This observation then immediately given the reduced tensor elements of the exchange interaction as +This observation then immediately gives the reduced tensor elements of the exchange interaction as ```@raw html
SUN_fusiontrees
``` @@ -981,7 +980,7 @@ f.vertices ``` !!! note - While we have given an explicit example using ``\mathsf{SU}_3`` with the adoint irrep on the physical level, the same construction holds for the general ``\mathsf{SU}_N`` with arbitrary physical irreps. + While we have given an explicit example using ``\mathsf{SU}_3`` with the adjoint irrep on the physical level, the same construction holds for the general ``\mathsf{SU}_N`` with arbitrary physical irreps. All we require is the expression for the eigenvalues of the quadratic Casimir in each irrep. diff --git a/docs/src/lib/fusiontrees.md b/docs/src/lib/fusiontrees.md index 8e037af93..57033eca6 100644 --- a/docs/src/lib/fusiontrees.md +++ b/docs/src/lib/fusiontrees.md @@ -2,6 +2,7 @@ ```@meta CurrentModule = TensorKit +CollapsedDocStrings = true ``` # Type hierarchy diff --git a/docs/src/lib/sectors.md b/docs/src/lib/sectors.md index f56980bf0..de3b1b702 100644 --- a/docs/src/lib/sectors.md +++ b/docs/src/lib/sectors.md @@ -2,6 +2,7 @@ ```@meta CurrentModule = TensorKit +CollapsedDocStrings = true ``` ## Type hierarchy @@ -126,7 +127,7 @@ Because we sometimes want to customize the string representation of our sector t TensorKitSectors.type_repr ``` -Finally, we provide functionality to compile all revelant methods for a sector: +Finally, we provide functionality to compile all relevant methods for a sector: ```@docs TensorKitSectors.precompile_sector diff --git a/docs/src/lib/spaces.md b/docs/src/lib/spaces.md index e5705fe3e..fbeaacb6d 100644 --- a/docs/src/lib/spaces.md +++ b/docs/src/lib/spaces.md @@ -2,6 +2,7 @@ ```@meta CurrentModule = TensorKit +CollapsedDocStrings = true ``` ## Type hierarchy @@ -58,7 +59,7 @@ const SU₂Space = SU2Space ## [Methods](@id s_spacemethods) -Methods often apply similar to e.g. spaces and corresponding tensors or tensor maps, e.g.: +Methods often apply similarly to e.g. spaces and corresponding tensors or tensor maps, e.g.: ```@docs field diff --git a/docs/src/lib/tensors.md b/docs/src/lib/tensors.md index 1ef3da0fb..217c6640a 100644 --- a/docs/src/lib/tensors.md +++ b/docs/src/lib/tensors.md @@ -147,17 +147,17 @@ Random.randexp! The operations that can be performed on an `AbstractTensorMap` can be organized into the following categories: -* *vector operations*: these do not change the `space` or index structure of a tensor and can be straightforwardly implemented on on the full data. +* *vector operations*: these do not change the `space` or index structure of a tensor and can be straightforwardly implemented on the full data. All the methods described in [VectorInterface.jl](https://github.com/Jutho/VectorInterface.jl) are supported. For compatibility reasons, we also provide implementations for equivalent methods from LinearAlgebra.jl, such as `axpy!`, `axpby!`. -* *index manipulations*: these change (permute) the index structure of a tensor, which affects the data in a way that is fully determined by the categorical data of the `sectortype` of the tensor . +* *index manipulations*: these change (permute) the index structure of a tensor, which affects the data in a way that is fully determined by the categorical data of the `sectortype` of the tensor. * *(planar) contractions* and *(planar) traces* (i.e., contractions with identity tensors). Tensor contractions correspond to a combination of some index manipulations followed by a composition or multiplication of the tensors in their role as linear maps. Tensor contractions are however of such importance and frequency that they require a dedicated implementation. -* *tensor factorizations*, which relies on their identification of tensors with linear maps between tensor spaces. +* *tensor factorizations*, which rely on their identification of tensors with linear maps between tensor spaces. The factorizations are applied as ordinary matrix factorizations to the matrix blocks associated with the coupled charges. ### Index manipulations diff --git a/docs/src/man/contractions.md b/docs/src/man/contractions.md new file mode 100644 index 000000000..cc916556b --- /dev/null +++ b/docs/src/man/contractions.md @@ -0,0 +1,319 @@ +# [Tensor contractions and tensor networks](@id ss_tensor_contraction) + +One of the most important operations with tensor maps is to compose them, more generally known as contracting them. +As mentioned in the section on [category theory](@ref s_categories), a typical composition of maps in a ribbon category can graphically be represented as a planar arrangement of the morphisms (i.e. tensor maps, boxes with lines emanating from top and bottom, corresponding to source and target, i.e. domain and codomain), where the lines connecting the source and targets of the different morphisms should be thought of as ribbons, that can braid over or underneath each other, and that can twist. +Technically, we can embed this diagram in ``ℝ × [0,1]`` and attach all the unconnected line endings corresponding to objects in the source at some position ``(x,0)`` for ``x∈ℝ``, and all line endings corresponding to objects in the target at some position ``(x,1)``. +The resulting morphism is then invariant under what is known as *framed three-dimensional isotopy*, i.e. three-dimensional rearrangements of the morphism that respect the rules of boxes connected by ribbons whose open endings are kept fixed. +Such a two-dimensional diagram cannot easily be encoded in a single line of code. + +However, things simplify when the braiding is symmetric (such that over- and under- crossings become equivalent, i.e. just crossings), and when twists, i.e. self-crossings in this case, are trivial. +This amounts to `BraidingStyle(I) == Bosonic()` in the language of TensorKit.jl, and is true for any subcategory of ``\mathbf{Vect}``, i.e. ordinary tensors, possibly with some symmetry constraint. +The case of ``\mathbf{SVect}`` and its subcategories, and more general categories, are discussed below. + +In the case of trivial twists, we can deform the diagram such that we first combine every morphism with a number of coevaluations ``η`` so as to represent it as a tensor, i.e. with a trivial domain. +We can then rearrange the morphism to be all lined up horizontally, where the original morphism compositions are now being performed by evaluations ``ϵ``. +This process will generate a number of crossings and twists, where the latter can be omitted because they act trivially. +Similarly, double crossings can also be omitted. +As a consequence, the diagram, or the morphism it represents, is completely specified by the tensors it is composed of, and which indices between the different tensors are connected, via the evaluation ``ϵ``, and which indices make up the source and target of the resulting morphism. +If we also compose the resulting morphisms with coevaluations so that it has a trivial domain, we just have one type of unconnected lines, henceforth called open indices. +We sketch such a rearrangement in the following picture + +```@raw html +tensor unitary +``` + +Hence, we can now specify such a tensor diagram, henceforth called a tensor contraction or also tensor network, using a one-dimensional syntax that mimics [abstract index notation](https://en.wikipedia.org/wiki/Abstract_index_notation) and specifies which indices are connected by the evaluation map using Einstein's summation convention. +Indeed, for `BraidingStyle(I) == Bosonic()`, such a tensor contraction can take the same format as if all tensors were just multi-dimensional arrays. +For this, we rely on the interface provided by the package [TensorOperations.jl](https://github.com/QuantumKitHub/TensorOperations.jl). + +The above picture would be encoded as +```julia +@tensor E[a, b, c, d, e] := A[v, w, d, x] * B[y, z, c, x] * C[v, e, y, b] * D[a, w, z] +``` +or +```julia +@tensor E[:] := A[1, 2, -4, 3] * B[4, 5, -3, 3] * C[1, -5, 4, -2] * D[-1, 2, 5] +``` +where the latter syntax is known as NCON-style, and labels the unconnected or outgoing indices with negative integers, and the contracted indices with positive integers. + +A number of remarks are in order. +TensorOperations.jl accepts both integers and any valid variable name as dummy label for indices, and everything in between `[ ]` is not resolved in the current context but interpreted as a dummy label. +Here, we label the indices of a `TensorMap`, like `A::TensorMap{T, S, N₁, N₂}`, in a linear fashion, where the first position corresponds to the first space in `codomain(A)`, and so forth, up to position `N₁`. +Index `N₁ + 1` then corresponds to the first space in `domain(A)`. +However, because we have applied the coevaluation ``η``, it actually corresponds to the corresponding dual space, in accordance with the interface of [`space(A, i)`](@ref), and as indicated by the dotted box around ``A`` in the above picture. +The same holds for the other tensor maps. +Note that our convention also requires that we braid indices that we brought from the domain to the codomain, and so this is only unambiguous for a symmetric braiding, where there is a unique way to permute the indices. + +With the current syntax, we create a new object `E` because we use the definition operator `:=`. +Furthermore, with the current syntax, it will be a `Tensor`, i.e. it will have a trivial domain, and correspond to the dotted box in the picture above, rather than the actual morphism `E`. +We can also directly define `E` with the correct codomain and domain by rather using +```julia +@tensor E[a b c;d e] := A[v, w, d, x] * B[y, z, c, x] * C[v, e, y, b] * D[a, w, z] +``` +or +```julia +@tensor E[(a, b, c);(d, e)] := A[v, w, d, x] * B[y, z, c, x] * C[v, e, y, b] * D[a, w, z] +``` +where the latter syntax can also be used when the codomain is empty. +When using the assignment operator `=`, the `TensorMap` `E` is assumed to exist and the contents will be written to the currently allocated memory. +Note that for existing tensors, both on the left hand side and right hand side, trying to specify the indices in the domain and the codomain separately using the above syntax, has no effect, as the bipartition of indices is already fixed by the existing object. +Hence, if `E` has been created by the previous line of code, all of the following lines are now equivalent +```julia +@tensor E[(a, b, c);(d, e)] = A[v, w, d, x] * B[y, z, c, x] * C[v, e, y, b] * D[a, w, z] +@tensor E[a, b, c, d, e] = A[v w d; x] * B[(y, z, c); (x, )] * C[v e y; b] * D[a, w, z] +@tensor E[a b; c d e] = A[v; w d x] * B[y, z, c, x] * C[v, e, y, b] * D[a w; z] +``` +and none of those will or can change the partition of the indices of `E` into its codomain and its domain. + +Two final remarks are in order. +Firstly, the order of the tensors appearing on the right hand side is irrelevant, as we can reorder them by using the allowed moves of the Penrose graphical calculus, which yields some crossings and a twist. +As the latter is trivial, it can be omitted, and we just use the same rules to evaluate the newly ordered tensor network. +For the particular case of matrix-matrix multiplication, which also captures more general settings by appropriately combining spaces into a single line, we indeed find + +```@raw html +tensor contraction reorder +``` + +or thus, the following two lines of code yield the same result +```julia +@tensor C[i, j] := B[i, k] * A[k, j] +@tensor C[i, j] := A[k, j] * B[i, k] +``` +Reordering of tensors can be used internally by the `@tensor` macro to evaluate the contraction in a more efficient manner. +In particular, the NCON-style of specifying the contraction gives the user control over the order, and there are other macros, such as `@tensoropt`, that try to automate this process. +There is also an `@ncon` macro and `ncon` function, and we recommend reading the [manual of TensorOperations.jl](https://quantumkithub.github.io/TensorOperations.jl/stable/) to learn more about the possibilities and how they work. + +A final remark involves the use of adjoints of tensors. +The current framework is such that the user should not be too worried about the actual bipartition into codomain and domain of a given `TensorMap` instance. +Indeed, for tensor contractions the `@tensor` macro figures out the correct manipulations automatically. +However, when wanting to use the `adjoint` of an instance `t::TensorMap{T, S, N₁, N₂}`, the resulting `adjoint(t)` is an `AbstractTensorMap{T, S, N₂, N₁}` and one needs to know the values of `N₁` and `N₂` to know exactly where the `i`th index of `t` will end up in `adjoint(t)`, and hence the index order of `t'`. +Within the `@tensor` macro, one can instead use `conj()` on the whole index expression so as to be able to use the original index ordering of `t`. +For example, for `TensorMap{T, S, 1, 1}` instances, this yields exactly the equivalence one expects, namely one between the following two expressions: + +```julia +@tensor C[i, j] := B'[i, k] * A[k, j] +@tensor C[i, j] := conj(B[k, i]) * A[k, j] +``` + +For e.g. an instance `A::TensorMap{T, S, 3, 2}`, the following two syntaxes have the same effect within an `@tensor` expression: `conj(A[a, b, c, d, e])` and `A'[d, e, a, b, c]`. + +## Fermionic tensor contractions + +Whenever `BraidingStyle(i) == Fermionic()`, some complications come up. +The most important distinction from the `Bosonic()` case is that twists are no longer trivial, such that we must be careful about how we can manipulate network diagrams. + +To illustrate these complications, we take a look at a concrete example first, and study the following tensor network: + +```@raw html +fermionic contraction example +``` + +```@example fermioncontraction +using TensorKit # hide +V₁ = Vect[FermionParity](0 => 1, 1 => 1) +V₂ = Vect[FermionParity](0 => 2, 1 => 2) +A = rand(V₁ ← V₁ ⊗ V₂) +X = rand(V₁ ← V₁) +B = rand(V₁ ⊗ V₂ ← V₁) +``` + +We can expand this into binary contractions, by first contracting `X` with `A`, and then contracting the result with `B`: + +```@example fermioncontraction +AX = repartition(A, 2, 1) * X +AXB = repartition(AX, 1, 2) * B +``` + +Alternatively, we could decide that we first wish to contract `A` with `B`, and only then contract the result with `X`: + +```@example fermioncontraction +AB = permute(A, ((1, 2), (3,))) * permute(B, ((2,), (1, 3))) +ABX = repartition(permute(AB, ((1, 4), (2, 3))) * repartition(X, 2, 0), 1, 1) +``` + +This is where the issue becomes clearer, as the results are no longer equal: + +```@example fermioncontraction +AXB ≈ ABX +``` + +### Trivializing the twist + +So what happened? +If we carefully inspect what we actually computed here, we can show that in order to deform one diagram into the other, we have to introduce a self-crossing, which then altered the result. +While the example here is still simple to follow, in general we would like that the result of `@tensor` expressions does not depend on the input order of the tensors. +This is especially true for larger expressions where we wish to dynamically compute the optimal contraction order, as this would alter the order in a very non-transparent manner. + +The way out of this effectively consists of absorbing this twist in the coevaluation map ``η``. +This modified map ``η′ := η ∘ θ`` where ``θ`` represents the twist ensures that the result no longer depends on the order of evaluation. +In particular, one can show that any time two tensors would swap places, we would simultaneously exchange one evaluation map ``ϵ`` for a coevaluation ``̃η``, while also incurring a twist ``θ`` such that both cancel out. +To make this concrete, we show how our previous example now leads to a unique result: + +```@example fermioncontraction +function fermion_mul(A, B) + return A * twist(B, findall(isdual, codomain(B).spaces)) +end + +# order I: +AX = fermion_mul(repartition(A, 2, 1), X) +AXB = fermion_mul(repartition(AX, 1, 2) , B) + +# order II: +AB = fermion_mul(permute(A, ((1, 2), (3,))), permute(B, ((2,), (1, 3)))) +ABX = repartition(fermion_mul(permute(AB, ((1, 4), (2, 3))), repartition(X, 2, 0)), 1, 1) + +AXB ≈ ABX +``` + +This is the so-called **supertrace** formalism, and is effectively what `@tensor` ends up implementing for fermionic contractions. +For more details about this formalism, we refer to [^Mortier]. + +```@example fermioncontraction +# @tensor +@tensor result[-1; -2] := A[-1; 1 2] * X[1; 3] * B[3 2; -2] + +AXB ≈ result +``` + +### (Non)-unitarity + +While this modified ``̃η`` solves the issues related to contractions, it does come at a cost. +The main issue is that this map does not constitute a positive definite map, and in particular is at odds with a positive inner product. +Such a positive inner product is however required to properly define (orthogonal) factorizations, non-negative norms, etc. + +Therefore, we reserve the supertrace formalism exclusively for tensor contractions. +For matrix-like operations such as factorizations, matrix functions, norms, etc, we retain the positive definite inner product. +It is also always possible to manually emulate one or the other, by inserting appropriate calls to `twist`. +In what follows, we simply showcase some noteworthy differences between the two formalisms, as these can be a common source of errors. +Throughout, we use the following simple fermionic tensor as a running example: + +```@example fermionnorm +using TensorKit # hide +V = Vect[FermionParity](0 => 1, 1 => 1) +t = ones(V' ← V') +``` + +- Computing a norm via a contraction: + the squared norm of `t`, computed via the supertrace contraction, no longer agrees with `norm(t)^2`. + In particular, the `@tensor` self-contraction can even vanish for a manifestly non-zero tensor: + +```@example fermionnorm +norm(t)^2, @tensor conj(t[a; b]) * t[a; b] +``` + +Inserting a `twist` on the contracted codomain index cancels the twist that `@tensor` automatically introduces, and recovers the trace-formalism result: + +```@example fermionnorm +norm(t)^2 ≈ @tensor conj(t[a; b]) * twist(t, 1)[a; b] +``` + +- Using unitarity to simplify `U * U' ≈ I`: + the factor `U` returned by `svd_compact` is left-isometric in the *trace* sense, i.e. `U' * U ≈ id(domain(U))` as a matrix product, but this identity no longer holds when the same product is written as a tensor contraction: + +```@example fermionnorm +U, S, Vᴴ = svd_compact(t) +@tensor UdU[i; j] := conj(U[k; i]) * U[k; j] +U' * U ≈ id(domain(U)), UdU ≈ id(domain(U)) +``` + +The matrix-mul version satisfies orthogonality, but the `@tensor` version differs by the fermionic twist on the contracted index. +This is a common pitfall whenever an isometry obtained from a factorization is fed straight into an `@tensor` expression. + +- Computing a matrix function through a manual Taylor expansion: + matrix functions such as `exp`, `log`, `sqrt` are defined through the matrix product (trace formalism) and therefore have no immediate counterpart in terms of `@tensor` expressions. + In particular, replacing each matrix power by an `@tensor` self-contraction yields a different result, even at low order: + +```@example fermionnorm +function exp_via_tensor(t, order) + out = id(domain(t)) + tn = id(domain(t)) + for n in 1:order + @tensor next[a; b] := tn[a; c] * t[c; b] + tn = next + out += tn / factorial(n) + end + return out +end +exp(t) ≈ exp_via_tensor(t, 10) +``` + +The same Taylor expansion written with the matrix product instead does reproduce `exp(t)`, confirming that the discrepancy is in the contraction step rather than the truncation order: + +```@example fermionnorm +exp_via_mul(t, order) = sum(t^n / factorial(n) for n in 0:order) +exp(t) ≈ exp_via_mul(t, 10) +``` + +!!! note + Both the supertrace and the trace formalism constitute valid, consistent frameworks, each with their own advantages and disadvantages. + For practical applications, it can be convenient to select one or the other, and to take special care when trying to use properties of one framework in the other. + In general, each case must be carefully evaluated to check which framework is correct, but a good rule of thumb is to be careful when using properties of orthogonality in combination with `@tensor` expressions. + + +## Anyonic tensor contractions + +When `BraidingStyle(I) == Anyonic()`, the situation is more restrictive still. +The relevant group describing the exchange of two lines is no longer the permutation group but the full braid group, so even a double crossing is non-trivial and there is no preferred way to reorder lines in a diagram. +As a consequence, the implicit reordering that `@tensor` performs is no longer well-defined, and attempting an anyonic contraction with `@tensor` raises a `SectorMismatch` error. + +```@example anyoncontraction +using TensorKit # hide +V = Vect[FibonacciAnyon](:I => 1, :τ => 1) +A = randn(ComplexF64, V ← V ⊗ V) +B = randn(ComplexF64, V ⊗ V ← V) +try + @tensor C[i; j] := A[i; k l] * B[k l; j] +catch err + err +end +``` + +The way out is to write the contraction as a literal *planar* diagram, in which every required crossing is made explicit through a braiding tensor. +This is what the `@planar` macro provides. + +### The `@planar` macro + +The surface syntax of `@planar` is identical to that of `@tensor`, but with a number of additional restrictions. + +A diagram is *planar* in this context when it can be drawn on a sheet of paper without any of its lines crossing, and additionally with all open legs ending on the exterior of the diagram. +The second condition rules out arrangements in which an open leg is enclosed by contracted ones, even if the resulting diagram itself contains no crossings. + +For the macro to recognise this layout unambiguously, the codomain–domain separator `;` must be present in every index list. +It fixes which legs sit on the top (codomain) and which on the bottom (domain) of each tensor box, and changing the partition can change whether a given index pattern is planar. + +Planarity is moreover enforced for each binary contraction, not only for the overall expression. +The pairwise contraction order can therefore matter: an expression whose final layout is planar may still be rejected when an intermediate contraction produces a non-planar subdiagram. +Manually controlling the order, for instance via parentheses, NCON-style numbering, or the `order=...` keyword, is still supported, but must be done with care. + +Finally, the name `τ` is reserved for the braiding tensor: every literal crossing must be written out as a `τ[a b; c d]` factor, with its adjoint `τ'[a b; c d]` representing the inverse (under-)crossing. +The `BraidingTensor` itself does not need to be constructed by the user; the macro figures out the appropriate spaces from the surrounding contraction. +Any layout the macro cannot identify as planar is rejected at parse time with `ArgumentError("not a planar diagram expression: ...")`. + +To make this concrete, consider the contraction `A * B` for two anyonic tensors, written in a manifestly planar fashion: + +```@example anyoncontraction +@planar C1[i; j] := A[i; k l] * B[k l; j] +``` + +Inserting an explicit braiding tensor on the contracted legs gives a genuinely different result, reflecting the non-trivial R-symbols of the anyon braiding: + +```@example anyoncontraction +@planar C2[i; j] := A[i; k l] * τ[k l; m n] * B[m n; j] +C1 ≈ C2 +``` + +Both expressions correspond to valid, but distinct, tensor network diagrams, and the choice between them must be made explicit by the user. + +### The `@plansor` macro + +For code that should work uniformly across braiding styles, TensorKit provides the `@plansor` macro. +It inspects the `BraidingStyle` of the first non-braiding tensor in the expression and dispatches to `@tensor` for `Bosonic` sectors, and to `@planar` otherwise. +Any explicit `τ` factors that appear in the expression are silently removed in the bosonic case, where braidings are trivial, and faithfully evaluated otherwise. +This makes `@plansor` the natural choice for generic library code that wishes to remain correct regardless of the underlying symmetry. + +One important thing to note here is that in the specific case of `BraidingStyle(I) == Fermionic()`, the `@planar` macro is part of the trace formalism, and not the supertrace formalism. +Looking back at the examples in [Fermionic tensor contractions](@ref), these all consist of planar diagrams, so we could also have used `@planar` to achieve the desired outcomes. +Since in this case `@tensor` and `@planar` yield different results, the `@plansor` macro will fall back to `@tensor` only when `BraidingStyle(I) == Bosonic()`, and not when it is `Fermionic()`. + + +[^Mortier]: Mortier, Q., Devos, L., Burgelman, L., et al. (2025). Fermionic Tensor Network Methods. SciPost Physics 18, no. 1. [10.21468/SciPostPhys.18.1.012](https://doi.org/10.21468/SciPostPhys.18.1.012). diff --git a/docs/src/man/factorizations.md b/docs/src/man/factorizations.md new file mode 100644 index 000000000..ee10a251d --- /dev/null +++ b/docs/src/man/factorizations.md @@ -0,0 +1,46 @@ +# [Tensor factorizations](@id ss_tensor_factorization) + +```@setup tensors +using TensorKit +using LinearAlgebra +``` + +As tensors are linear maps, they support various kinds of factorizations. +These functions all interpret the provided `AbstractTensorMap` instances as a map from `domain` to `codomain`, which can be thought of as reshaping the tensor into a matrix according to the current bipartition of the indices. + +TensorKit's factorizations are provided by [MatrixAlgebraKit.jl](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl), which is used to supply both the interface, as well as the implementation of the various operations on the blocks of data. +For specific details on the provided functionality, we refer to its [documentation page](https://quantumkithub.github.io/MatrixAlgebraKit.jl/stable/user_interface/decompositions/). + +Finally, note that each of the factorizations takes the current partition of `domain` and `codomain` as the *axis* along which to matricize and perform the factorization. +In order to obtain factorizations according to a different bipartition of the indices, we can use any of the previously mentioned [index manipulations](@ref s_indexmanipulations) before the factorization. + +Some examples to conclude this section +```@repl tensors +V1 = SU₂Space(0 => 2, 1/2 => 1) +V2 = SU₂Space(0 => 1, 1/2 => 1, 1 => 1) + +t = randn(V1 ⊗ V1, V2); +U, S, Vh = svd_compact(t); +t ≈ U * S * Vh +D, V = eigh_full(t' * t); +D ≈ S * S +U' * U ≈ id(domain(U)) +S + +Q, R = left_orth(t; alg = :svd); +Q' * Q ≈ id(domain(Q)) +t ≈ Q * R + +U2, S2, Vh2, ε = svd_trunc(t; trunc = truncspace(V1)); +Vh2 * Vh2' ≈ id(codomain(Vh2)) +S2 +ε ≈ norm(block(S, Irrep[SU₂](1))) * sqrt(dim(Irrep[SU₂](1))) + +L, Q = right_orth(permute(t, ((1,), (2, 3)))); +codomain(L), domain(L), domain(Q) +Q * Q' +P = Q' * Q; +P ≈ P * P +t′ = permute(t, ((1,), (2, 3))); +t′ ≈ t′ * P +``` diff --git a/docs/src/man/fusiontrees.md b/docs/src/man/fusiontrees.md index 57fc7aa3b..56b6b2618 100644 --- a/docs/src/man/fusiontrees.md +++ b/docs/src/man/fusiontrees.md @@ -210,7 +210,7 @@ The main ingredient that we need is summarized in ``` We will only need the B-symbol and not the A-symbol. -Applying the left evaluation on the second sector of a splitting tensor thus yields a linear combination of fusion tensors (when `FusionStyle(I) == GenericFusion()`, or just a scalar times the corresponding fusion tensor otherwise), with corresponding ``Z`` ismorphism. +Applying the left evaluation on the second sector of a splitting tensor thus yields a linear combination of fusion tensors (when `FusionStyle(I) == GenericFusion()`, or just a scalar times the corresponding fusion tensor otherwise), with corresponding ``Z`` isomorphism. Taking the adjoint of this relation yields the required relation to transform a fusion tensor into a splitting tensor with an added ``Z^†`` isomorphism. However, we have to be careful if we bend a line on which a ``Z`` isomorphism (or its adjoint) is already present. diff --git a/docs/src/man/gradedspaces.md b/docs/src/man/gradedspaces.md index f2388023d..69342dcba 100644 --- a/docs/src/man/gradedspaces.md +++ b/docs/src/man/gradedspaces.md @@ -76,7 +76,7 @@ const SU₂Space = SU2Space To create specific instances of those types, one can e.g. just use `V = GradedSpace(a => n_a, b => n_b, c => n_c)` or `V = GradedSpace(iterator)` where `iterator` is any iterator (e.g. a dictionary or a generator) that yields `Pair{I, Int}` instances. With those constructions, `I` is inferred from the type of sectors. -However, it is often more convenient to specify the sector type explicitly (using one of the many alias provided), since then the sectors are automatically converted to the correct type. +However, it is often more convenient to specify the sector type explicitly (using one of the many aliases provided), since then the sectors are automatically converted to the correct type. Thereto, one can use `Vect[I]`, or when `I` corresponds to the irreducible representations of a group, `Rep[G]`. Some examples: @@ -101,8 +101,8 @@ The degeneracy dimensions `n_a` can be extracted as `dim(V, a)`, it properly ret With [`hassector(V, a)`](@ref) one can check if `V` contains a sector `a` with `dim(V, a) > 0`. Finally, `dim(V)` returns the total dimension of the space `V`, i.e. ``∑_a n_a d_a`` or thus `dim(V) = sum(dim(V, a) * dim(a) for a in sectors(V))`. Note that a representation space `V` has certain sectors `a` with dimensions `n_a`, then its dual `V'` will report to have sectors `dual(a)`, and `dim(V', dual(a)) == n_a`. -There is a subtelty regarding the difference between the dual of a representation space ``R_a^*``, on which the conjugate representation acts, and the representation space of the irrep `dual(a) == conj(a)` that is isomorphic to the conjugate representation, i.e. ``R_{\overline{a}} ≂ R_a^*`` but they are not equal. -We return to this in the section on [fusion trees](@ref sss_fusiontrees). +There is a subtlety regarding the difference between the dual of a representation space ``R_a^*``, on which the conjugate representation acts, and the representation space of the irrep `dual(a) == conj(a)` that is isomorphic to the conjugate representation, i.e. ``R_{\overline{a}} ≂ R_a^*`` but they are not equal. +We return to this in the section on [fusion trees](@ref s_fusiontrees). This is true also in more general fusion categories beyond the representation categories of groups. Other methods for `ElementarySpace`, such as [`dual`](@ref), [`fuse`](@ref) and [`flip`](@ref) also work. @@ -118,7 +118,7 @@ The function [`dims(W, as)`](@ref) returns the corresponding tuple with degenera [`hassector(W, as)`](@ref) is equivalent to `dim(W, as) > 0`. Finally, there is the function [`blocksectors(W)`](@ref) which returns a list (of type `Vector`) with all possible "block sectors" or total/coupled sectors that can result from fusing the individual uncoupled sectors in `W`. Correspondingly, [`blockdim(W, a)`](@ref) counts the total degeneracy dimension of the coupled sector `a` in `W`. -The machinery for computing this is the topic of the next section on [Fusion trees](@ref sss_fusiontrees), but first, it's time for some examples. +The machinery for computing this is the topic of the next section on [Fusion trees](@ref s_fusiontrees), but first, it's time for some examples. ## Examples diff --git a/docs/src/man/img/tensor-fermioniccontraction.svg b/docs/src/man/img/tensor-fermioniccontraction.svg new file mode 100644 index 000000000..c0ec84cf9 --- /dev/null +++ b/docs/src/man/img/tensor-fermioniccontraction.svg @@ -0,0 +1,317 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + B + A + X + + diff --git a/docs/src/man/indexmanipulations.md b/docs/src/man/indexmanipulations.md new file mode 100644 index 000000000..344dbe404 --- /dev/null +++ b/docs/src/man/indexmanipulations.md @@ -0,0 +1,126 @@ +# [Index manipulations](@id s_indexmanipulations) + +```@meta +CollapsedDocStrings = true +``` + +```@setup indexmanip +using TensorKit +using LinearAlgebra +``` + +A `TensorMap{T, S, N₁, N₂}` is a linear map from a domain (a `ProductSpace{S, N₂}`) to a codomain (a `ProductSpace{S, N₁}`). +In practice, the bipartition of the `N₁ + N₂` indices between domain and codomain rarely remains fixed: algorithms typically need to reshuffle indices between the two sides, reorder them, or change the arrow direction on individual indices before passing a tensor to a factorization or contraction. + +Index manipulations cover all such operations. +They act on the structure of the tensor data in a way that is fully determined by the categorical data of the `sectortype`, such that TensorKit automatically manipulates the tensor entries accordingly. +The operations fall into three groups, which mirror the structure of the source file: + +* **Reweighting**: [`flip`](@ref) and [`twist`](@ref) apply local isomorphisms to individual indices without changing the index structure. +* **Space insertion/removal**: [`insertleftunit`](@ref), [`insertrightunit`](@ref) and [`removeunit`](@ref) add or remove trivial (scalar) index factors. +* **Index rearrangements**: [`permute`](@ref), [`braid`](@ref), [`transpose`](@ref) and [`repartition`](@ref) reorder indices and/or move them between domain and codomain. + +Throughout this page, new index positions are specified using `Index2Tuple{N₁, N₂}`, i.e. a pair `(p₁, p₂)` of index tuples. +The indices listed in `p₁` form the new codomain and those in `p₂` form the new domain. +The following helpers retrieve the current index structure of a tensor: + +```@docs; canonical=false +numout +numin +numind +codomainind +domainind +allind +``` + +## Reweighting + +Reweighting operations modify the entries of a tensor by applying local isomorphisms to individual indices, without changing the number of indices or their partition between domain and codomain. +In particular, [`twist`](@ref) applies the topological spin (monoidal twist) to selected indices; this operation preserves the space of the indices and is completely trivial for `BraidingStyle(I) == Bosonic()`. +In contrast, [`flip`](@ref) changes the arrow direction on selected indices by applying a (non-canonical!) isomorphism between the index space and its dual. + +```@docs; canonical=false +twist(::AbstractTensorMap, ::Int) +twist! +flip(t::AbstractTensorMap, I) +``` + +## Inserting and removing unit spaces + +The next set of functions add or remove a trivial tensor product factor at a specified index position, without affecting any other indices. +We distinguish between [`insertleftunit`](@ref), which inserts a unit index before index `i` (the unit index becoming index `i`), +and [`insertrightunit`](@ref), which inserts after index `i` (the unit index becoming index `i + 1`); +[`removeunit`](@ref) undoes either insertion. + +For tensors `t` with `UnitStyle(sectortype(t)) = SimpleUnit()`, the only relevant difference between `insertleftunit(t, i + 1)` and `insertrightunit(t, i)` is that `insertleftunit(t, numout(t) + 1)` inserts the unit index as first index in the domain, whereas `insertrightunit(t, numout(t))` will insert the unit index as last index in the codomain. + +Passing `Val(i)` instead of an integer `i` for the position may improve type stability. + +```@docs; canonical=false +insertleftunit(::AbstractTensorMap, ::Val{i}) where {i} +insertrightunit(::AbstractTensorMap, ::Val{i}) where {i} +removeunit(::AbstractTensorMap, ::Val{i}) where {i} +``` + +## Index rearrangements + +These operations reorder indices and/or move them between domain and codomain by applying the transposing or braiding isomorphisms of the underlying category. +They form a hierarchy from most general to most restricted: + +- [`braid`](@ref) is the most general: it accepts any permutation and requires a `levels` argument — a tuple of heights, one per index — that determines whether each index crosses over or under the others it has to pass. +- [`permute`](@ref) is a simpler interface for sector types with a symmetric braiding (`BraidingStyle(I) isa SymmetricBraiding`), where over- and under-crossings are equivalent and `levels` is therefore not needed. +- [`transpose`](@ref) is restricted to *cyclic* permutations (indices do not cross). +- [`repartition`](@ref) only moves the codomain/domain boundary without reordering the indices at all. + +For plain tensors (`sectortype(t) == Trivial`), `permute` and `braid` act like `permutedims` on the underlying array: + +```@repl indexmanip +V = ℂ^2; +t = randn(V ⊗ V ← V ⊗ V); +ta = convert(Array, t); +t′ = permute(t, ((4, 2, 3), (1,))); +convert(Array, t′) ≈ permutedims(ta, (4, 2, 3, 1)) +``` + +```@docs; canonical=false +braid(::AbstractTensorMap, ::Index2Tuple, ::IndexTuple) +braid! +permute(::AbstractTensorMap, ::Index2Tuple) +permute!(::AbstractTensorMap, ::AbstractTensorMap, ::Index2Tuple) +transpose(::AbstractTensorMap, ::Index2Tuple) +transpose! +repartition(::AbstractTensorMap, ::Int, ::Int) +repartition! +``` + +## Fusing and splitting indices + +There is no dedicated functionality for fusing or splitting indices. +In the general case there is no canonical embedding of `V1 ⊗ V2` into the fused space `V = fuse(V1 ⊗ V2)`: any two such embeddings differ by a basis transform, i.e. there is a gauge freedom. +TensorKit resolves this by requiring the user to construct an explicit isomorphism — the *fuser* — and contract it with the tensor. +One particular isomorphism can be constructed using the [`unitary`](@ref) function. +It preserves norms and inner products, and has an inverse given by its adjoint. +For a plain tensor (`sectortype(t) == Trivial`), applying this particular `unitary` is equivalent to `reshape` on the underlying array. + +Fusing index `i` and `j = i+1` of a tensor `t` is then accomplished as + +```@repl indexmanip +t = randn(ℂ^2 ⊗ ℂ^3 ⊗ ℂ^4); +F = unitary(fuse(space(t, 2) ⊗ space(t, 3)), space(t, 2) ⊗ space(t, 3)); +@tensor t_fused[a, c] := F[c, b₁, b₂] * t[a, b₁, b₂] +``` + +The fusion is undone by contracting with the adjoint fuser, which is its inverse: + +```@repl indexmanip +@tensor t_split[a, b₁, b₂] := t_fused[a, c] * conj(F[c, b₁, b₂]); +t_split ≈ t +``` + +The resulting `unitary` is a dense `TensorMap`, and this fusion and splitting approach is not optimized for maximal performance. +However, because most tensor operations including tensor factorizations (SVD, QR, etc.) can be applied without needing any fusion, we do not expect fusion and splitting to be an essential part of performance critical parts of typical tensor algorithms. + +!!! warning + For `BraidingStyle(I) == Fermionic()`, special care has to be taken with these operations. + The definition of `unitary` is such that `F * F' ≈ I`, which may differ from the equivalent `@tensor` call depending on the duality of the spaces. + See also the section on [Fermionic tensor contractions](@ref). diff --git a/docs/src/man/intro.md b/docs/src/man/intro.md index 40f31fa8e..7e8e505fd 100644 --- a/docs/src/man/intro.md +++ b/docs/src/man/intro.md @@ -58,8 +58,8 @@ This simple example introduces two new concepts. The implementation of all of this is discussed in [Vector spaces](@ref s_spaces). 2. In the generic case, the identification between maps ``W → V`` and tensors in ``V ⊗ W^*`` is not an equivalence but an isomorphism, which needs to be defined. - Similarly, there is an isomorphism between between ``V ⊗ W`` and ``W ⊗ V`` that can be non-trivial (e.g. in the case of fermions / super vector spaces). - The correct formalism here is provided by theory of monoidal categories, the details of which are explained in the appendix. + Similarly, there is an isomorphism between ``V ⊗ W`` and ``W ⊗ V`` that can be non-trivial (e.g. in the case of fermions / super vector spaces). + The correct formalism here is provided by the theory of monoidal categories, the details of which are explained in the appendix. Nonetheless, we try to hide these canonical isomorphisms from the user wherever possible, and one does not need to know category theory to be able to use this package. This brings us to our final (yet formal) definition diff --git a/docs/src/man/linearalgebra.md b/docs/src/man/linearalgebra.md new file mode 100644 index 000000000..cd545d65a --- /dev/null +++ b/docs/src/man/linearalgebra.md @@ -0,0 +1,85 @@ +# [Basic linear algebra](@id ss_tensor_linalg) + +```@setup tensors +using TensorKit +using LinearAlgebra +``` + +`AbstractTensorMap` instances `t` represent linear maps, i.e. homomorphisms in a `𝕜`-linear category, just like matrices. +To a large extent, they follow the interface of `Matrix` in Julia's `LinearAlgebra` standard library. +Many methods from `LinearAlgebra` are (re)exported by TensorKit.jl, and can then be used without `using LinearAlgebra` explicitly. +In all of the following methods, the implementation acts directly on the underlying matrix blocks (typically using the same method) and never needs to perform any basis transforms. + +In particular, `AbstractTensorMap` instances can be composed, provided the domain of the first object coincides with the codomain of the second. +Composing tensor maps uses the regular multiplication symbol as in `t = t1 * t2`, which is also used for matrix multiplication. +TensorKit.jl also supports (and exports) the mutating method `mul!(t, t1, t2)`. +We can then also try to invert a tensor map using `inv(t)`, though this can only exist if the domain and codomain are isomorphic, which can e.g. be checked as `fuse(codomain(t)) == fuse(domain(t))`. +If the inverse is composed with another tensor `t2`, we can use the syntax `t1 \ t2` or `t2 / t1`. +However, this syntax also accepts instances `t1` whose domain and codomain are not isomorphic, and then amounts to `pinv(t1)`, the Moore-Penrose pseudoinverse. +This, however, is only really justified as minimizing the least squares problem if `InnerProductStyle(t) <: EuclideanProduct`. + +`AbstractTensorMap` instances behave themselves as vectors (i.e. they are `𝕜`-linear) and so they can be multiplied by scalars and, if they live in the same space, i.e. have the same domain and codomain, they can be added to each other. +There is also a `zero(t)`, the additive identity, which produces a zero tensor with the same domain and codomain as `t`. +In addition, `TensorMap` supports basic Julia methods such as `fill!` and `copy!`, as well as `copy(t)` to create a copy with independent data. +Aside from basic `+` and `*` operations, TensorKit.jl reexports a number of efficient in-place methods from `LinearAlgebra`, such as `axpy!` (for `y ← α * x + y`), `axpby!` (for `y ← α * x + β * y`), `lmul!` and `rmul!` (for `y ← α * y` and `y ← y * α`, which is typically the same) and `mul!`, which can also be used for out-of-place scalar multiplication `y ← α * x`. + +For `S = spacetype(t)` where `InnerProductStyle(S) <: EuclideanProduct`, we can compute `norm(t)`, and for two such instances, the inner product `dot(t1, t2)`, provided `t1` and `t2` have the same domain and codomain. +Furthermore, there is `normalize(t)` and `normalize!(t)` to return a scaled version of `t` with unit norm. +These operations should also exist for `InnerProductStyle(S) <: HasInnerProduct`, but require an interface for defining a custom inner product in these spaces. +Currently, there is no concrete subtype of `HasInnerProduct` that is not an `EuclideanProduct`. +In particular, `CartesianSpace`, `ComplexSpace` and `GradedSpace` all have `InnerProductStyle(S) <: EuclideanProduct`. + +With tensors that have `InnerProductStyle(t) <: EuclideanProduct` there is associated an adjoint operation, given by `adjoint(t)` or simply `t'`, such that `domain(t') == codomain(t)` and `codomain(t') == domain(t)`. +Note that for an instance `t::TensorMap{S, N₁, N₂}`, `t'` is simply stored in a wrapper called `AdjointTensorMap{S, N₂, N₁}`, which is another subtype of `AbstractTensorMap`. +This should be mostly invisible to the user, as all methods should work for this type as well. +It can be hard to reason about the index order of `t'`, i.e. index `i` of `t` appears in `t'` at index position `j = TensorKit.adjointtensorindex(t, i)`, where the latter method is typically not necessary and hence unexported. +There is also a plural `TensorKit.adjointtensorindices` to convert multiple indices at once. +Note that, because the adjoint interchanges domain and codomain, we have `space(t', j) == space(t, i)'`. + +`AbstractTensorMap` instances can furthermore be tested for exact (`t1 == t2`) or approximate (`t1 ≈ t2`) equality, though the latter requires that `norm` can be computed. + +When tensor map instances are endomorphisms, i.e. they have the same domain and codomain, there is a multiplicative identity which can be obtained as `one(t)` or `one!(t)`, where the latter overwrites the contents of `t`. +The multiplicative identity on a space `V` can also be obtained using `id(A, V)` as discussed [above](@ref ss_tensor_construction), such that for a general homomorphism `t′`, we have `t′ == id(codomain(t′)) * t′ == t′ * id(domain(t′))`. +Returning to the case of endomorphisms `t`, we can compute the trace via `tr(t)` and exponentiate them using `exp(t)`, or if the contents of `t` can be destroyed in the process, `exp!(t)`. +Furthermore, there are a number of tensor factorizations for both endomorphisms and general homomorphisms that we discuss on the [Tensor factorizations](@ref ss_tensor_factorization) page. + +Finally, there are a number of operations that also belong in this paragraph because of their analogy to common matrix operations. +The tensor product of two `TensorMap` instances `t1` and `t2` is obtained as `t1 ⊗ t2` and results in a new `TensorMap` with `codomain(t1 ⊗ t2) = codomain(t1) ⊗ codomain(t2)` and `domain(t1 ⊗ t2) = domain(t1) ⊗ domain(t2)`. +If we have two `TensorMap{T, S, N, 1}` instances `t1` and `t2` with the same codomain, we can combine them in a way that is analogous to `hcat`, i.e. we stack them such that the new tensor `catdomain(t1, t2)` has also the same codomain, but has a domain which is `domain(t1) ⊕ domain(t2)`. +Similarly, if `t1` and `t2` are of type `TensorMap{T, S, 1, N}` and have the same domain, the operation `catcodomain(t1, t2)` results in a new tensor with the same domain and a codomain given by `codomain(t1) ⊕ codomain(t2)`, which is the analogy of `vcat`. +Note that the direct sum only makes sense between `ElementarySpace` objects, i.e. there is no way to give a tensor product meaning to a direct sum of tensor product spaces. + +Time for some more examples: +```@repl tensors +using TensorKit # hide +V1 = ℂ^2 +t = randn(V1 ← V1 ⊗ V1 ⊗ V1) +t == t + zero(t) == t * id(domain(t)) == id(codomain(t)) * t +t2 = randn(ComplexF64, codomain(t), domain(t)); +dot(t2, t) +tr(t2' * t) +dot(t2, t) ≈ dot(t', t2') +dot(t2, t2) +norm(t2)^2 +t3 = copy!(similar(t, ComplexF64), t); +t3 == t +rmul!(t3, 0.8); +t3 ≈ 0.8 * t +axpby!(0.5, t2, 1.3im, t3); +t3 ≈ 0.5 * t2 + 0.8 * 1.3im * t +t4 = randn(fuse(codomain(t)), codomain(t)); +t5 = TensorMap{Float64}(undef, fuse(codomain(t)), domain(t)); +mul!(t5, t4, t) == t4 * t +inv(t4) * t4 ≈ id(codomain(t)) +t4 * inv(t4) ≈ id(fuse(codomain(t))) +t4 \ (t4 * t) ≈ t +t6 = randn(ComplexF64, V1, codomain(t)); +numout(t4) == numout(t6) == 1 +t7 = catcodomain(t4, t6); +foreach(println, (codomain(t4), codomain(t6), codomain(t7))) +norm(t7) ≈ sqrt(norm(t4)^2 + norm(t6)^2) +t8 = t4 ⊗ t6; +foreach(println, (codomain(t4), codomain(t6), codomain(t8))) +foreach(println, (domain(t4), domain(t6), domain(t8))) +norm(t8) ≈ norm(t4)*norm(t6) +``` diff --git a/docs/src/man/sectors.md b/docs/src/man/sectors.md index 9d5a26c2e..af0621285 100644 --- a/docs/src/man/sectors.md +++ b/docs/src/man/sectors.md @@ -33,10 +33,9 @@ The minimal data to completely specify a type of sector closely matches the [top * If the category is braided (see below), the R-symbol ``R^{ab}_c``; implemented as the function [`Rsymbol(a, b, c)`](@ref). Furthermore, sectors should provide information about the structure of their fusion rules. -For irreps of Abelian groups, we have that for every ``a`` and ``b``, there exists a unique ``c`` such that ``a ⊗ b = c``, i. -. there is only a single fusion channel. +For irreps of Abelian groups, we have that for every ``a`` and ``b``, there exists a unique ``c`` such that ``a ⊗ b = c``, i.e. there is only a single fusion channel. This follows simply from the fact that all irreps are one-dimensional. -In all other cases, there is at least one pair of (``a``, ``b``) exists such that ``a ⊗ b`` has multiple fusion outputs. +In all other cases, there is at least one pair of (``a``, ``b``) such that ``a ⊗ b`` has multiple fusion outputs. This is often referred to as non-abelian fusion, and is the case for the irreps of a non-abelian group or some more general fusion category. We however still distinguish between the case where all entries of ``N^{ab}_c ≦ 1``, i.e. they are zero or one. In that case, ``[F^{abc}_{d}]^f_e`` and ``R^{ab}_c`` are scalars. @@ -114,7 +113,7 @@ Hereto, the following methods are defined: * `findindex(::SectorValues{I}, c::I)`: reverse mapping that associates an index `i::Integer ∈ 1:length(values(I))` to a given sector `c::I`. The fallback implementation simply searches linearly through the `values(I)` iterator. -Note that `findindex` acts similar to `Base.indexin`, but with the order of the arguments reversed (so that is more similar to `getindex`), and returns an `Int` rather than an `Array{0, Union{Int, Nothing}}`. +Note that `findindex` acts similar to `Base.indexin`, but with the order of the arguments reversed (so that it is more similar to `getindex`), and returns an `Int` rather than an `Array{Union{Int, Nothing}}`. Secondly, it is often useful to know the scalar type in which the topological data in the F- and R-symbols are expressed. For this, the method [`sectorscalartype(I::Type{<:Sector})`](@ref) is provided, which has a default implementation that uses type inference on the return values of `Fsymbol` and `Rsymbol`. @@ -229,7 +228,7 @@ BraidingStyle(::Type{<:AbstractIrrep}) = Bosonic() We will need different data structures to represent irreps of different groups, but it would be convenient to easily obtain the relevant structure for a given group `G` in a uniform manner. Hereto, we define a singleton type `IrrepTable` with an associated exported constant `Irrep = IrrepTable()` as the only instance. -When a concrete type for representing the the irreps of a certain group `G` is iplemented, this type can the be "discovered" or obtained as `Irrep[G]`, provided it was registered by defining `Base.getindex(::IrrepTable, ::Type{G})` to return the concrete type. +When a concrete type for representing the irreps of a certain group `G` is implemented, this type can then be "discovered" or obtained as `Irrep[G]`, provided it was registered by defining `Base.getindex(::IrrepTable, ::Type{G})` to return the concrete type. Furthermore, we combine the more common functionality for irreps of abelian groups ```julia @@ -323,7 +322,7 @@ end Base.getindex(::IrrepTable, ::Type{ℤ{N}}) where {N} = N ≤ SMALL_ZN_CUTOFF ? ZNIrrep{N} : LargeZNIrrep{N} ... ``` -and continues along simular lines of the `U1Irrep` implementation above, by replacing the arithmetic with modulo `N` arithmetic. +and continues along similar lines of the `U1Irrep` implementation above, by replacing the arithmetic with modulo `N` arithmetic. The storage benefits for small `N` are not only due to a smaller integer type in the sector itself, but emerges as a result of the following distinction in the iterator size: ```julia @@ -460,9 +459,9 @@ FusionStyle(::Type{CU1Irrep}) = SimpleFusion() ... ``` The rest of the implementation can be read in the source code, but is rather long due to all the different cases for the arguments of `Fsymbol`. -For the dihedrial groups ``\mathsf{D}_N``, which can be intepreted as the semidirect product ``\mathbb{Z}_N ⋉ ℤ_2``, the representation theory is obtained quite similarly, and is implmented as the type [`DNIrrep{N}`](@ref). +For the dihedral groups ``\mathsf{D}_N``, which can be interpreted as the semidirect product ``\mathbb{Z}_N ⋉ ℤ_2``, the representation theory is obtained quite similarly, and is implemented as the type [`DNIrrep{N}`](@ref). -Of the aforementioned groups, only ``\mathsf{A}_4`` has a representation theory for which `FusionStyle(I) == GenericFusion()`, i.e. where fusion mulitplicities are required. +Of the aforementioned groups, only ``\mathsf{A}_4`` has a representation theory for which `FusionStyle(I) == GenericFusion()`, i.e. where fusion multiplicities are required. Another example where this does appear is for the irreps of `SU{N}` for ``N > 2``. Such sectors are supported through [SUNRepresentations.jl](https://github.com/QuantumKitHub/SUNRepresentations.jl), which implements numerical routines to compute the topological data of the representation theory of these groups, as no general analytic formula is available. @@ -541,7 +540,7 @@ TensorKit.BraidingStyle(::Type{I}) = ... # NoBraiding(), Bosonic(), Fermionic(), TensorKit.Rsymbol(a::I, b::I, c::I) = ... # only if BraidingStyle(I) != NoBraiding() Base.iterate(::TensorKit.SectorValues{I}[, state]) = ... -Base.IteratorSize(::Type{TensorKit.SectorValues{I}}) = ... # HasLenght() or IsInfinite() +Base.IteratorSize(::Type{TensorKit.SectorValues{I}}) = ... # HasLength() or IsInfinite() # if previous function returns HasLength(): Base.length(::TensorKit.SectorValues{I}) = ... # optional, but recommended if IteratorSize returns HasLength(): @@ -649,14 +648,14 @@ Note in particular how the `Rsymbol` values have opposite signs to the bosonic c ## Anyons Both `Bosonic` and `Fermionic` braiding styles are `SymmetricBraiding` styles, which means that exchanging two sectors twice is equivalent to the identity operation. -In tensor network diagrams, this implies that lines that cross twice are equivalent to them not crossing at all, or also, that there is no distinction betweeen a line crossing "above" or "below" another line. +In tensor network diagrams, this implies that lines that cross twice are equivalent to them not crossing at all, or also, that there is no distinction between a line crossing "above" or "below" another line. More technically, the relevant group describing the exchange processes is the permutation group, whereas in more general cases it would be the braid group. This more general case is denoted as the `Anyonic` braiding style in TensorKit.jl, because examples of this behaviour appear in the context of anyons in topological phases of matter. -There are currently two well-known sector types with `Anyonic` braiding style implemented in TensorKitSectors. -l, namely [`FibonacciAnyon`](@ref) and [`IsingAnyon`](@ref). Their values represent the (equivalence classes of) simple objects of the well-known Fibonicci and Ising fusion categories. -As an example, we illustrate below the Fibonacci anyons, which has only two distinct sectors, namely the unit sector `𝟙` and one non-trivial sector denoted as `τ`. +There are currently two well-known sector types with `Anyonic` braiding style implemented in TensorKitSectors.jl, namely [`FibonacciAnyon`](@ref) and [`IsingAnyon`](@ref). +Their values represent the (equivalence classes of) simple objects of the well-known Fibonacci and Ising fusion categories. +As an example, we illustrate below the Fibonacci anyons, which have only two distinct sectors, namely the unit sector `𝟙` and one non-trivial sector denoted as `τ`. The fusion rules are given by `τ ⊗ τ = 𝟙 ⊕ τ`, and the topological data is summarized by the following code ```@repl sectors diff --git a/docs/src/man/spaces.md b/docs/src/man/spaces.md index 52e2b2f92..63362cd9e 100644 --- a/docs/src/man/spaces.md +++ b/docs/src/man/spaces.md @@ -94,7 +94,7 @@ to denote for a vector space `V` whether it has an inner product and thus a cano This mapping is provided by the metric, but no further support for working with vector spaces with general metrics is currently implemented. A number of concrete elementary spaces are implemented in TensorKit.jl. -There is concrete type `GeneralSpace` which is completely characterized by its field `𝔽`, its dimension and whether its the dual and/or complex conjugate of $𝔽^d$. +There is a concrete type `GeneralSpace` which is completely characterized by its field `𝔽`, its dimension and whether it's the dual and/or complex conjugate of $𝔽^d$. ```@docs; canonical=false GeneralSpace @@ -138,7 +138,7 @@ InnerProductStyle(ℂ^5) This means that even for `ℂ^n`, arrows matter in the diagrammatic notation for categories or for tensors, and in particular that a contraction between two tensor indices will check that one is living in the space and the other in the dual space. This is in contrast with several other software packages, especially in the context of tensor networks, where arrows are only introduced when discussing symmetries. We believe that our more puristic approach can be useful to detect errors (e.g. unintended contractions). - Only with `ℝ^n` will their be no distinction between a space and its dual. + Only with `ℝ^n` will there be no distinction between a space and its dual. When creating tensors with indices in `ℝ^n` that have complex data, a one-time warning will be printed, but most operations should continue to work nonetheless. One more important concrete implementation of `ElementarySpace` with a `EuclideanInnerProduct()` is the [`GradedSpace`](@ref) type, which is used to represent a graded complex vector space, where the grading is provided by the irreducible representations of a group, or more generally, the simple objects of a unitary fusion category. @@ -253,7 +253,7 @@ Other scientific fields, like general relativity, might also benefit from tensor Vector spaces of the same `spacetype` can be given a partial order, based on whether there exist injective morphisms (a.k.a *monomorphisms*) or surjective morphisms (a.k.a. *epimorphisms*) between them. In particular, we define `ismonomorphic(V1, V2)`, with Unicode synonym `V1 ≾ V2` (obtained as `\precsim+TAB`), to express whether there exist monomorphisms in `V1 → V2`. Similarly, we define `isepimorphic(V1, V2)`, with Unicode synonym `V1 ≿ V2` (obtained as `\succsim+TAB`), to express whether there exist epimorphisms in `V1 → V2`. -Finally, we define `isisomorphic(V1, V2)`, with Unicode alternative `V1 ≅ V2` (obtained as `\cong+TAB`), to express whether there exist isomorphism in `V1 → V2`. +Finally, we define `isisomorphic(V1, V2)`, with Unicode alternative `V1 ≅ V2` (obtained as `\cong+TAB`), to express whether there exists an isomorphism in `V1 → V2`. In particular `V1 ≅ V2` if and only if `V1 ≾ V2 && V1 ≿ V2`. For completeness, we also export the strict comparison operators `≺` and `≻` (`\prec+TAB` and `\succ+TAB`), with definitions @@ -263,7 +263,7 @@ For completeness, we also export the strict comparison operators `≺` and `≻` ``` However, as we expect these to be less commonly used, no ASCII alternative is provided. -In the context of `InnerProductStyle(V) <: EuclideanInnerProduct`, `V1 ≾ V2` implies that there exists isometries ``W : V1 → V2`` such that ``W^† ∘ W = \mathrm{id}_{V1}``, while `V1 ≅ V2` implies that there exist unitaries ``U : V1 → V2`` such that ``U^† ∘ U = \mathrm{id}_{V1}`` and ``U ∘ U^† = \mathrm{id}_{V2}``. +In the context of `InnerProductStyle(V) <: EuclideanInnerProduct`, `V1 ≾ V2` implies that there exist isometries ``W : V1 → V2`` such that ``W^† ∘ W = \mathrm{id}_{V1}``, while `V1 ≅ V2` implies that there exist unitaries ``U : V1 → V2`` such that ``U^† ∘ U = \mathrm{id}_{V1}`` and ``U ∘ U^† = \mathrm{id}_{V2}``. Note that spaces that are isomorphic are not necessarily equal. One can be a dual space, and the other a normal space, or one can be an instance of `ProductSpace`, while the other is an `ElementarySpace`. diff --git a/docs/src/man/symmetries.md b/docs/src/man/symmetries.md index 21145b41b..61e03c9d3 100644 --- a/docs/src/man/symmetries.md +++ b/docs/src/man/symmetries.md @@ -53,7 +53,7 @@ We refer to the appendix on [categories](@ref s_categories), and in particular t Let the different irreps or sectors be labeled as ``a``, ``b``, ``c``, … First and foremost, we need to specify the *fusion rules* ``a ⊗ b = ⨁ N^{ab}_{c} c`` with ``N^{ab}_{c}`` some non-negative integers. The meaning of the fusion rules is that the space of covariant maps ``R_a ⊗ R_b → R_c`` (or vice versa ``R_c → R_a ⊗ R_b``) has dimension ``N^{ab}_c``. -In particular, there should always exists a unique trivial sector ``u`` (called the identity object ``I`` or ``1`` in the language of categories) such that ``a ⊗ u = a = u ⊗ a`` for every other sector ``a``. +In particular, there should always exist a unique trivial sector ``u`` (called the identity object ``I`` or ``1`` in the language of categories) such that ``a ⊗ u = a = u ⊗ a`` for every other sector ``a``. Furthermore, with respect to every sector ``a`` there should exist a unique sector ``\bar{a}`` such that ``N^{a\bar{a}}_{u} = 1``, whereas for all ``b \neq \bar{a}``, ``N^{ab}_{u} = 0``. For irreps of groups, ``\bar{a}`` corresponds to the complex conjugate of the representation ``a``, or some representation isomorphic to it. For example, for the representations of ``\mathsf{SU}_2``, the trivial sector corresponds to spin zero and all irreps are self-dual (i.e. ``a = \bar{a}``), meaning that the conjugate representation is isomorphic to the non-conjugated one (they are however not equal but related by a similarity transform). diff --git a/docs/src/man/tensormanipulations.md b/docs/src/man/tensormanipulations.md deleted file mode 100644 index 15285fd78..000000000 --- a/docs/src/man/tensormanipulations.md +++ /dev/null @@ -1,337 +0,0 @@ -# [Manipulating tensors](@id s_tensormanipulations) - -## [Vector space and linear algebra operations](@id ss_tensor_linalg) - -`AbstractTensorMap` instances `t` represent linear maps, i.e. homomorphisms in a `𝕜`-linear category, just like matrices. -To a large extent, they follow the interface of `Matrix` in Julia's `LinearAlgebra` standard library. -Many methods from `LinearAlgebra` are (re)exported by TensorKit.jl, and can then us be used without `using LinearAlgebra` explicitly. -In all of the following methods, the implementation acts directly on the underlying matrix blocks (typically using the same method) and never needs to perform any basis transforms. - -In particular, `AbstractTensorMap` instances can be composed, provided the domain of the first object coincides with the codomain of the second. -Composing tensor maps uses the regular multiplication symbol as in `t = t1 * t2`, which is also used for matrix multiplication. -TensorKit.jl also supports (and exports) the mutating method `mul!(t, t1, t2)`. -We can then also try to invert a tensor map using `inv(t)`, though this can only exist if the domain and codomain are isomorphic, which can e.g. be checked as `fuse(codomain(t)) == fuse(domain(t))`. -If the inverse is composed with another tensor `t2`, we can use the syntax `t1 \ t2` or `t2 / t1`. -However, this syntax also accepts instances `t1` whose domain and codomain are not isomorphic, and then amounts to `pinv(t1)`, the Moore-Penrose pseudoinverse. -This, however, is only really justified as minimizing the least squares problem if `InnerProductStyle(t) <: EuclideanProduct`. - -`AbstractTensorMap` instances behave themselves as vectors (i.e. they are `𝕜`-linear) and so they can be multiplied by scalars and, if they live in the same space, i.e. have the same domain and codomain, they can be added to each other. -There is also a `zero(t)`, the additive identity, which produces a zero tensor with the same domain and codomain as `t`. -In addition, `TensorMap` supports basic Julia methods such as `fill!` and `copy!`, as well as `copy(t)` to create a copy with independent data. -Aside from basic `+` and `*` operations, TensorKit.jl reexports a number of efficient in-place methods from `LinearAlgebra`, such as `axpy!` (for `y ← α * x + y`), `axpby!` (for `y ← α * x + β * y`), `lmul!` and `rmul!` (for `y ← α * y` and `y ← y * α`, which is typically the same) and `mul!`, which can also be used for out-of-place scalar multiplication `y ← α * x`. - -For `S = spacetype(t)` where `InnerProductStyle(S) <: EuclideanProduct`, we can compute `norm(t)`, and for two such instances, the inner product `dot(t1, t2)`, provided `t1` and `t2` have the same domain and codomain. -Furthermore, there is `normalize(t)` and `normalize!(t)` to return a scaled version of `t` with unit norm. -These operations should also exist for `InnerProductStyle(S) <: HasInnerProduct`, but require an interface for defining a custom inner product in these spaces. -Currently, there is no concrete subtype of `HasInnerProduct` that is not an `EuclideanProduct`. -In particular, `CartesianSpace`, `ComplexSpace` and `GradedSpace` all have `InnerProductStyle(S) <: EuclideanProduct`. - -With tensors that have `InnerProductStyle(t) <: EuclideanProduct` there is associated an adjoint operation, given by `adjoint(t)` or simply `t'`, such that `domain(t') == codomain(t)` and `codomain(t') == domain(t)`. -Note that for an instance `t::TensorMap{S, N₁, N₂}`, `t'` is simply stored in a wrapper called `AdjointTensorMap{S, N₂, N₁}`, which is another subtype of `AbstractTensorMap`. -This should be mostly invisible to the user, as all methods should work for this type as well. -It can be hard to reason about the index order of `t'`, i.e. index `i` of `t` appears in `t'` at index position `j = TensorKit.adjointtensorindex(t, i)`, where the latter method is typically not necessary and hence unexported. -There is also a plural `TensorKit.adjointtensorindices` to convert multiple indices at once. -Note that, because the adjoint interchanges domain and codomain, we have `space(t', j) == space(t, i)'`. - -`AbstractTensorMap` instances can furthermore be tested for exact (`t1 == t2`) or approximate (`t1 ≈ t2`) equality, though the latter requires that `norm` can be computed. - -When tensor map instances are endomorphisms, i.e. they have the same domain and codomain, there is a multiplicative identity which can be obtained as `one(t)` or `one!(t)`, where the latter overwrites the contents of `t`. -The multiplicative identity on a space `V` can also be obtained using `id(A, V)` as discussed [above](@ref ss_tensor_construction), such that for a general homomorphism `t′`, we have `t′ == id(codomain(t′)) * t′ == t′ * id(domain(t′))`. -Returning to the case of endomorphisms `t`, we can compute the trace via `tr(t)` and exponentiate them using `exp(t)`, or if the contents of `t` can be destroyed in the process, `exp!(t)`. -Furthermore, there are a number of tensor factorizations for both endomorphisms and general homomorphism that we discuss below. - -Finally, there are a number of operations that also belong in this paragraph because of their analogy to common matrix operations. -The tensor product of two `TensorMap` instances `t1` and `t2` is obtained as `t1 ⊗ t2` and results in a new `TensorMap` with `codomain(t1 ⊗ t2) = codomain(t1) ⊗ codomain(t2)` and `domain(t1 ⊗ t2) = domain(t1) ⊗ domain(t2)`. -If we have two `TensorMap{T, S, N, 1}` instances `t1` and `t2` with the same codomain, we can combine them in a way that is analogous to `hcat`, i.e. we stack them such that the new tensor `catdomain(t1, t2)` has also the same codomain, but has a domain which is `domain(t1) ⊕ domain(t2)`. -Similarly, if `t1` and `t2` are of type `TensorMap{T, S, 1, N}` and have the same domain, the operation `catcodomain(t1, t2)` results in a new tensor with the same domain and a codomain given by `codomain(t1) ⊕ codomain(t2)`, which is the analogy of `vcat`. -Note that direct sum only makes sense between `ElementarySpace` objects, i.e. there is no way to give a tensor product meaning to a direct sum of tensor product spaces. - -Time for some more examples: -```@repl tensors -using TensorKit # hide -V1 = ℂ^2 -t = randn(V1 ← V1 ⊗ V1 ⊗ V1) -t == t + zero(t) == t * id(domain(t)) == id(codomain(t)) * t -t2 = randn(ComplexF64, codomain(t), domain(t)); -dot(t2, t) -tr(t2' * t) -dot(t2, t) ≈ dot(t', t2') -dot(t2, t2) -norm(t2)^2 -t3 = copy!(similar(t, ComplexF64), t); -t3 == t -rmul!(t3, 0.8); -t3 ≈ 0.8 * t -axpby!(0.5, t2, 1.3im, t3); -t3 ≈ 0.5 * t2 + 0.8 * 1.3im * t -t4 = randn(fuse(codomain(t)), codomain(t)); -t5 = TensorMap{Float64}(undef, fuse(codomain(t)), domain(t)); -mul!(t5, t4, t) == t4 * t -inv(t4) * t4 ≈ id(codomain(t)) -t4 * inv(t4) ≈ id(fuse(codomain(t))) -t4 \ (t4 * t) ≈ t -t6 = randn(ComplexF64, V1, codomain(t)); -numout(t4) == numout(t6) == 1 -t7 = catcodomain(t4, t6); -foreach(println, (codomain(t4), codomain(t6), codomain(t7))) -norm(t7) ≈ sqrt(norm(t4)^2 + norm(t6)^2) -t8 = t4 ⊗ t6; -foreach(println, (codomain(t4), codomain(t6), codomain(t8))) -foreach(println, (domain(t4), domain(t6), domain(t8))) -norm(t8) ≈ norm(t4)*norm(t6) -``` - -## [Index manipulations](@id ss_indexmanipulation) - -In many cases, the bipartition of tensor indices (i.e. `ElementarySpace` instances) between the codomain and domain is not fixed throughout the different operations that need to be performed on that tensor map, i.e. we want to use the duality to move spaces from domain to codomain and vice versa. -Furthermore, we want to use the braiding to reshuffle the order of the indices. - -For this, we use an interface that is closely related to that for manipulating splitting- fusion tree pairs, namely [`braid`](@ref) and [`permute`](@ref), with the interface - -```julia -braid(t::AbstractTensorMap{T,S,N₁,N₂}, (p1, p2)::Index2Tuple{N₁′,N₂′}, levels::IndexTuple{N₁+N₂,Int}) -``` - -and - -```julia -permute(t::AbstractTensorMap{T,S,N₁,N₂}, (p1, p2)::Index2Tuple{N₁′,N₂′}; copy = false) -``` - -both of which return an instance of `AbstractTensorMap{T, S, N₁′, N₂′}`. - -In these methods, `p1` and `p2` specify which of the original tensor indices ranging from `1` to `N₁ + N₂` make up the new codomain (with `N₁′` spaces) and new domain (with `N₂′` spaces). -Hence, `(p1..., p2...)` should be a valid permutation of `1:(N₁ + N₂)`. -Note that, throughout TensorKit.jl, permutations are always specified using tuples of `Int`s, for reasons of type stability. -For `braid`, we also need to specify `levels` or depths for each of the indices of the original tensor, which determine whether indices will braid over or underneath each other (use the braiding or its inverse). -We refer to the section on [manipulating fusion trees](@ref ss_fusiontrees) for more details. - -When `BraidingStyle(sectortype(t)) isa SymmetricBraiding`, we can use the simpler interface of `permute`, which does not require the argument `levels`. -`permute` accepts a keyword argument `copy`. -When `copy == true`, the result will be a tensor with newly allocated data that can independently be modified from that of the input tensor `t`. -When `copy` takes the default value `false`, `permute` can try to return the result in a way that it shares its data with the input tensor `t`, though this is only possible in specific cases (e.g. when `sectortype(S) == Trivial` and `(p1..., p2...) = (1:(N₁+N₂)...)`). - -Both `braid` and `permute` come in a version where the result is stored in an already existing tensor, i.e. [`braid!(tdst, tsrc, (p1, p2), levels)`](@ref) and [`permute!(tdst, tsrc, (p1, p2))`](@ref). - -Another operation that belongs under index manipulations is taking the `transpose` of a tensor, i.e. `LinearAlgebra.transpose(t)` and `LinearAlgebra.transpose!(tdst, tsrc)`, both of which are reexported by TensorKit.jl. -Note that `transpose(t)` is not simply equal to reshuffling domain and codomain with `braid(t, (1:(N₁+N₂)...), reverse(domainind(tsrc)), reverse(codomainind(tsrc))))`. -Indeed, the graphical representation (where we draw the codomain and domain as a single object), makes clear that this introduces an additional (inverse) twist, which is then compensated in the `transpose` implementation. - -```@raw html -transpose -``` - -In categorical language, the reason for this extra twist is that we use the left coevaluation ``η``, but the right evaluation ``\tilde{ϵ}``, when repartitioning the indices between domain and codomain. - -There are a number of other index related manipulations. -We can apply a twist (or inverse twist) to one of the tensor map indices via [`twist(t, i; inv = false)`](@ref) or [`twist!(t, i; inv = false)`](@ref). -Note that the latter method does not store the result in a new destination tensor, but just modifies the tensor `t` in place. -Twisting several indices simultaneously can be obtained by using the defining property - -```math -θ_{V⊗W} = τ_{W,V} ∘ (θ_W ⊗ θ_V) ∘ τ_{V,W} = (θ_V ⊗ θ_W) ∘ τ_{W,V} ∘ τ_{V,W}, -``` - -but is currently not implemented explicitly. - -For all sector types `I` with `BraidingStyle(I) == Bosonic()`, all twists are `1` and thus have no effect. -Let us start with some examples, in which we illustrate that, albeit `permute` might act highly non-trivial on the fusion trees and on the corresponding data, after conversion to a regular `Array` (when possible), it just acts like `permutedims` - -```@repl tensors -domain(t) → codomain(t) -ta = convert(Array, t); -t′ = permute(t, (1, 2, 3, 4)); -domain(t′) → codomain(t′) -convert(Array, t′) ≈ ta -t′′ = permute(t, ((4, 2, 3), (1,))); -domain(t′′) → codomain(t′′) -convert(Array, t′′) ≈ permutedims(ta, (4, 2, 3, 1)) -transpose(t) -convert(Array, transpose(t)) ≈ permutedims(ta, (4, 3, 2, 1)) -dot(t2, t) ≈ dot(transpose(t2), transpose(t)) -transpose(transpose(t)) ≈ t -twist(t, 3) ≈ t -``` - -Note that `transpose` acts like one would expect on a `TensorMap{T, S, 1, 1}`. -On a `TensorMap{T, S, N₁, N₂}`, because `transpose` replaces the codomain with the dual of the domain, which has its tensor product operation reversed, this in the end amounts in a complete reversal of all tensor indices when representing it as a plain multi-dimensional `Array`. -Also, note that we have not defined the conjugation of `TensorMap` instances. -One definition that one could think of is `conj(t) = adjoint(transpose(t))`. -However note that `codomain(adjoint(tranpose(t))) == domain(transpose(t)) == dual(codomain(t))` and similarly `domain(adjoint(tranpose(t))) == dual(domain(t))`, where `dual` of a `ProductSpace` is composed of the dual of the `ElementarySpace` instances, in reverse order of tensor product. -This might be very confusing, and as such we leave tensor conjugation undefined. -However, note that we have a conjugation syntax within the context of [tensor contractions](@ref ss_tensor_contraction). - -To show the effect of `twist`, we now consider a type of sector `I` for which `BraidingStyle(I) != Bosonic()`. -In particular, we use `FibonacciAnyon`. -We cannot convert the resulting `TensorMap` to an `Array`, so we have to rely on indirect tests to verify our results. - -```@repl tensors -V1 = GradedSpace{FibonacciAnyon}(:I => 3, :τ => 2) -V2 = GradedSpace{FibonacciAnyon}(:I => 2, :τ => 1) -m = randn(Float32, V1, V2) -transpose(m) -twist(braid(m, ((2,), (1,)), (1, 2)), 1) -t1 = randn(V1 * V2', V2 * V1); -t2 = randn(ComplexF64, V1 * V2', V2 * V1); -dot(t1, t2) ≈ dot(transpose(t1), transpose(t2)) -transpose(transpose(t1)) ≈ t1 -``` - -A final operation that one might expect in this section is to fuse or join indices, and its inverse, to split a given index into two or more indices. -For a plain tensor (i.e. with `sectortype(t) == Trivial`) amount to the equivalent of `reshape` on the multidimensional data. -However, this represents only one possibility, as there is no canonically unique way to embed the tensor product of two spaces `V1 ⊗ V2` in a new space `V = fuse(V1 ⊗ V2)`. -Such a mapping can always be accompagnied by a basis transform. -However, one particular choice is created by the function `isomorphism`, or for `EuclideanProduct` spaces, `unitary`. -Hence, we can join or fuse two indices of a tensor by first constructing `u = unitary(fuse(space(t, i) ⊗ space(t, j)), space(t, i) ⊗ space(t, j))` and then contracting this map with indices `i` and `j` of `t`, as explained in the section on [contracting tensors](@ref ss_tensor_contraction). -Note, however, that a typical algorithm is not expected to often need to fuse and split indices, as e.g. tensor factorizations can easily be applied without needing to `reshape` or fuse indices first, as explained in the next section. - -## [Tensor factorizations](@id ss_tensor_factorization) - -As tensors are linear maps, they suport various kinds of factorizations. -These functions all interpret the provided `AbstractTensorMap` instances as a map from `domain` to `codomain`, which can be thought of as reshaping the tensor into a matrix according to the current bipartition of the indices. - -TensorKit's factorizations are provided by [MatrixAlgebraKit.jl](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl), which is used to supply both the interface, as well as the implementation of the various operations on the blocks of data. -For specific details on the provided functionality, we refer to its [documentation page](https://quantumkithub.github.io/MatrixAlgebraKit.jl/stable/user_interface/decompositions/). - -Finally, note that each of the factorizations takes the current partition of `domain` and `codomain` as the *axis* along which to matricize and perform the factorization. -In order to obtain factorizations according to a different bipartition of the indices, we can use any of the previously mentioned [index manipulations](@ref ss_indexmanipulation) before the factorization. - -Some examples to conclude this section -```@repl tensors -V1 = SU₂Space(0 => 2, 1/2 => 1) -V2 = SU₂Space(0 => 1, 1/2 => 1, 1 => 1) - -t = randn(V1 ⊗ V1, V2); -U, S, Vh = svd_compact(t); -t ≈ U * S * Vh -D, V = eigh_full(t' * t); -D ≈ S * S -U' * U ≈ id(domain(U)) -S - -Q, R = left_orth(t; alg = :svd); -Q' * Q ≈ id(domain(Q)) -t ≈ Q * R - -U2, S2, Vh2, ε = svd_trunc(t; trunc = truncspace(V1)); -Vh2 * Vh2' ≈ id(codomain(Vh2)) -S2 -ε ≈ norm(block(S, Irrep[SU₂](1))) * sqrt(dim(Irrep[SU₂](1))) - -L, Q = right_orth(permute(t, ((1,), (2, 3)))); -codomain(L), domain(L), domain(Q) -Q * Q' -P = Q' * Q; -P ≈ P * P -t′ = permute(t, ((1,), (2, 3))); -t′ ≈ t′ * P -``` - -## [Bosonic tensor contractions and tensor networks](@id ss_tensor_contraction) - -One of the most important operation with tensor maps is to compose them, more generally known as contracting them. -As mentioned in the section on [category theory](@ref s_categories), a typical composition of maps in a ribbon category can graphically be represented as a planar arrangement of the morphisms (i.e. tensor maps, boxes with lines eminating from top and bottom, corresponding to source and target, i.e. domain and codomain), where the lines connecting the source and targets of the different morphisms should be thought of as ribbons, that can braid over or underneath each other, and that can twist. -Technically, we can embed this diagram in ``ℝ × [0,1]`` and attach all the unconnected line endings corresponding objects in the source at some position ``(x,0)`` for ``x∈ℝ``, and all line endings corresponding to objects in the target at some position ``(x,1)``. -The resulting morphism is then invariant under what is known as *framed three-dimensional isotopy*, i.e. three-dimensional rearrangements of the morphism that respect the rules of boxes connected by ribbons whose open endings are kept fixed. -Such a two-dimensional diagram cannot easily be encoded in a single line of code. - -However, things simplify when the braiding is symmetric (such that over- and under- crossings become equivalent, i.e. just crossings), and when twists, i.e. self-crossings in this case, are trivial. -This amounts to `BraidingStyle(I) == Bosonic()` in the language of TensorKit.jl, and is true for any subcategory of ``\mathbf{Vect}``, i.e. ordinary tensors, possibly with some symmetry constraint. -The case of ``\mathbf{SVect}`` and its subcategories, and more general categories, are discussed below. - -In the case of trivial twists, we can deform the diagram such that we first combine every morphism with a number of coevaluations ``η`` so as to represent it as a tensor, i.e. with a trivial domain. -We can then rearrange the morphism to be all ligned up horizontally, where the original morphism compositions are now being performed by evaluations ``ϵ``. -This process will generate a number of crossings and twists, where the latter can be omitted because they act trivially. -Similarly, double crossings can also be omitted. -As a consequence, the diagram, or the morphism it represents, is completely specified by the tensors it is composed of, and which indices between the different tensors are connect, via the evaluation ``ϵ``, and which indices make up the source and target of the resulting morphism. -If we also compose the resulting morphisms with coevaluations so that it has a trivial domain, we just have one type of unconnected lines, henceforth called open indices. -We sketch such a rearrangement in the following picture - -```@raw html -tensor unitary -``` - -Hence, we can now specify such a tensor diagram, henceforth called a tensor contraction or also tensor network, using a one-dimensional syntax that mimicks [abstract index notation](https://en.wikipedia.org/wiki/Abstract_index_notation) and specifies which indices are connected by the evaluation map using Einstein's summation conventation. -Indeed, for `BraidingStyle(I) == Bosonic()`, such a tensor contraction can take the same format as if all tensors were just multi-dimensional arrays. -For this, we rely on the interface provided by the package [TensorOperations.jl](https://github.com/QuantumKitHub/TensorOperations.jl). - -The above picture would be encoded as -```julia -@tensor E[a, b, c, d, e] := A[v, w, d, x] * B[y, z, c, x] * C[v, e, y, b] * D[a, w, z] -``` -or -```julia -@tensor E[:] := A[1, 2, -4, 3] * B[4, 5, -3, 3] * C[1, -5, 4, -2] * D[-1, 2, 5] -``` -where the latter syntax is known as NCON-style, and labels the unconnected or outgoing indices with negative integers, and the contracted indices with positive integers. - -A number of remarks are in order. -TensorOperations.jl accepts both integers and any valid variable name as dummy label for indices, and everything in between `[ ]` is not resolved in the current context but interpreted as a dummy label. -Here, we label the indices of a `TensorMap`, like `A::TensorMap{T, S, N₁, N₂}`, in a linear fashion, where the first position corresponds to the first space in `codomain(A)`, and so forth, up to position `N₁`. -Index `N₁ + 1` then corresponds to the first space in `domain(A)`. -However, because we have applied the coevaluation ``η``, it actually corresponds to the corresponding dual space, in accordance with the interface of [`space(A, i)`](@ref) that we introduced [above](@ref ss_tensor_properties), and as indiated by the dotted box around ``A`` in the above picture. -The same holds for the other tensor maps. -Note that our convention also requires that we braid indices that we brought from the domain to the codomain, and so this is only unambiguous for a symmetric braiding, where there is a unique way to permute the indices. - -With the current syntax, we create a new object `E` because we use the definition operator `:=`. -Furthermore, with the current syntax, it will be a `Tensor`, i.e. it will have a trivial domain, and correspond to the dotted box in the picture above, rather than the actual morphism `E`. -We can also directly define `E` with the correct codomain and domain by rather using -```julia -@tensor E[a b c;d e] := A[v, w, d, x] * B[y, z, c, x] * C[v, e, y, b] * D[a, w, z] -``` -or -```julia -@tensor E[(a, b, c);(d, e)] := A[v, w, d, x] * B[y, z, c, x] * C[v, e, y, b] * D[a, w, z] -``` -where the latter syntax can also be used when the codomain is empty. -When using the assignment operator `=`, the `TensorMap` `E` is assumed to exist and the contents will be written to the currently allocated memory. -Note that for existing tensors, both on the left hand side and right hand side, trying to specify the indices in the domain and the codomain seperately using the above syntax, has no effect, as the bipartition of indices are already fixed by the existing object. -Hence, if `E` has been created by the previous line of code, all of the following lines are now equivalent -```julia -@tensor E[(a, b, c);(d, e)] = A[v, w, d, x] * B[y, z, c, x] * C[v, e, y, b] * D[a, w, z] -@tensor E[a, b, c, d, e] = A[v w d; x] * B[(y, z, c); (x, )] * C[v e y; b] * D[a, w, z] -@tensor E[a b; c d e] = A[v; w d x] * B[y, z, c, x] * C[v, e, y, b] * D[a w; z] -``` -and none of those will or can change the partition of the indices of `E` into its codomain and its domain. - -Two final remarks are in order. -Firstly, the order of the tensors appearing on the right hand side is irrelevant, as we can reorder them by using the allowed moves of the Penrose graphical calculus, which yields some crossings and a twist. -As the latter is trivial, it can be omitted, and we just use the same rules to evaluate the newly ordered tensor network. -For the particular case of matrix-matrix multiplication, which also captures more general settings by appropriotely combining spaces into a single line, we indeed find - -```@raw html -tensor contraction reorder -``` - -or thus, the following two lines of code yield the same result -```julia -@tensor C[i, j] := B[i, k] * A[k, j] -@tensor C[i, j] := A[k, j] * B[i, k] -``` -Reordering of tensors can be used internally by the `@tensor` macro to evaluate the contraction in a more efficient manner. -In particular, the NCON-style of specifying the contraction gives the user control over the order, and there are other macros, such as `@tensoropt`, that try to automate this process. -There is also an `@ncon` macro and `ncon` function, an we recommend reading the [manual of TensorOperations.jl](https://quantumkithub.github.io/TensorOperations.jl/stable/) to learn more about the possibilities and how they work. - -A final remark involves the use of adjoints of tensors. -The current framework is such that the user should not be too worried about the actual bipartition into codomain and domain of a given `TensorMap` instance. -Indeed, for tensor contractions the `@tensor` macro figures out the correct manipulations automatically. -However, when wanting to use the `adjoint` of an instance `t::TensorMap{T, S, N₁, N₂}`, the resulting `adjoint(t)` is an `AbstractTensorMap{T, S, N₂, N₁}` and one needs to know the values of `N₁` and `N₂` to know exactly where the `i`th index of `t` will end up in `adjoint(t)`, and hence the index order of `t'`. -Within the `@tensor` macro, one can instead use `conj()` on the whole index expression so as to be able to use the original index ordering of `t`. -For example, for `TensorMap{T, S, 1, 1}` instances, this yields exactly the equivalence one expects, namely one between the following two expressions: - -```julia -@tensor C[i, j] := B'[i, k] * A[k, j] -@tensor C[i, j] := conj(B[k, i]) * A[k, j] -``` - -For e.g. an instance `A::TensorMap{T, S, 3, 2}`, the following two syntaxes have the same effect within an `@tensor` expression: `conj(A[a, b, c, d, e])` and `A'[d, e, a, b, c]`. - -Some examples: - -## Fermionic tensor contractions - -TODO - -## Anyonic tensor contractions - -TODO diff --git a/docs/src/man/tensors.md b/docs/src/man/tensors.md index 2921e5f1b..3482ddc00 100644 --- a/docs/src/man/tensors.md +++ b/docs/src/man/tensors.md @@ -5,7 +5,7 @@ using TensorKit using LinearAlgebra ``` -This last page explains how to create and manipulate tensors in TensorKit.jl. +This page explains how to construct and access tensors in TensorKit.jl. As this is probably the most important part of the manual, we will also focus more strongly on the usage and interface, and less so on the underlying implementation. The only aspect of the implementation that we will address is the storage of the tensor data, as this is important to know how to create and initialize a tensor, but will in fact also shed light on how some of the methods work. @@ -13,7 +13,7 @@ As mentioned, all tensors in TensorKit.jl are interpreted as linear maps (morphi The overall type for all such tensor maps is `AbstractTensorMap{T, S, N₁, N₂}`. Note that we place information about the codomain before that of the domain. Indeed, we have already encountered the constructor for the concrete parametric type `TensorMap` in the form `TensorMap(..., codomain, domain)`. -This convention is opposite to the mathematical notation, e.g. ``\mathrm{Hom}(W, V)`` or ``f : W → V``, but originates from the fact that a normal matrix is also denoted as having size `m × n` or is constructed in Julia as `Array(..., (m, n))`, where the first integer `m` refers to the codomain being `m`- dimensional, and the seond integer `n` to the domain being `n`-dimensional. +This convention is opposite to the mathematical notation, e.g. ``\mathrm{Hom}(W, V)`` or ``f : W → V``, but originates from the fact that a normal matrix is also denoted as having size `m × n` or is constructed in Julia as `Array(..., (m, n))`, where the first integer `m` refers to the codomain being `m`- dimensional, and the second integer `n` to the domain being `n`-dimensional. This also explains why we have consistently used the symbol ``W`` for spaces in the domain and ``V`` for spaces in the codomain. A tensor map ``t : (W_1 ⊗ … ⊗ W_{N_2}) → (V_1 ⊗ … ⊗ V_{N_1})`` will be created in Julia as `TensorMap(..., V1 ⊗ ... ⊗ VN₁, W1 ⊗ ... ⊗ WN₂)`. @@ -22,12 +22,12 @@ Furthermore, the abstract type `AbstractTensor{T, S, N}` is just a synonym for ` Currently, `AbstractTensorMap` has three subtypes. `TensorMap` provides the actual implementation, where the data of the tensor is stored in a `DenseArray` (more specifically a `DenseMatrix` as will be explained below). `AdjointTensorMap` is a simple wrapper type to denote the adjoint of an existing `TensorMap` object. -`DiagonalTensorMap` provides an efficient representations of diagonal tensor maps. +`DiagonalTensorMap` provides an efficient representation of diagonal tensor maps. In the future, additional types could be defined, to deal with sparse data, static data, etc... ## [Storage of tensor data](@id ss_tensor_storage) -Before discussion how to construct and initalize a `TensorMap`, let us discuss what is meant by 'tensor data' and how it can efficiently and compactly be stored. +Before discussing how to construct and initialize a `TensorMap`, let us discuss what is meant by 'tensor data' and how it can efficiently and compactly be stored. Let us first discuss the case `sectortype(S) == Trivial` sector, i.e. the case of no symmetries. In that case the data of a tensor `t = TensorMap(..., V1 ⊗ ... ⊗ VN₁, W₁ ⊗ ... ⊗ WN₂)` can just be represented as a multidimensional array of size @@ -35,7 +35,7 @@ In that case the data of a tensor `t = TensorMap(..., V1 ⊗ ... ⊗ VN₁, W₁ (dim(V1), dim(V2), …, dim(VN₁), dim(W1), …, dim(WN₂)) ``` -which can also be reshaped into matrix of size +which can also be reshaped into a matrix of size ```julia (dim(V1) * dim(V2) * … * dim(VN₁), dim(W1) * dim(W2) * … * dim(WN₂)) @@ -46,12 +46,12 @@ In particular, given a second tensor `t2` whose domain matches with the codomain Similarly, tensor factorizations such as the singular value decomposition, which we discuss below, can act directly on this matrix representation. !!! note - One might wonder if it would not have been more natural to represent the tensor data as `(dim(V1), dim(V2), …, dim(VN₁), dim(WN₂), …, dim(W1))` given how employing the duality naturally reverses the tensor product, as encountered with the interface of [`repartition`](@ref) for [fusion trees](@ref ss_fusiontrees). + One might wonder if it would not have been more natural to represent the tensor data as `(dim(V1), dim(V2), …, dim(VN₁), dim(WN₂), …, dim(W1))` given how employing the duality naturally reverses the tensor product, as encountered with the interface of [`repartition`](@ref) for [fusion trees](@ref s_fusiontrees). However, such a representation, when plainly `reshape`d to a matrix, would not have the above properties and would thus not constitute the matrix representation of the tensor in a compatible basis. Now consider the case where `sectortype(S) == I` for some `I` which has `FusionStyle(I) == UniqueFusion()`, i.e. the representations of an Abelian group, e.g. `I == Irrep[ℤ₂]` or `I == Irrep[U₁]`. In this case, the tensor data is associated with sectors `(a1, a2, …, aN₁) ∈ sectors(V1 ⊗ V2 ⊗ … ⊗ VN₁)` and `(b1, …, bN₂) ∈ sectors(W1 ⊗ … ⊗ WN₂)` such that they fuse to a same common charge, i.e. `(c = first(⊗(a1, …, aN₁))) == first(⊗(b1, …, bN₂))`. -The data associated with this takes the form of a multidimensional array with size `(dim(V1, a1), …, dim(VN₁, aN₁), dim(W1, b1), …, dim(WN₂, bN₂))`, or equivalently, a matrix of with row size `dim(V1, a1) * … * dim(VN₁, aN₁) == dim(codomain, (a1, …, aN₁))` and column size `dim(W1, b1) * … * dim(WN₂, aN₂) == dim(domain, (b1, …, bN₂))`. +The data associated with this takes the form of a multidimensional array with size `(dim(V1, a1), …, dim(VN₁, aN₁), dim(W1, b1), …, dim(WN₂, bN₂))`, or equivalently, a matrix with row size `dim(V1, a1) * … * dim(VN₁, aN₁) == dim(codomain, (a1, …, aN₁))` and column size `dim(W1, b1) * … * dim(WN₂, aN₂) == dim(domain, (b1, …, bN₂))`. However, there are multiple combinations of `(a1, …, aN₁)` giving rise to the same `c`, and so there is data associated with all of these, as well as all possible combinations of `(b1, …, bN₂)`. Stacking all matrices for different `(a1, …)` and a fixed value of `(b1, …)` underneath each other, and for fixed value of `(a1, …)` and different values of `(b1, …)` next to each other, gives rise to a larger block matrix of all data associated with the central sector `c`. @@ -63,13 +63,13 @@ For example, composing two tensor maps amounts to multiplying the matrices corre We directly concatenate these blocks as consecutive entries in a single larger `DenseVector`, together with metadata to retrieve a block by using the corresponding coupled sector `c` as key. For a given tensor `t`, we can access a specific block as `block(t, c)`, whereas `blocks(t)` yields an iterator over pairs `c => block(t, c)`. -The subblocks corresponding to a particular combination of sectors then correspond to a particular view for some range of the rows and some range of the colums, i.e. `view(block(t, c), m₁:m₂, n₁:n₂)` where the ranges `m₁:m₂` associated with `(a1, …, aN₁)` and `n₁:n₂` associated with `(b₁, …, bN₂)` are stored within the fields of the instance `t` of type `TensorMap`. +The subblocks corresponding to a particular combination of sectors then correspond to a particular view for some range of the rows and some range of the columns, i.e. `view(block(t, c), m₁:m₂, n₁:n₂)` where the ranges `m₁:m₂` associated with `(a1, …, aN₁)` and `n₁:n₂` associated with `(b₁, …, bN₂)` are stored within the fields of the instance `t` of type `TensorMap`. This `view` can then lazily be reshaped to a multidimensional array, for which we rely on the package [Strided.jl](https://github.com/Jutho/Strided.jl). Indeed, the data in this `view` is not contiguous, because the stride between the different columns is larger than the length of the columns. Nonetheless, this does not pose a problem and even as multidimensional array there is still a definite stride associated with each dimension. When `FusionStyle(I) isa MultipleFusion`, things become slightly more complicated. -Not only do `(a1, …, aN₁)` give rise to different coupled sectors `c`, there can be multiply ways in which they fuse to `c`. +Not only do `(a1, …, aN₁)` give rise to different coupled sectors `c`, there can be multiple ways in which they fuse to `c`. These different possibilities are enumerated by the iterator `fusiontrees((a1, …, aN₁), c)` and `fusiontrees((b1, …, bN₂), c)`, and with each of those, there is tensor data that takes the form of a multidimensional array, or, after reshaping, a matrix of size `(dim(codomain, (a1, …, aN₁)), dim(domain, (b1, …, bN₂))))`. Again, we can stack all such matrices with the same value of `f₁ ∈ fusiontrees((a1, …, aN₁), c)` horizontally (as they all have the same number of rows), and with the same value of `f₂ ∈ fusiontrees((b1, …, bN₂), c)` vertically (as they have the same number of columns). What emerges is a large matrix of size `(blockdim(codomain, c), blockdim(domain, c))` containing all the tensor data associated with the coupled sector `c`, where `blockdim(P, c) = sum(dim(P, s) * length(fusiontrees(s, c)) for s in sectors(P))` for some instance `P` of `ProductSpace`. @@ -81,7 +81,7 @@ Schur's lemma now tells that there is a unitary basis transform which makes this The reason for this extra identity is that the group representation is recoupled to act as ``⨁_{c} 𝟙 ⊗ u_c(g)`` for all ``g ∈ \mathsf{I}``, with ``u_c(g)`` the matrix representation of group element ``g`` according to the irrep ``c``. In the abelian case, `dim(c) == 1`, i.e. all irreducible representations are one-dimensional and Schur's lemma only dictates that all off-diagonal blocks are zero. However, in this case the basis transform to the block diagonal representation is not simply a permutation matrix, but a more general unitary matrix composed of the different fusion trees. -Indeed, let us denote the fusion trees `f₁ ∈ fusiontrees((a1, …, aN₁), c)` as ``X^{a_1, …, a_{N₁}}_{c,α}`` where ``α = (e_1, …, e_{N_1-2}; μ₁, …, μ_{N_1-1})`` is a collective label for the internal sectors `e` and the vertex degeneracy labels `μ` of a generic fusion tree, as discussed in the [corresponding section](@ref ss_fusiontrees). +Indeed, let us denote the fusion trees `f₁ ∈ fusiontrees((a1, …, aN₁), c)` as ``X^{a_1, …, a_{N₁}}_{c,α}`` where ``α = (e_1, …, e_{N_1-2}; μ₁, …, μ_{N_1-1})`` is a collective label for the internal sectors `e` and the vertex degeneracy labels `μ` of a generic fusion tree, as discussed in the [corresponding section](@ref s_fusiontrees). The tensor is then represented as ```@raw html @@ -119,7 +119,7 @@ The thick vertical line represents the separation between the two different coup Dashed vertical lines represent different ways of reaching the coupled sector, corresponding to different `α`. In this example, the first sector ``(a…)`` has one fusion tree to ``c``, labeled by ``c,α``, and two fusion trees to ``c'``, labeled by ``c',α`` and ``c',α'``. The second sector has only a fusion tree to ``c``, labeled by ``c,α'``. -The third sector only has a fusion tree to ``c'``, labeld by ``c', α''``. +The third sector only has a fusion tree to ``c'``, labeled by ``c', α''``. Finally then, because the fusion trees do not act on the spaces ``ℂ^{n_{a_1} × … n_{a_{N_1}}}``, the dotted lines which represent the different ``n_{(a…)} = n_{a_1} × … n_{a_{N_1}}`` dimensions are also drawn vertically. In particular, for a given sector ``(a…)`` and a specific fusion tree ``X^{(a…)}_{c,α} : R_{(a…)}→R_c``, the action is ``X^{(a…)}_{c,α} ⊗ 𝟙_{n_{(a…)}}``, which corresponds to the diagonal green blocks in this drawing where the same matrix ``X^{(a…)}_{c,α}`` (the fusion tree) is repeated along the diagonal. Note that the fusion tree is not a vector or single column, but a matrix with number of rows equal to ``\mathrm{dim}(R_{(a\ldots)}) = d_{a_1} d_{a_2} … d_{a_{N_1}} `` and number of columns equal to ``d_c``. @@ -155,7 +155,7 @@ Tensor{eltype::Type{<:Number}}(undef, codomain) ``` Here, `f` is any of the typical functions from Base that normally create arrays, namely `zeros`, `ones`, `rand`, `randn` and `Random.randexp`. Remember that `one(codomain)` is the empty `ProductSpace{S, 0}()`. -The third and fourth calling syntax use the `UndefInitializer` from Julia Base and generates a `TensorMap` with unitialized data, which can thus contain `NaN`s. +The third and fourth calling syntax use the `UndefInitializer` from Julia Base and generates a `TensorMap` with uninitialized data, which can thus contain `NaN`s. In all of these constructors, the last two arguments can be replaced by `domain → codomain` or `codomain ← domain`, where the arrows are obtained as `\rightarrow+TAB` and `\leftarrow+TAB` and create a `HomSpace` as explained in the section on [Spaces of morphisms](@ref ss_homspaces). Some examples are perhaps in order @@ -259,7 +259,7 @@ for (c, b) in blocks(t4) println() end ``` -While this indeed does not require considering the internal structure of the representation spaces, it still requires knowing the precise row and column indices corresponding to each set of uncoupled sectors in the codmain and domain respectively to correctly assign the nonzero entries in each block. +While this indeed does not require considering the internal structure of the representation spaces, it still requires knowing the precise row and column indices corresponding to each set of uncoupled sectors in the codomain and domain respectively to correctly assign the nonzero entries in each block. Perhaps the most natural way of constructing a particular `TensorMap` is to directly assign the data slices for each splitting - fusion tree pair using the `fusiontrees(::TensorMap)` method. This returns an iterator over all tuples `(f₁, f₂)` of splitting - fusion tree pairs corresponding to all ways in which the set of domain uncoupled sectors can fuse to a coupled sector and split back into the set of codomain uncoupled sectors. diff --git a/docs/src/man/tutorial.md b/docs/src/man/tutorial.md index 3039a584c..4f4f19890 100644 --- a/docs/src/man/tutorial.md +++ b/docs/src/man/tutorial.md @@ -57,7 +57,7 @@ scalarAA = dot(A, A) normA² = norm(A)^2 ``` -More generally, our tensor objects implement the full interface layed out in [VectorInterface.jl](https://github.com/Jutho/VectorInterface.jl). +More generally, our tensor objects implement the full interface laid out in [VectorInterface.jl](https://github.com/Jutho/VectorInterface.jl). If two tensors live on different spaces, these operations have no meaning and are thus not allowed @@ -100,7 +100,7 @@ U ``` Note that the `svd_compact` routine returns the decomposition of the linear map as three factors, `U`, `S` and `Vd`, each of them a `TensorMap`, such that `Vd` is already what is commonly called `V'`. -Furthermore, observe that `U` is printed differently then `A`, i.e. as a `TensorMap((ℝ^3 ⊗ ℝ^4) ← ProductSpace(ℝ^2))`. +Furthermore, observe that `U` is printed differently than `A`, i.e. as a `TensorMap((ℝ^3 ⊗ ℝ^4) ← ProductSpace(ℝ^2))`. What this means is that tensors (or more appropriately, `TensorMap` instances) in TensorKit.jl are always considered to be linear maps between two `ProductSpace` instances, with ```@repl tutorial @@ -363,5 +363,5 @@ Instead, recoupling relations are used to symbolically manipulate the basis of f In fact, this formalism extends beyond the case of group representations on vector spaces, and can also deal with super vector spaces (to describe fermions) and more general (unitary) fusion categories. Support for all of these generalizations is present in TensorKit.jl. Indeed, all of these concepts will be explained throughout the remainder of this manual, including several details regarding their implementation. -However, to just use tensors and their manipulations (contractions, factorizations, ...) in higher level algorithms (e.g. tensoer network algorithms), one does not need to know or understand most of these details, and one can immediately refer to the general interface of the `TensorMap` type, discussed on the [last page](@ref s_tensors). +However, to just use tensors and their manipulations (contractions, factorizations, ...) in higher level algorithms (e.g. tensor network algorithms), one does not need to know or understand most of these details, and one can immediately refer to the general interface of the `TensorMap` type, discussed on the [last page](@ref s_tensors). Adhering to this interface should yield code and algorithms that are oblivious to the underlying symmetries and can thus work with arbitrary symmetric tensors. diff --git a/src/tensors/indexmanipulations.jl b/src/tensors/indexmanipulations.jl index 5fa236caa..58f8985ed 100644 --- a/src/tensors/indexmanipulations.jl +++ b/src/tensors/indexmanipulations.jl @@ -6,7 +6,7 @@ # flip # ------ """ - flip(t::AbstractTensorMap, I) -> t′::AbstractTensorMap + flip(t::AbstractTensorMap, I; inv::Bool = false) -> t′::AbstractTensorMap Return a new tensor that is isomorphic to `t` but where the arrows on the indices `i` that satisfy `i ∈ I` are flipped, i.e. `space(t′, i) = flip(space(t, i))`.