Tutorial

Before discussing at length all aspects of this package, both its usage and implementation, we start with a short tutorial to sketch the main capabilities. Thereto, we start by loading TensorLabXD.jl

julia> using TensorLabXD

Cartesian tensors

The most important objects in TensorLabXD.jl are tensors, which we now create with random normally distributed entries in the following manner

julia> A = Tensor(randn, ℝ^3 ⊗ ℝ^2 ⊗ ℝ^4)TensorMap((ℝ^3 ⊗ ℝ^2 ⊗ ℝ^4) ← ProductSpace{CartesianSpace, 0}()):
[:, :, 1] =
  0.49551113670302604  -1.662317488281836
 -0.6690449521339833   -0.42809769631463307
  0.34475265000491656  -0.3649841738811947

[:, :, 2] =
 -0.4464613269477426   -1.8645875877480311
  0.24394971977128915  -1.2834697401922586
  0.22990069701659682   0.9423239060770728

[:, :, 3] =
  1.0469619498724498    0.45271122620040294
 -0.28162889191365664  -1.0264470429888852
 -1.356762293429346    -0.7097909606117588

[:, :, 4] =
 -0.9720873909668194  -0.34271368610106023
 -0.8039115975185488  -1.3346943257664887
 -0.9508234614697689   0.15451823844243748

The tensor is created by specifying the vector space associated to each of the tensor indices, in this case ℝ^n (\bbR+TAB). The tensor lives in the tensor product of the index spaces, which can be obtained by typing ⊗ (\otimes+TAB) or *. The tensor A is printed as an instance of a parametric type TensorMap.

Briefly sidetracking into the nature of ℝ^n:

julia> V = ℝ^3ℝ^3
julia> typeof(V)CartesianSpace
julia> V == CartesianSpace(3)true
julia> supertype(CartesianSpace)EuclideanSpace{ℝ}
julia> supertype(EuclideanSpace)InnerProductSpace{𝕜} where 𝕜
julia> supertype(InnerProductSpace)ElementarySpace{𝕜} where 𝕜
julia> supertype(ElementarySpace)VectorSpace

We have seen that ℝ^n can also be created using the longer syntax CartesianSpace(n). It is subtype of EuclideanSpace{ℝ}, a space with a standard Euclidean inner product over the real numbers. Furthermore,

julia> W = ℝ^3 ⊗ ℝ^2 ⊗ ℝ^4(ℝ^3 ⊗ ℝ^2 ⊗ ℝ^4)
julia> typeof(W)ProductSpace{CartesianSpace,3}
julia> supertype(ProductSpace)CompositeSpace{S} where S<:ElementarySpace
julia> supertype(CompositeSpace)VectorSpace

The tensor product of a number of CartesianSpaces is some generic parametric ProductSpace type, specifically ProductSpace{CartesianSpace,N} for the tensor product of N instances of CartesianSpace.

Tensors behave like vectors (but not VectorSpace instance), so we can compute linear combinations provided they live in the same space.

julia> B = Tensor(randn, ℝ^3 * ℝ^2 * ℝ^4);
julia> C = 0.5*A + 2.5*BTensorMap((ℝ^3 ⊗ ℝ^2 ⊗ ℝ^4) ← ProductSpace{CartesianSpace, 0}()):
[:, :, 1] =
 -2.005599669053141   -0.5220323830149305
 -1.2088560851352586  -0.5331522026139472
 -3.863864028252894   -1.3782257163333647

[:, :, 2] =
 1.6665785284159282  -6.895537035586623
 1.124339434027722   -1.433443904441217
 2.418716557532789    0.7125710512968026

[:, :, 3] =
  2.8149733964937034   0.8962786705039005
  0.8129560226722392  -2.4668245628167607
 -2.494564629175893   -1.2529887661349441

[:, :, 4] =
  1.8422789055051239  -0.18057011007941234
  0.7562930741367853   1.6935849737294448
 -5.079598834463508    1.993939015007102

Tensors also have inner product and norm, which they inherit from the Euclidean inner product on the individual ℝ^n spaces:

julia> scalarBA = dot(B,A)6.520499951721173
julia> scalarAA = dot(A,A)19.412059423278077
julia> normA² = norm(A)^219.41205942327808

If two tensors live in different spaces, these operations have no meaning and are thus not allowed

julia> B′ = Tensor(randn, ℝ^4 * ℝ^2 * ℝ^3);
julia> space(B′) == space(A)false
julia> C′ = 0.5*A + 2.5*B′ERROR: SpaceMismatch()
julia> scalarBA′ = dot(B′,A)ERROR: SpaceMismatch()

However, in this particular case, we can reorder the indices of B′ to match space of A, using the routine permute (we deliberately choose not to overload permutedims from Julia Base):

julia> space(permute(B′,(3,2,1))) == space(A)true

We can contract two tensors using Einstein summation convention, which takes the interface from TensorContractionsXS.jl. TensorLabXD.jl reexports the @tensor macro

julia> @tensor D[a,b,c,d] := A[a,b,e]*B[d,c,e]TensorMap((ℝ^3 ⊗ ℝ^2 ⊗ ℝ^2 ⊗ ℝ^3) ← ProductSpace{CartesianSpace, 0}()):
[:, :, 1, 1] =
 -0.7298068013214084    0.18460535160575847
 -0.21940065739121142  -2.7682205242594717
 -2.266086402876952     0.5346156315702271

[:, :, 2, 1] =
  1.4103492385383414   4.364624950605062
 -0.7371255941340689   2.7383849299996323
 -0.8658191195604291  -2.4836247368803757

[:, :, 1, 2] =
 -0.4032467174884391   -0.15229665467308542
 -0.14809759094202574  -1.3748429885574813
 -0.9865258101789941    0.3062651796830693

[:, :, 2, 2] =
 -1.6580120465150159   0.12524841885363527
 -0.5309731750515478   0.002754914134600318
  0.04548655077459111  0.44875070363693353

[:, :, 1, 3] =
 -0.1817389340174916  1.267860132069076
  2.990111402626114   2.712195818580729
  2.3920125142203004  1.6886950246523993

[:, :, 2, 3] =
 -1.4014916052608384    0.18964446818155992
 -0.1716083797269435   -0.5737148176485254
 -0.38426292972925613   0.639011965024643
julia> @tensor d = A[a,b,c]*A[a,b,c]19.412059423278077
julia> d ≈ scalarAA ≈ normA²true

The := is to create a new tensor D. The = write the contraction result in an existing tensor d, which would yield an error if no tensor d exists. If the contraction yields a scalar, regular assignment with = can be used.

We can also factorize a tensor. With a plain Julia Array, one would apply permutedims and reshape to cast the array into a matrix before applying e.g. the singular value decomposition. With TensorLabXD.jl, one just specifies which indices go to the left (rows) and right (columns)

julia> U, S, Vd = tsvd(A, (1,3), (2,));
julia> @tensor A′[a,b,c] := U[a,c,d]*S[d,e]*Vd[e,b];
julia> A ≈ A′true
julia> UTensorMap((ℝ^3 ⊗ ℝ^4) ← ProductSpace(ℝ^2)):
[:, :, 1] =
 -0.35956495300097246  -0.5054247252546143   …  -0.18357121773543753
 -0.17375989904646988  -0.29174240806757007     -0.4111312265561198
 -0.05490873534868604   0.2558660291335364      -0.05875199589220592

[:, :, 2] =
  0.4713725853203608   0.12905358805138395  …  -0.33119750995406727
 -0.1962903485515649   0.30866235417931775     -0.10075869314809352
  0.1975889923700912  -0.06351889034985764     -0.4046100702516337

The tsvd routine returns the decomposition of the linear map as three factors, U, S and Vd, each of them a TensorMap, such that Vd is what is commonly called V'.

Notice that U is printed as TensorMap((ℝ^3 ⊗ ℝ^4) ← ProductSpace(ℝ^2)), which is a linear map between two ProductSpace instances, with

julia> codomain(U)(ℝ^3 ⊗ ℝ^4)
julia> domain(U)ProductSpace(ℝ^2)
julia> codomain(A)(ℝ^3 ⊗ ℝ^2 ⊗ ℝ^4)
julia> domain(A)ProductSpace{CartesianSpace, 0}()

Hence, a Tensor instance such as A is just a specific case of TensorMap with an empty domain, i.e. a ProductSpace{CartesianSpace,0} instance. For example, we can represent a vector v and matrix m as

julia> v = Tensor(randn, ℝ^3)TensorMap(ProductSpace(ℝ^3) ← ProductSpace{CartesianSpace, 0}()):
 -0.02725393758662904
  1.304537315609476
 -0.42187201496482446
julia> m1 = TensorMap(randn, ℝ^4, ℝ^3)TensorMap(ProductSpace(ℝ^4) ← ProductSpace(ℝ^3)):
 -1.8395305450131687    1.0782623278366903    -1.1348632331915538
  1.1438372476071392    0.028282818104834323  -0.7641240904139287
  2.208173927277968    -0.8311354625119804     0.3338476040791284
  0.41017418420382823  -1.3213397479144595     0.5176749041505545
julia> m2 = TensorMap(randn, ℝ^4 → ℝ^2) # alternative syntax for TensorMap(randn, ℝ^2, ℝ^4)TensorMap(ProductSpace(ℝ^2) ← ProductSpace(ℝ^4)):
  0.34352172056703006  -0.603966351819401     …  -1.049581426497123
 -0.692647301427436     0.012863553654253802     -0.007205891105239502
julia> w = m1 * v # matrix vector productTensorMap(ProductSpace(ℝ^4) ← ProductSpace{CartesianSpace, 0}()):
  1.9355349322374036
  0.3280844923588929
 -1.2852696209917154
 -1.953308424279085
julia> m3 = m2 * m1 # matrix matrix productTensorMap(ProductSpace(ℝ^2) ← ProductSpace(ℝ^3)):
 -1.6128264391840086   1.687317043371392   -0.4504537628850122
  5.00222530465704    -2.1357578938203026   1.334360457724273

Note that for the construction of m1, in accordance with how one specifies the dimensions of a matrix (e.g. randn(4,3)), the first space is the codomain and the second the domain. This is somewhat opposite to the general notation for a function f:domain→codomain, so that we also support this more mathematical notation, as illustrated in the construction of m2. There is a third syntax which mixes both like TensorMap(randn, codomain←domain).

This 'matrix vector' or 'matrix matrix' product can be computed between any two TensorMap instances for which the domain of the first matches with the codomain of the second, e.g.

julia> v′ = v ⊗ vTensorMap((ℝ^3 ⊗ ℝ^3) ← ProductSpace{CartesianSpace, 0}()):
  0.0007427771139758711  -0.03555377857904925   0.011497673565396757
 -0.03555377857904925     1.7018176078175775   -0.5503477859329727
  0.011497673565396757   -0.5503477859329727    0.17797599701048109
julia> m1′ = m1 ⊗ m1TensorMap((ℝ^4 ⊗ ℝ^4) ← (ℝ^3 ⊗ ℝ^3)):
[:, :, 1, 1] =
  3.3838726260364456  -2.104123555497124   …  -0.7545279406188
 -2.104123555497124    1.308363649013476       0.4691725098992106
 -4.06200338792951     2.5257915872154784      0.905735939201404
 -0.7545279406188      0.4691725098992106      0.168242861387276

[:, :, 2, 1] =
 -1.983496487592595    1.2333566132711866    …   0.4422753706781352
 -0.0520271078028942   0.032350940815607056      0.011600881843135681
  1.5288990703344354  -0.9506836998283904       -0.34091031029872304
  2.4306448266286487  -1.5113976204083865       -0.5419794531569054

[:, :, 3, 1] =
  2.0876155818682656  -1.298098837064366   …  -0.46549160085726443
  1.4056296044968262  -0.8740335964093771     -0.3134239754160255
 -0.6141228650830196   0.3818673245701082      0.13693566865155912
 -0.9522787985717094   0.5921358375188599      0.21233688149274865

[:, :, 1, 2] =
 -1.983496487592595   -0.0520271078028942    …   2.4306448266286487
  1.2333566132711866   0.032350940815607056     -1.5113976204083865
  2.3809907590950283   0.062453381529040425     -2.9177479804207525
  0.4422753706781352   0.011600881843135681     -0.5419794531569054

[:, :, 2, 2] =
  1.1626496476317982     0.030496297287500345   …  -1.4247508724493907
  0.030496297287500345   0.0007999177999511442     -0.037371211744952296
 -0.8961820585557922    -0.02350685310670369        1.098212322518348
 -1.4247508724493907    -0.037371211744952296       1.7459387294186475

[:, :, 3, 2] =
 -1.2236802715973976   -0.032097130398220895  …   1.4995398984627162
 -0.8239262204858163   -0.02161158265869912       1.0096675330029061
  0.35997529471706274   0.009442151060904734     -0.44112610901576177
  0.5581893472120124    0.014641305151527676     -0.6840244273519357

[:, :, 1, 3] =
  2.0876155818682656    1.4056296044968262  …  -0.9522787985717094
 -1.298098837064366    -0.8740335964093771      0.5921358375188599
 -2.505975402559966    -1.6873188936570302      1.1431162261513756
 -0.46549160085726443  -0.3134239754160255      0.21233688149274865

[:, :, 2, 3] =
 -1.2236802715973976    -0.8239262204858163   …   0.5581893472120124
 -0.032097130398220895  -0.02161158265869912      0.014641305151527676
  0.9432250782065036     0.635090629302727       -0.4302579708920163
  1.4995398984627162     1.0096675330029061      -0.6840244273519357

[:, :, 3, 3] =
  1.2879145580499871   0.8671763358067063   …  -0.587490215466426
  0.8671763358067063   0.5838856255509138      -0.3955678652641602
 -0.3788713713584934  -0.25510099680383336      0.17282452644255505
 -0.587490215466426   -0.3955678652641602       0.2679873063872858
julia> w′ = m1′ * v′TensorMap((ℝ^4 ⊗ ℝ^4) ← ProductSpace{CartesianSpace, 0}()):
  3.7462954739112497   0.6350189956860127   …  -3.7806966886257682
  0.6350189956860127   0.10763943412639247     -0.6408502027999526
 -2.487684248772993   -0.42167703114737365      2.5105279781531045
 -3.7806966886257682  -0.6408502027999526       3.815413800359643
julia> w′ ≈ w ⊗ wtrue

Another example involves checking that U from the singular value decomposition is a left isometric tensor

julia> codomain(U)(ℝ^3 ⊗ ℝ^4)
julia> domain(U)ProductSpace(ℝ^2)
julia> space(U)(ℝ^3 ⊗ ℝ^4) ← ℝ^2
julia> U'*U # should be the identity on the corresponding domain = codomainTensorMap(ProductSpace(ℝ^2) ← ProductSpace(ℝ^2)):
  1.0                     -1.2785673698038158e-16
 -1.2785673698038158e-16   0.9999999999999999
julia> U'*U ≈ one(U'*U)true
julia> P = U*U' # should be a projectorTensorMap((ℝ^3 ⊗ ℝ^4) ← (ℝ^3 ⊗ ℝ^4)):
[:, :, 1, 1] =
  0.35147906961819236   0.2425653410263545   …  -0.09011165024134074
 -0.03004791913608823   0.2503953171257802       0.10033349447116559
  0.11288129100932806  -0.12194152030083837     -0.16959693621947103

[:, :, 2, 1] =
 -0.03004791913608823   0.06249057545539533   …  0.09690819092984142
  0.06872240345113424  -0.009894309713209861     0.09121607941888217
 -0.02924387587132068  -0.031991110265282216     0.08962979259217466

[:, :, 3, 1] =
  0.11288129100932806   0.053251800902505766  …  -0.05536131885503624
 -0.02924387587132068   0.07700749021944217       0.0026658870608916524
  0.04205637912342     -0.02659991361910722      -0.0767204982901599

[:, :, 1, 2] =
 0.2425653410263545     0.2721089814876387   …   0.05003920527533216
 0.06249057545539533    0.18728781074589923      0.1947926163475896
 0.053251800902505766  -0.13751835818549835     -0.022521669945717446

[:, :, 2, 2] =
  0.2503953171257802     0.18728781074589923  …  -0.04867269400671766
 -0.009894309713209861   0.18038608155258318      0.08884399863613407
  0.07700749021944217   -0.09425286171234995     -0.10774744804816026

[:, :, 3, 2] =
 -0.12194152030083837   -0.13751835818549835  …  -0.02593234022625593
 -0.031991110265282216  -0.09425286171234995     -0.09879443400984599
 -0.02659991361910722    0.06950207429584093      0.010667742794153145

[:, :, 1, 3] =
  0.08314744710621719  -0.06609309428212204   …  -0.15366018039965887
 -0.10525309188126912   0.042155588978652576     -0.12431740929090027
  0.05577265794979442   0.03408074468525326      -0.15158413765707354

[:, :, 2, 3] =
 0.1281050374644951     0.14979578056244544  …   0.03295219758111409
 0.037841941888406234   0.09975362788490413      0.11017479594879213
 0.02669040385560122   -0.07572967263743255     -0.006398544332168311

[:, :, 3, 3] =
 -0.08727277395639435   0.10361711128315597  …  0.19798813844920712
  0.13769479373883303  -0.03951039903434825     0.1715465074999822
 -0.06660181858549412  -0.05322411669771836     0.19001836428043942

[:, :, 1, 4] =
 -0.09011165024134074   0.05003920527533216  …  0.1433901825806459
  0.09690819092984142  -0.04867269400671766     0.10884288818484547
 -0.05536131885503624  -0.02593234022625593     0.14479102320000098

[:, :, 2, 4] =
 0.10033349447116559     0.1947926163475896   …  0.10884288818484547
 0.09121607941888217     0.08884399863613407     0.17918119969445115
 0.0026658870608916524  -0.09879443400984599     0.06492276204689565

[:, :, 3, 4] =
 -0.16959693621947103  -0.022521669945717446  …  0.14479102320000098
  0.08962979259217466  -0.10774744804816026      0.06492276204689565
 -0.0767204982901599    0.010667742794153145     0.16716110597034975
julia> P*P ≈ Ptrue

Here, the adjoint of a TensorMap results in a new tensor map (actually a simple wrapper of type AdjointTensorMap <: AbstractTensorMap) with domain and codomain interchanged.

Our original tensor A living in ℝ^4 * ℝ^2 * ℝ^3 is isomorphic to a linear map ℝ^3 → ℝ^4 * ℝ^2. This is where the full power of permute emerges. It allows to specify a permutation where some indices go to the codomain, and others go to the domain, as

julia> A2 = permute(A,(1,2),(3,))TensorMap((ℝ^3 ⊗ ℝ^2) ← ProductSpace(ℝ^4)):
[:, :, 1] =
  0.49551113670302604  -1.662317488281836
 -0.6690449521339833   -0.42809769631463307
  0.34475265000491656  -0.3649841738811947

[:, :, 2] =
 -0.4464613269477426   -1.8645875877480311
  0.24394971977128915  -1.2834697401922586
  0.22990069701659682   0.9423239060770728

[:, :, 3] =
  1.0469619498724498    0.45271122620040294
 -0.28162889191365664  -1.0264470429888852
 -1.356762293429346    -0.7097909606117588

[:, :, 4] =
 -0.9720873909668194  -0.34271368610106023
 -0.8039115975185488  -1.3346943257664887
 -0.9508234614697689   0.15451823844243748
julia> codomain(A2)(ℝ^3 ⊗ ℝ^2)
julia> domain(A2)ProductSpace(ℝ^4)

In fact, tsvd(A, (1,3),(2,)) is a shorthand for tsvd(permute(A,(1,3),(2,))), where tsvd(A::TensorMap) will just compute the singular value decomposition according to the given codomain and domain of A.

The @tensor macro treats all indices at the same footing and thus does not distinguish between codomain and domain. The linear numbering is first all indices in the codomain, followed by all indices in the domain. However, when @tensor creates a new tensor (using :=), the default syntax creates a Tensor, i.e. with all indices in the codomain.

julia> @tensor A′[a,b,c] := U[a,c,d]*S[d,e]*Vd[e,b];
julia> codomain(A′)(ℝ^3 ⊗ ℝ^2 ⊗ ℝ^4)
julia> domain(A′)ProductSpace{CartesianSpace, 0}()
julia> @tensor A2′[(a,b);(c,)] := U[a,c,d]*S[d,e]*Vd[e,b];
julia> codomain(A2′)(ℝ^3 ⊗ ℝ^2)
julia> domain(A2′)ProductSpace(ℝ^4)
julia> @tensor A2′′[a b; c] := U[a,c,d]*S[d,e]*Vd[e,b];
julia> A2 ≈ A2′ == A2′′true

As illustrated for A2′ and A2′′, additional syntax is available that enables one to immediately specify the desired codomain and domain indices.

Complex tensors

To create a complex tensor, we work with complex vector spaces

julia> A = Tensor(randn, ComplexF64, ℂ^3 ⊗ ℂ^2 ⊗ ℂ^4)TensorMap((ℂ^3 ⊗ ℂ^2 ⊗ ℂ^4) ← ProductSpace{ComplexSpace, 0}()):
[:, :, 1] =
  0.3026296594997782 - 0.3280663810672145im  …  0.055624942853064976 - 0.20643959687043806im
 0.37673110224326967 + 1.1716266056947116im      0.44424650986199526 - 0.21687236195161047im
  -1.925837880815753 - 0.4506180892947755im       0.8998321551257689 - 0.37399521776642997im

[:, :, 2] =
    -1.21900301328605 - 0.04718860517331444im   …  -0.19932867756972422 - 0.570318936133236im
 -0.20838748635286428 - 0.3367165830795446im        0.08109709898712866 + 1.6832769568782144im
  -0.5870099373021512 + 0.003215350188337075im      0.06222451386036586 + 0.4212655977639447im

[:, :, 3] =
 -0.26529852880882826 + 0.8126221491097887im  …   0.20699436210486735 - 1.5170059036944157im
  -0.3560716360433608 + 0.5679982354876247im      -0.5878733085402784 - 0.18225552972137954im
    0.553177169265601 - 1.4738214754847463im     0.011879176478060504 - 0.38202465776345823im

[:, :, 4] =
 -1.4901987352734223 + 0.9700448409150365im   …  1.1444790105281866 - 0.614504475436288im
 0.20919108149481502 - 0.10095389671411428im      0.620791681816405 + 0.2178943113783455im
 -0.5626618449763482 + 0.30528150203663507im     1.0345286217534222 + 1.01258135907421im

where ℂ is obtained as \bbC+TAB and we also have the non-Unicode alternative ℂ^n == ComplexSpace(n). Most functionality works exactly the same with real tensors

julia> B = Tensor(randn, ℂ^3 * ℂ^2 * ℂ^4);
julia> C = im*A + (2.5-0.8im)*BTensorMap((ℂ^3 ⊗ ℂ^2 ⊗ ℂ^4) ← ProductSpace{ComplexSpace, 0}()):
[:, :, 1] =
  -2.456822131851595 + 1.1937939836337972im   …   -1.0196603058655642 + 0.4479769117285857im
 -1.1761179255743668 + 0.37816832460475935im     -0.19298347630007778 + 0.5754003781025355im
 -0.2868087961695454 - 1.6898612774671702im        1.1328454091333109 + 0.6570000938883669im

[:, :, 2] =
 -0.9283288912331696 - 0.9068374144359752im   …   2.2756275604215066 - 0.7450274373419707im
  0.4300740618089774 - 0.23826187954628275im     -6.6252498806599505 + 1.6625284345972844im
  1.5958803892967772 - 1.098720573937388im        2.2655068961876297 - 0.7975426842041379im

[:, :, 3] =
 -1.245617872243075 - 0.1267398974061766im  …  -4.703484775368731 + 2.1975513794050743im
 1.9173647557619393 - 1.1513877932432213im     3.6148202334556947 - 1.6862940137352593im
   1.43347197324761 + 0.5660890099814846im     0.5149847470164078 - 0.030668052082883382im

[:, :, 4] =
 -0.7100919408786712 - 1.573383663285059im    …   0.4955177334092525 + 1.182554767976838im
  1.8926335130478436 - 0.3641463957319784im      -0.1916144847277488 + 0.612382137288214im
  -2.827665362682839 + 0.24450099043043694im       6.220001974398564 - 1.2798980449578659im
julia> scalarBA = dot(B,A)2.1354058803024554 + 3.896900584767086im
julia> scalarAA = dot(A,A)26.609992940409455 + 0.0im
julia> normA² = norm(A)^226.609992940409455
julia> U,S,Vd = tsvd(A,(1,3),(2,));
julia> @tensor A′[a,b,c] := U[a,c,d]*S[d,e]*Vd[e,b];
julia> A′ ≈ Atrue
julia> permute(A,(1,3),(2,)) ≈ U*S*Vdtrue

However, trying the following

julia> @tensor D[a,b,c,d] := A[a,b,e]*B[d,c,e]ERROR: SpaceMismatch()
julia> @tensor d = A[a,b,c]*A[a,b,c]ERROR: SpaceMismatch()

we obtain SpaceMismatch errors. The reason for this is that, with ComplexSpace, an index in a space ℂ^n can only be contracted with an index in the dual space dual(ℂ^n) == (ℂ^n)'. Because of the complex Euclidean inner product, the dual space is equivalent to the complex conjugate space, but not the the space itself.

julia> dual(ℂ^3) == conj(ℂ^3) == (ℂ^3)'true
julia> (ℂ^3)' == ℂ^3false
julia> @tensor d = conj(A[a,b,c])*A[a,b,c]26.609992940409455 + 0.0im
julia> d ≈ normA²true

This might seem overly strict or puristic, but we believe that it can help to catch errors, e.g. unintended contractions. In particular, contracting two indices both living in ℂ^n would represent an operation that is not invariant under arbitrary unitary basis changes.

It also makes clear the isomorphism between linear maps ℂ^n → ℂ^m and tensors in ℂ^m ⊗ (ℂ^n)':

julia> m = TensorMap(randn, ComplexF64, ℂ^3, ℂ^4)TensorMap(ProductSpace(ℂ^3) ← ProductSpace(ℂ^4)):
 -0.3859774622263471 - 0.007224842435603468im  …   -1.01277410070073 + 0.3779352074338083im
 -0.3550849309179006 - 0.2595794780600099im       1.5904800828081487 - 0.5071728841143391im
 0.42210303890259426 + 1.1271794560552317im       1.7356387495627128 - 0.8308269760847476im
julia> m2 = permute(m, (1,2), ())TensorMap((ℂ^3 ⊗ (ℂ^4)') ← ProductSpace{ComplexSpace, 0}()):
 -0.3859774622263471 - 0.007224842435603468im  …   -1.01277410070073 + 0.3779352074338083im
 -0.3550849309179006 - 0.2595794780600099im       1.5904800828081487 - 0.5071728841143391im
 0.42210303890259426 + 1.1271794560552317im       1.7356387495627128 - 0.8308269760847476im
julia> codomain(m2)(ℂ^3 ⊗ (ℂ^4)')
julia> space(m, 1)ℂ^3
julia> space(m, 2)(ℂ^4)'

Hence, spaces become their corresponding dual space if they are 'permuted' from the domain to the codomain or vice versa. Also, spaces in the domain are reported as their dual when probing them with space(A, i). Generalizing matrix vector and matrix matrix multiplication to arbitrary tensor contractions require that the two indices to be contracted have spaces which are each others dual. Knowing this, all the other functionality of tensors with CartesianSpace indices remains the same for tensors with ComplexSpace indices.

Symmetries

So far, the functionality that we have illustrated seems to be just a wrapper around dense multidimensional arrays, e.g. Julia's Base Array. More power becomes visible when involving symmetries. With symmetries, we imply that there is some symmetry action defined on every vector space associated with each of the indices of a TensorMap, and the TensorMap is then required to be equivariant, i.e. it acts as an intertwiner between the tensor product representation on the domain and that on the codomain. By Schur's lemma, this means that the tensor is block diagonal in some basis corresponding to the irreducible representations that can be coupled to by combining the different representations on the different spaces in the domain or codomain. For Abelian symmetries, this does not require a basis change and it just imposes that the tensor has some block sparsity. Let's clarify all of this with some examples.

We start with a simple $ℤ₂$ symmetry:

julia> V1 = ℤ₂Space(0=>3,1=>2)Rep[ℤ₂](0=>3, 1=>2)
julia> dim(V1)5
julia> V2 = ℤ₂Space(0=>1,1=>1)Rep[ℤ₂](0=>1, 1=>1)
julia> dim(V2)2
julia> A = Tensor(randn, V1*V1*V2')TensorMap((Rep[ℤ₂](0=>3, 1=>2) ⊗ Rep[ℤ₂](0=>3, 1=>2) ⊗ Rep[ℤ₂](0=>1, 1=>1)') ← ProductSpace{GradedSpace{ZNIrrep{2},Tuple{Int64,Int64}}, 0}()):
* Data for sector (Irrep[ℤ₂](1), Irrep[ℤ₂](1), Irrep[ℤ₂](0)) ← ():
[:, :, 1] =
  1.9524482123604343  -1.1025215403769222
 -0.8254262318339333  -1.5260135932702945
* Data for sector (Irrep[ℤ₂](0), Irrep[ℤ₂](1), Irrep[ℤ₂](1)) ← ():
[:, :, 1] =
 1.0325098716240877    0.019810372415070165
 0.26465816706362744  -0.22366055243340616
 0.16471880106280667   1.6809818250056574
* Data for sector (Irrep[ℤ₂](0), Irrep[ℤ₂](0), Irrep[ℤ₂](0)) ← ():
[:, :, 1] =
 0.27721256236870445  0.6506194726164967   -0.25978844599649953
 0.08534487596800691  0.40554168137394964   0.08922862230746569
 0.3871662871585982   0.5645977174422583    0.04284128761489019
* Data for sector (Irrep[ℤ₂](1), Irrep[ℤ₂](0), Irrep[ℤ₂](1)) ← ():
[:, :, 1] =
 2.2063643639462627  -0.5253273545171222   -0.04713971497989311
 0.3827843131592129   0.06647428107053462   0.9682632751642726
julia> convert(Array, A)5×5×2 Array{Float64,3}:
[:, :, 1] =
 0.277213   0.650619  -0.259788    0.0        0.0
 0.0853449  0.405542   0.0892286   0.0        0.0
 0.387166   0.564598   0.0428413   0.0        0.0
 0.0        0.0        0.0         1.95245   -1.10252
 0.0        0.0        0.0        -0.825426  -1.52601

[:, :, 2] =
 0.0        0.0         0.0        1.03251    0.0198104
 0.0        0.0         0.0        0.264658  -0.223661
 0.0        0.0         0.0        0.164719   1.68098
 2.20636   -0.525327   -0.0471397  0.0        0.0
 0.382784   0.0664743   0.968263   0.0        0.0

Here, we create a space 5-dimensional space V1, which has a three-dimensional subspace associated with charge 0 (the trivial irrep of $ℤ₂$) and a two-dimensional subspace with charge 1 (the non-trivial irrep). Similar for V2, where both subspaces are one- dimensional. Representing the tensor as a dense Array, we see that it is zero in those regions where the charges don't add to zero (modulo 2). The Tensor(Map) type in TensorLabXD.jl won't store these zero blocks, and only stores the non-zero information, which we can recognize in the full Array representation.

From there on, the resulting tensors support all of the same operations as the ones we encountered in the previous examples.

julia> B = Tensor(randn, V1'*V1*V2);
julia> @tensor C[a,b] := A[a,c,d]*B[c,b,d]TensorMap((Rep[ℤ₂](0=>3, 1=>2) ⊗ Rep[ℤ₂](0=>3, 1=>2)) ← ProductSpace{GradedSpace{ZNIrrep{2},Tuple{Int64,Int64}}, 0}()):
* Data for sector (Irrep[ℤ₂](1), Irrep[ℤ₂](1)) ← ():
 2.2185315361612723  -3.0278014305191423
 0.9645250390493137  -0.07204781990588607
* Data for sector (Irrep[ℤ₂](0), Irrep[ℤ₂](0)) ← ():
  0.430601370381284   1.0456786392403254  -0.2832174638270491
 -0.5560454897091337  0.4971552572751569   0.4530555643767742
  1.8242251340646758  2.586203937646578    0.11809779375733462
julia> U,S,V = tsvd(A,(1,3),(2,));
julia> U'*U # should be the identity on the corresponding domain = codomainTensorMap(ProductSpace(Rep[ℤ₂](0=>3, 1=>2)') ← ProductSpace(Rep[ℤ₂](0=>3, 1=>2)')):
* Data for sector (Irrep[ℤ₂](0),) ← (Irrep[ℤ₂](0),):
  1.0                     -4.404287470696228e-17   5.3287027541852015e-17
 -4.404287470696228e-17    1.0000000000000007      2.8718922486671986e-16
  5.3287027541852015e-17   2.8718922486671986e-16  1.0000000000000007
* Data for sector (Irrep[ℤ₂](1),) ← (Irrep[ℤ₂](1),):
 0.9999999999999998     2.363512647801699e-16
 2.363512647801699e-16  0.9999999999999997
julia> U'*U ≈ one(U'*U)true
julia> P = U*U' # should be a projectorTensorMap((Rep[ℤ₂](0=>3, 1=>2) ⊗ Rep[ℤ₂](0=>1, 1=>1)') ← (Rep[ℤ₂](0=>3, 1=>2) ⊗ Rep[ℤ₂](0=>1, 1=>1)')):
* Data for sector (Irrep[ℤ₂](0), Irrep[ℤ₂](0)) ← (Irrep[ℤ₂](0), Irrep[ℤ₂](0)):
[:, :, 1, 1] =
 0.5233439358015453
 0.24408928693144416
 0.4020695085073358

[:, :, 2, 1] =
 0.24408928693144416
 0.16701250044681235
 0.24296910314052214

[:, :, 3, 1] =
 0.4020695085073358
 0.24296910314052214
 0.37626866294733025
* Data for sector (Irrep[ℤ₂](1), Irrep[ℤ₂](1)) ← (Irrep[ℤ₂](0), Irrep[ℤ₂](0)):
[:, :, 1, 1] =
  0.008992844983991234
 -0.16773407445822175

[:, :, 2, 1] =
 -0.06336172342562373
  0.1284174912939829

[:, :, 3, 1] =
 0.029665725028167368
 0.11452784044707999
* Data for sector (Irrep[ℤ₂](0), Irrep[ℤ₂](0)) ← (Irrep[ℤ₂](1), Irrep[ℤ₂](1)):
[:, :, 1, 1] =
  0.008992844983991234
 -0.06336172342562373
  0.029665725028167368

[:, :, 2, 1] =
 -0.16773407445822175
  0.1284174912939829
  0.11452784044707999
* Data for sector (Irrep[ℤ₂](1), Irrep[ℤ₂](1)) ← (Irrep[ℤ₂](1), Irrep[ℤ₂](1)):
[:, :, 1, 1] =
 0.9949541005439055
 0.006693568327540397

[:, :, 2, 1] =
 0.006693568327540397
 0.9384208002604083
* Data for sector (Irrep[ℤ₂](1), Irrep[ℤ₂](0)) ← (Irrep[ℤ₂](1), Irrep[ℤ₂](0)):
[:, :, 1, 1] =
  0.7949940702065715
 -0.06088765103822543

[:, :, 2, 1] =
 -0.06088765103822543
  0.5348949793055442
* Data for sector (Irrep[ℤ₂](0), Irrep[ℤ₂](1)) ← (Irrep[ℤ₂](1), Irrep[ℤ₂](0)):
[:, :, 1, 1] =
  0.33713206725445405
  0.11801519406778958
 -0.17800442604577993

[:, :, 2, 1] =
 -0.18627908267151655
  0.010713844878355727
 -0.45854154801641883
* Data for sector (Irrep[ℤ₂](1), Irrep[ℤ₂](0)) ← (Irrep[ℤ₂](0), Irrep[ℤ₂](1)):
[:, :, 1, 1] =
  0.33713206725445405
 -0.18627908267151655

[:, :, 2, 1] =
 0.11801519406778958
 0.010713844878355727

[:, :, 3, 1] =
 -0.17800442604577993
 -0.45854154801641883
* Data for sector (Irrep[ℤ₂](0), Irrep[ℤ₂](1)) ← (Irrep[ℤ₂](0), Irrep[ℤ₂](1)):
[:, :, 1, 1] =
 0.19152505464992922
 0.04406905305148063
 0.06740325839462899

[:, :, 2, 1] =
  0.04406905305148063
  0.018254937619999783
 -0.04401409924899038

[:, :, 3, 1] =
  0.06740325839462899
 -0.04401409924899038
  0.4603309582179549
julia> P*P ≈ Ptrue

We also support other abelian symmetries, e.g.

julia> V = U₁Space(0=>2,1=>1,-1=>1)Rep[U₁](0=>2, 1=>1, -1=>1)
julia> dim(V)4
julia> A = TensorMap(randn, V*V, V)TensorMap((Rep[U₁](0=>2, 1=>1, -1=>1) ⊗ Rep[U₁](0=>2, 1=>1, -1=>1)) ← ProductSpace(Rep[U₁](0=>2, 1=>1, -1=>1))):
* Data for sector (Irrep[U₁](1), Irrep[U₁](-1)) ← (Irrep[U₁](0),):
[:, :, 1] =
 -0.7853320601334783

[:, :, 2] =
 -1.1084342559366906
* Data for sector (Irrep[U₁](-1), Irrep[U₁](1)) ← (Irrep[U₁](0),):
[:, :, 1] =
 0.6962407750871137

[:, :, 2] =
 -0.9226046248053619
* Data for sector (Irrep[U₁](0), Irrep[U₁](0)) ← (Irrep[U₁](0),):
[:, :, 1] =
  0.5023235670242905   0.21935756149093294
 -0.7707629856218604  -0.014179035139020968

[:, :, 2] =
 0.6786424588146966   0.9401921579979681
 0.02043744466604632  0.060103329980365154
* Data for sector (Irrep[U₁](1), Irrep[U₁](0)) ← (Irrep[U₁](1),):
[:, :, 1] =
 0.8055419568686918  -0.03837463141170912
* Data for sector (Irrep[U₁](0), Irrep[U₁](1)) ← (Irrep[U₁](1),):
[:, :, 1] =
 -0.6328447530726637
  0.25513365852555797
* Data for sector (Irrep[U₁](-1), Irrep[U₁](0)) ← (Irrep[U₁](-1),):
[:, :, 1] =
 -0.24167456514424587  0.3633740381448547
* Data for sector (Irrep[U₁](0), Irrep[U₁](-1)) ← (Irrep[U₁](-1),):
[:, :, 1] =
 0.3778822203187275
 0.4660190553275556
julia> dim(A)20
julia> convert(Array, A)4×4×4 Array{Float64,3}:
[:, :, 1] =
  0.502324   0.219358  0.0        0.0
 -0.770763  -0.014179  0.0        0.0
  0.0        0.0       0.0       -0.785332
  0.0        0.0       0.696241   0.0

[:, :, 2] =
 0.678642   0.940192    0.0        0.0
 0.0204374  0.0601033   0.0        0.0
 0.0        0.0         0.0       -1.10843
 0.0        0.0        -0.922605   0.0

[:, :, 3] =
 0.0        0.0        -0.632845  0.0
 0.0        0.0         0.255134  0.0
 0.805542  -0.0383746   0.0       0.0
 0.0        0.0         0.0       0.0

[:, :, 4] =
  0.0       0.0       0.0  0.377882
  0.0       0.0       0.0  0.466019
  0.0       0.0       0.0  0.0
 -0.241675  0.363374  0.0  0.0
julia> V = Rep[U₁ × ℤ₂]((0, 0)=>2, (-1, 0)=>1, (1, 1)=>1)Rep[U₁ × ℤ₂]((0, 0)=>2, (-1, 0)=>1, (1, 1)=>1)
julia> dim(V)4
julia> A = TensorMap(randn, V*V, V)TensorMap((Rep[U₁ × ℤ₂]((0, 0)=>2, (-1, 0)=>1, (1, 1)=>1) ⊗ Rep[U₁ × ℤ₂]((0, 0)=>2, (-1, 0)=>1, (1, 1)=>1)) ← ProductSpace(Rep[U₁ × ℤ₂]((0, 0)=>2, (-1, 0)=>1, (1, 1)=>1))):
* Data for sector ((Irrep[U₁](0) ⊠ Irrep[ℤ₂](0)), (Irrep[U₁](0) ⊠ Irrep[ℤ₂](0))) ← ((Irrep[U₁](0) ⊠ Irrep[ℤ₂](0)),):
[:, :, 1] =
 0.5609371166193737   0.7039926552673686
 0.34721320816812107  0.3335986792730675

[:, :, 2] =
 0.10332921153308108  -0.3415868875468802
 0.3456449800355501    0.3154865251454887
* Data for sector ((Irrep[U₁](-1) ⊠ Irrep[ℤ₂](0)), (Irrep[U₁](0) ⊠ Irrep[ℤ₂](0))) ← ((Irrep[U₁](-1) ⊠ Irrep[ℤ₂](0)),):
[:, :, 1] =
 0.3060771018762352  1.2937831623543794
* Data for sector ((Irrep[U₁](0) ⊠ Irrep[ℤ₂](0)), (Irrep[U₁](-1) ⊠ Irrep[ℤ₂](0))) ← ((Irrep[U₁](-1) ⊠ Irrep[ℤ₂](0)),):
[:, :, 1] =
 0.7496540302079333
 1.4226382351875397
* Data for sector ((Irrep[U₁](1) ⊠ Irrep[ℤ₂](1)), (Irrep[U₁](0) ⊠ Irrep[ℤ₂](0))) ← ((Irrep[U₁](1) ⊠ Irrep[ℤ₂](1)),):
[:, :, 1] =
 2.5987484323439047  -0.5441914220152319
* Data for sector ((Irrep[U₁](0) ⊠ Irrep[ℤ₂](0)), (Irrep[U₁](1) ⊠ Irrep[ℤ₂](1))) ← ((Irrep[U₁](1) ⊠ Irrep[ℤ₂](1)),):
[:, :, 1] =
 -0.5977526311681345
 -0.21870935540072767
julia> dim(A)16
julia> convert(Array, A)4×4×4 Array{Float64,3}:
[:, :, 1] =
 0.560937  0.703993  0.0  0.0
 0.347213  0.333599  0.0  0.0
 0.0       0.0       0.0  0.0
 0.0       0.0       0.0  0.0

[:, :, 2] =
 0.103329  -0.341587  0.0  0.0
 0.345645   0.315487  0.0  0.0
 0.0        0.0       0.0  0.0
 0.0        0.0       0.0  0.0

[:, :, 3] =
 0.0       0.0      0.749654  0.0
 0.0       0.0      1.42264   0.0
 0.306077  1.29378  0.0       0.0
 0.0       0.0      0.0       0.0

[:, :, 4] =
 0.0       0.0       0.0  -0.597753
 0.0       0.0       0.0  -0.218709
 0.0       0.0       0.0   0.0
 2.59875  -0.544191  0.0   0.0

Here, the dim of a TensorMap returns the number of linearly independent components, i.e. the number of non-zero entries in the case of an abelian symmetry. Note that we can use × (obtained as \times+TAB) to combine different symmetries.

The general space associated with symmetries is a GradedSpace. The concrete type of GradedSpace can be obtained as Vect[I], or if I == Irrep[G] for some G<:Group, as Rep[G]. The ℤ₂Space (or Z2Space) and U₁Space (or U1Space) are just convenient synonyms, e.g.

julia> Rep[U₁](0=>3,1=>2,-1=>1) == U1Space(-1=>1,1=>2,0=>3)true
julia> V = U₁Space(1=>2,0=>3,-1=>1)Rep[U₁](0=>3, 1=>2, -1=>1)
julia> for s in sectors(V)
         @show s, dim(V, s)
       end(s, dim(V, s)) = (Irrep[U₁](0), 3)
(s, dim(V, s)) = (Irrep[U₁](1), 2)
(s, dim(V, s)) = (Irrep[U₁](-1), 1)
julia> U₁Space(-1=>1,0=>3,1=>2) == GradedSpace(Irrep[U₁](1)=>2, Irrep[U₁](0)=>3, Irrep[U₁](-1)=>1)true
julia> supertype(GradedSpace)EuclideanSpace{ℂ}

Generally, GradedSpace supports a grading that is derived from the fusion ring of a (unitary) pre-fusion category. The order in which the charges and their corresponding subspace dimensionality are specified is irrelevant. We can probe the subspace dimension of a certain sector s in a space V with dim(V, s).

The GradedSpace is a subtype of EuclideanSpace{ℂ}, i.e., it has the standard Euclidean inner product and we assume all representations to be unitary.

TensorLabXD.jl also allows for non-abelian symmetries such as SU₂. In this case, the vector space is characterized via the spin quantum number (i.e. the irrep label of SU₂) for each of its subspaces, and is created using SU₂Space (or SU2Space or Rep[SU₂])

julia> V = SU₂Space(0=>2,1/2=>1,1=>1)Rep[SU₂](0=>2, 1/2=>1, 1=>1)
julia> dim(V)7
julia> V == Rep[SU₂](0=>2, 1=>1, 1//2=>1)true

V has a two-dimensional subspace with spin zero, and two one-dimensional subspaces with spin 1/2 and spin 1. A subspace with spin j has an additional 2j+1 dimensional degeneracy on which the irreducible representation acts. This brings the total dimension to 2*1 + 1*2 + 1*3. Creating a tensor with SU₂ symmetry yields

julia> A = TensorMap(randn, V*V, V)TensorMap((Rep[SU₂](0=>2, 1/2=>1, 1=>1) ⊗ Rep[SU₂](0=>2, 1/2=>1, 1=>1)) ← ProductSpace(Rep[SU₂](0=>2, 1/2=>1, 1=>1))):
* Data for fusiontree FusionTree{Irrep[SU₂]}((1/2, 1/2), 0, (false, false), ()) ← FusionTree{Irrep[SU₂]}((0,), 0, (false,), ()):
[:, :, 1] =
 -0.47025617648598567

[:, :, 2] =
 0.20124762449649497
* Data for fusiontree FusionTree{Irrep[SU₂]}((1, 1), 0, (false, false), ()) ← FusionTree{Irrep[SU₂]}((0,), 0, (false,), ()):
[:, :, 1] =
 -1.6572201750296665

[:, :, 2] =
 0.03278245840893683
* Data for fusiontree FusionTree{Irrep[SU₂]}((0, 0), 0, (false, false), ()) ← FusionTree{Irrep[SU₂]}((0,), 0, (false,), ()):
[:, :, 1] =
  1.0139156165927035  -0.2525719982005983
 -0.5773645584479153  -0.4635440009245654

[:, :, 2] =
 -1.1415824679289446   -0.9362196431651272
 -0.12407643282531783   0.005293262286217654
* Data for fusiontree FusionTree{Irrep[SU₂]}((1/2, 0), 1/2, (false, false), ()) ← FusionTree{Irrep[SU₂]}((1/2,), 1/2, (false,), ()):
[:, :, 1] =
 -0.29300477776549266  -0.9511456138756826
* Data for fusiontree FusionTree{Irrep[SU₂]}((1/2, 1), 1/2, (false, false), ()) ← FusionTree{Irrep[SU₂]}((1/2,), 1/2, (false,), ()):
[:, :, 1] =
 0.42649192460510976
* Data for fusiontree FusionTree{Irrep[SU₂]}((0, 1/2), 1/2, (false, false), ()) ← FusionTree{Irrep[SU₂]}((1/2,), 1/2, (false,), ()):
[:, :, 1] =
 -0.6012338018498872
  0.45821972748167505
* Data for fusiontree FusionTree{Irrep[SU₂]}((1, 1/2), 1/2, (false, false), ()) ← FusionTree{Irrep[SU₂]}((1/2,), 1/2, (false,), ()):
[:, :, 1] =
 1.4576765898085842
* Data for fusiontree FusionTree{Irrep[SU₂]}((1, 0), 1, (false, false), ()) ← FusionTree{Irrep[SU₂]}((1,), 1, (false,), ()):
[:, :, 1] =
 0.3051908606526768  -0.08262342597048705
* Data for fusiontree FusionTree{Irrep[SU₂]}((0, 1), 1, (false, false), ()) ← FusionTree{Irrep[SU₂]}((1,), 1, (false,), ()):
[:, :, 1] =
 -0.2732108647893583
 -0.7079743661334351
* Data for fusiontree FusionTree{Irrep[SU₂]}((1, 1), 1, (false, false), ()) ← FusionTree{Irrep[SU₂]}((1,), 1, (false,), ()):
[:, :, 1] =
 1.9999827599340714
* Data for fusiontree FusionTree{Irrep[SU₂]}((1/2, 1/2), 1, (false, false), ()) ← FusionTree{Irrep[SU₂]}((1,), 1, (false,), ()):
[:, :, 1] =
 -0.15640845467119363
julia> dim(A)24
julia> convert(Array, A)7×7×7 Array{Float64,3}:
[:, :, 1] =
  1.01392   -0.252572  0.0        0.0        0.0       0.0        0.0
 -0.577365  -0.463544  0.0        0.0        0.0       0.0        0.0
  0.0        0.0       0.0       -0.332521   0.0       0.0        0.0
  0.0        0.0       0.332521   0.0        0.0       0.0        0.0
  0.0        0.0       0.0        0.0        0.0       0.0       -0.956797
  0.0        0.0       0.0        0.0        0.0       0.956797   0.0
  0.0        0.0       0.0        0.0       -0.956797  0.0        0.0

[:, :, 2] =
 -1.14158   -0.93622      0.0       0.0       0.0        0.0       0.0
 -0.124076   0.00529326   0.0       0.0       0.0        0.0       0.0
  0.0        0.0          0.0       0.142304  0.0        0.0       0.0
  0.0        0.0         -0.142304  0.0       0.0        0.0       0.0
  0.0        0.0          0.0       0.0       0.0        0.0       0.018927
  0.0        0.0          0.0       0.0       0.0       -0.018927  0.0
  0.0        0.0          0.0       0.0       0.018927   0.0       0.0

[:, :, 3] =
  0.0        0.0       -0.601234  0.0       0.0       0.0       0.0
  0.0        0.0        0.45822   0.0       0.0       0.0       0.0
 -0.293005  -0.951146   0.0       0.0       0.0       0.246235  0.0
  0.0        0.0        0.0       0.0      -0.348229  0.0       0.0
  0.0        0.0        0.0       1.19019   0.0       0.0       0.0
  0.0        0.0       -0.84159   0.0       0.0       0.0       0.0
  0.0        0.0        0.0       0.0       0.0       0.0       0.0

[:, :, 4] =
  0.0        0.0        0.0      -0.601234  0.0   0.0       0.0
  0.0        0.0        0.0       0.45822   0.0   0.0       0.0
  0.0        0.0        0.0       0.0       0.0   0.0       0.348229
 -0.293005  -0.951146   0.0       0.0       0.0  -0.246235  0.0
  0.0        0.0        0.0       0.0       0.0   0.0       0.0
  0.0        0.0        0.0       0.84159   0.0   0.0       0.0
  0.0        0.0       -1.19019   0.0       0.0   0.0       0.0

[:, :, 5] =
 0.0        0.0         0.0       0.0  -0.273211  0.0     0.0
 0.0        0.0         0.0       0.0  -0.707974  0.0     0.0
 0.0        0.0        -0.156408  0.0   0.0       0.0     0.0
 0.0        0.0         0.0       0.0   0.0       0.0     0.0
 0.305191  -0.0826234   0.0       0.0   0.0       1.4142  0.0
 0.0        0.0         0.0       0.0  -1.4142    0.0     0.0
 0.0        0.0         0.0       0.0   0.0       0.0     0.0

[:, :, 6] =
 0.0        0.0         0.0        0.0        0.0     -0.273211  0.0
 0.0        0.0         0.0        0.0        0.0     -0.707974  0.0
 0.0        0.0         0.0       -0.110597   0.0      0.0       0.0
 0.0        0.0        -0.110597   0.0        0.0      0.0       0.0
 0.0        0.0         0.0        0.0        0.0      0.0       1.4142
 0.305191  -0.0826234   0.0        0.0        0.0      0.0       0.0
 0.0        0.0         0.0        0.0       -1.4142   0.0       0.0

[:, :, 7] =
 0.0        0.0        0.0   0.0       0.0   0.0     -0.273211
 0.0        0.0        0.0   0.0       0.0   0.0     -0.707974
 0.0        0.0        0.0   0.0       0.0   0.0      0.0
 0.0        0.0        0.0  -0.156408  0.0   0.0      0.0
 0.0        0.0        0.0   0.0       0.0   0.0      0.0
 0.0        0.0        0.0   0.0       0.0   0.0      1.4142
 0.305191  -0.0826234  0.0   0.0       0.0  -1.4142   0.0
julia> norm(A) ≈ norm(convert(Array, A))true

In this case, the full Array representation of the tensor has again many zeros, but it is less obvious to recognize the dense blocks, as there are additional zeros and the numbers in the original tensor data do not match with those in the Array. The reason is of course that the original tensor data now needs to be transformed with a construction known as fusion trees, which are made up out of the Clebsch-Gordan coefficients of the group. Indeed, note that the non-zero blocks are also no longer labeled by a list of sectors, but by pair of fusion trees. This will be explained further in the manual. However, the Clebsch-Gordan coefficients of the group are only needed to actually convert a tensor to an Array. For working with tensors with SU₂Space indices, e.g. contracting or factorizing them, the Clebsch-Gordan coefficients are never needed explicitly. Instead, recoupling relations are used to symbolically manipulate the basis of fusion trees, and this only requires what is known as the topological data of the group (or its representation theory).

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. Preliminary support for these generalizations is present in TensorLabXD.jl and will be extended in the near future.

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. Adhering to this interface should yield code and algorithms that are oblivious to the underlying symmetries and can thus work with arbitrary symmetric tensors.