First draft of FusionTensor

[ogauthe]

This PR provides a first draft for FusionTensors. It defines a generic interface to implement non-abelian symmetries. It currently supports O(2), SU(2), any abelian group defined in SymmetrySectors. It allows to construct, permute and contract tensors.

Work to be done in future PR:

• remove FusedAxes, currently only used in FusionTensor initialization
• handle dual in a more explicit way when fusing codomain and domain fusion trees
• store unitary as a sparse matrix
• improve caching of structural data
• maybe improve detection of a charge block size
• maybe moving Clebsch-Gordan tensor definition to SymmetrySectors

Work depending on other packages

• implement linear algebra decomposition such as SVD (depends on [Enhancement] Add singular value decomposition BlockSparseArrays.jl#16)
• rewrite TensorAlgebra interface once contract preserves domain and codomain information
• move nsymbol to SymmetrySectors (define nsymbol SymmetrySectors.jl#4)

[ogauthe]

I have issues with building test both on local and on GitHub actions. Locally, I get an error

  LoadError: ArgumentError: Package LinearAlgebra not found in current path.

although the package is included in the Project.toml. I guess there was an issue with SafeTestsets configuration when I merged the new skeleton. @lkdvos do you have a clue?

[ogauthe]

So I added test dependencies and I am now able to pass tests locally with (FusionTensors) pkg> test FusionTensors. However they still fail on GitHub actions. Am I missing something? I cannot find any clear error message. I wonder if I need to specify somewhere to get packages from ITensorRegistry?

[mtfishman]

So I added test dependencies and I am now able to pass tests locally with (FusionTensors) pkg> test FusionTensors. However they still fail on GitHub actions. Am I missing something? I cannot find any clear error message. I wonder if I need to specify somewhere to get packages from ITensorRegistry?

I fixed it, the URL for BlockSparseArrays.jl was wrong in the ITensorRegistry. Now the only failure is an Aqua failure in Julia 1.10, we could just disable that for now.

[mtfishman]

A naming convention I have been using is defining functions that provide default values for a keyword argument k as default_k(...), so here it would be default_tol(...).

Also, here's a suggestion for making it a bit more compact and general:

default_tol(a::AbstractArray) = default_tol(eltype(a))
default_tol(arrayt::Type{<:AbstractArray}) = default_tol(eltype(arrayt))
default_tol(elt::Type{<:Number}) = 10 * eps(real(float(elt)))

[lkdvos]

It might be nice to specify atol or rtol too, as this can become relevant for SU(N), or even larger irreps of SU(2), where the CGCs can become tiny

[ogauthe]

Base Julia uses rtoldefault, looks good to me. I also set atol to 0, could be higher.

[lkdvos]

Given that this is adapted from TensorKit anyways, how would you feel about:

function clebsch_gordan_tensor(s1, s2, s3)
    s1_tk, s2_tk, s3_tk = convert.(TensorKitSectors.Sector, (s1, s2, s3))
    cgt = TensorKitSectors.fusiontensor(s1_tk, s2_tk, s3_tk)
    return reshape(cgt, size(cgt, 1), size(cgt, 2), size(cgt, 3)) # discard trailing singleton dimension
end   

I don't think there is a lot of value in having the implementation copied here, it's definitely more maintenance and error-prone, so it would be cool if we could immediately share this.
(This is a registered light-weight package that only contains the symmetry data, so the dependency is quite light, and this automatically inherits many of the symmetries)
https://github.com/QuantumKitHub/TensorKitSectors.jl

[mtfishman]

That sounds reasonable to me.

[lkdvos]

If you use symbol_1j, it might be more internally consistent to then also use the symbol_3j and symbol_6j notations instead of clebsch_gordan_tensor? In general, clebsch gordan seems to be very often related specifically to SU(N) anyways.

[ogauthe]

I do not have a strong opinion here. We can have both symbol_1/3/6j and flipper, clebsch_gordan and f_symbol interface, with one implementation and one converter.

[lkdvos]

I think it would be nice to be consistent with the fieldnames: either add _sector(s) everywhere or nowhere? Given that sector is already in the typename, it might not give that much additional information to add it in the fields as well, so it could maybe be leaves, arrows, root, branches and vertices/multiplicities/... I usually get a bit confused by outer_multiplicity_indices because these don't actually live on the "outer" charges, but only inside of the tree. Ultimately I don't feel too strongly about this though, just sharing my opinion.

[mtfishman]

Those changes all sound reasonable to me. I agree the name outer_multiplicity_indices is a bit confusing.

[ogauthe]

The problem here is we are using outer for two different meaning. We use it to label the outer axes, the ones we fuse with the fusion tree. However in outer_multiplicity_indices, outer refers to the fusion ring outer multiplicity, a well defined concept widely used in the literature. Since we plan to use nsymbol for this, maybe we should use nsymbol_indices?

[lkdvos]

If deciding to keep this, definitely should make sure that this is properly documented :)

[lkdvos]

probably don't want to hard-code that here, since this is not always the case (could also be Complex, and for abelian symmetries it is even convenient to use Bool). eltype(fusiontensortype) or something similar might work?

[ogauthe]

I propose fusiontree_eltype(::Type<:AbstractSector}) as an interface. Currently it always return Float64, we may use another type for abelian groups.

[mtfishman]

Related to the discussion elsewhere, probably we should redefine length(::AbstractGradedUnitRange) to output the number of basis vectors in the vector space as opposed to the number of independent parameters, so then this can be defined as: Base.size(ft::FusionTensor) = length.(axes(ft)) (though I think we won't need to define this overload since Base defines size automatically for AbstractArrays as long as you define axes).

[ogauthe]

From reading the comments above, I see two problems with the current design:

• length(::AbstractGradedUnitRange) does not match the size of the fusion tensor
• data_matrix does not make explicitly the pairing between a sector and a matrix block

Independently, there is a question whether axes(::FusionTensor) should return a flat tuple or a tuple of 2 tuples.

I think length(::AbstractGradedUnitRange) should return the quantum dimension and not the number of multiplets. This means rewriting some parts of GradedUnitRanges, hopefully the change can be easily implemented now that AbstractGradedUnitRange is no more an alias for BlockedUnitRange. Then we would automatically have length.(axes(ft)) == size(ft). This could also play well with ITensor/SymmetrySectors.jl#4, where we already discussed changing blocklengths(::AbstractGradedUnitRange) behavior.

Concerning data_matrix, to me the cleanest would be to have a custom type to take into account the label. Something like

struct SymmetricDataMatrix <: AbstractBlockSparseMatrix
    matrix_blocks::Vector{Pair{AbstractSector,AbstractMatrix}}
    row_axis::AbstractGradedUnitRange
    column_axis::AbstractGradedUnitRange
end

which would act very similarly to the current data_matrix, but the quantum dimension of each block would be easier to access.