`IncompatibleBases` as a static versus dynamic check, i.e. should `Basis` be `isbits`?

[akirakyle]

The deeper I dig into the basis system as part of #32, the more I seem to see reasons for an overhaul 😅. This is of course very related to #27 and #33 but different enough that I thought it helpful to have a dedicated issue. Currently it seems there's two methods at tension for checking whether bases are compatible: the samebases machinery versus the typesystem. This is easily illustrated with an example

using QuantumOpticsBase
k1 = basisstate(SpinBasis(1), 1)
k2 = basisstate(SpinBasis(3//2), 2)
k1 + k2

throws ERROR: QuantumInterface.IncompatibleBases()
meanwhile

using QuantumOpticsBase
k1 = basisstate(NLevelBasis(3), 1)
k2 = basisstate(NLevelBasis(4), 2)
k1 + k2

throws ERROR: DimensionMismatch: a has size (2,), b has size (3,), mismatch at dim 1

The relevant part isn't so much the error messages themselves, we can easily (and probably should) make the latter example throw IncompatibleBases, but in how these two situations are caught due to using the type system

+(a::Ket{B}, b::Ket{B}) where {B} = Ket(a.basis, a.data+b.data)
+(a::Ket, b::Ket) = throw(IncompatibleBases())

When differences in bases are encoded within their types as currently happens with SpinBasis when called created with an isbits type, then this pattern catches incompatible bases. However, if incompatible bases are not evident at the type level, for example when using BigInt which has isbits false, then a dynamic samebases check is required.

I think making Basis isbits could potentially to unlock performance benefits since this might allow StaticArrays to be used! Then downside is that this means one has to carefully subtype Basis so that equality of bases is fully encoded in type parameters for which isbits is true. I think this is certainly doable for all the bases collected or proposed in #33. Of course this would also be a major breaking change.

In either case, I think it makes sense to pick only one pattern, document it here, and try to stick to it in dependent packages.

[Krastanov]

"Breaking change" can mean two different things in the Julia ecosystem.

If it is breaking to some other libraries that rely on internal details of QuantumInterface, then that is not "technically/semver" breaking.

If it is breaking the public API, then it is "technically/semver" breaking.

I suspect much of what you are suggesting here is not "technically" breaking because it is about the internal layout of structs. We can still keep their old constructor API even after they are modified to be isbits. Can we do this in small increments? What is the smallest commit that can be made in this direction? Maybe making just one of the basis be isbits? Doing it in such small increments might be possible to do 90% of the work without any breaking.

[akirakyle]

Yes I've been keeping the constructors the same as I've been experimenting on this, that's easy. The question is whether changes in parametric typing is considered breaking? This is also the same question about changing from AbstractOperator{BL,BL} to AbstractOperator{B<:OperatorBasis} since anything that method dispatches based on the parametric types will break.

[Krastanov]

I have not checked too carefully, but I suspect changes in the type parameters is not semver/technically breaking.

Also, as we are going through these changes, especially through fixes, we should see whether we can simplify the dispatch signatures so that they are less susceptible to being broken by such internal changes.

[amilsted]

I have been thinking about related issues for a while and I would like to offer my two cents here. There is a place where relying on the type system for basis-equality checks is, in my opinion, problematic: CompositeBasis.

A lot of functionality relies on type checking to exclude invalid basis combinations. This sometimes misses incompatibilities (as you describe above). It is also sometimes too strict (as in #27) and, in the case of CompositeBasis in particular, leads to a combinatorial explosion of possible basis types, because the concrete type of a CompositeBasis includes a tuple of the component basis: tensor(FockBasis(d), GenericBasis(2)) gives us a basis with a different type to tensor(GenericBasis(2), FockBasis(d)) for example.

Of course, this is by design since we rely on types to test basis equivalence, but it can lead to excessive compilation overhead. Consider the basis reduce(tensor, (GenericBasis(2) for _ in 1:N)) - we will end up compiling any high-level function for every value of N encountered (and even for fixed N will probably compile a version of tensor on every application).

This makes me consider whether we should rely much less on the type system for basis checking and have runtime (even runtime-dispatched) basis checking by default. We can choose to turn off these checks (e.g. via a global variable) for things like integration loops in time evolution. Interested to hear your thoughts @akirakyle and @Krastanov.

One strategy might be to introduce runtime checks for CompositeBasis in particular (and maybe transition CompositeBasis to use a vector of component bases, rather than a tuple).

[Krastanov]

I am also leaning towards using much more runtime tools. The incompatible basis checks are not needed in hot loops, so there is no engineering need to make them static.

[akirakyle]

In terms of pure programming style I think the using the type system for basis checks is much more elegant. As for the performance cost (https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type), I'm not sure I necessarily agree that the overhead is unmanageable since it really depends on the usage patterns? Ideally both approaches could be benchmarked, but that requires a bit of work. I think one thing to keep in mind is that the number of elements in a CompositeBasis, at least for anything using explicit data types, is practically quite limited given that the total Hilbert dimension grows exponentially.

and maybe transition CompositeBasis to use a vector of component bases, rather than a tuple

I'm curious why you think a vector is preferable over a tuple? I would think that regardless of how checks are implemented, we should treat bases as immutable objects?

[amilsted]

The reason to use a vector is again to reduce "type load". The type of a tuple is specialized to the types of its entries, whereas a vector can have a generic element type. Any function that takes the tuple as an argument must be compiled again for each basis combination.

Regarding your other point, it is quite possible to have say, 20-30 element composite bases. Say you only allow your local basis to either be FockBasis(1) or GenericBasis(2). That's still millions of possible combinations. Of course, we won't likely be hitting them all in practice, but it still seems wasteful to a new function for each.

The compilation overhead can be significant because the bases pop up all over QO in type parameters. Any function that takes an operator or state argument specializes on the basis type, since the basis type is a type parameter of the state or operator, and hence needs to be recompiled for every new basis type encountered. Even might be worth removing the basis type parameters on operator and state types for this reason.

[Krastanov]

Here is an example where a dynamic setup would be faster than the static type check:

for operator in large_list_of_operators
    apply!(state, operator, affected_qubits(operator))
end

If large_list_of_operators has a ton of different isbits basis types, this will break Julia's union splitting and will require a runtime dispatch for each apply! operation.

If we use more dynamism in the basis checks, then large_list_of_operators will not be as diversely typed and everything will be statically dispatched.

This issue has been coming up often in circuit simulations with QuantumClifford, and is something I need from QuantumOptics, but I have not made it possible yet.

[akirakyle]

Alright I think I'm convinced that at least the potential explosion in types with making bases isbits isn't going to be worth it for how bad the performance pitfalls could be. As an experiment I tried

abstract type Basis end

struct GenericBasis{N} <: Basis
    GenericBasis(N) = new{N}()
end

struct CompositeBasis{T<:Tuple{Vararg{Basis}}} <: Basis
    bases
    CompositeBasis(x) = new{typeof(x)}(x)
end

bases(b::CompositeBasis) = b.bases
nsubsystems(b::Basis) = 1
nsubsystems(b::CompositeBasis) = length(bases(b))

function reduced(b::CompositeBasis, indices)
    if length(indices)==0
        throw(ArgumentError("At least one subsystem must be specified in reduced."))
    elseif length(indices)==1
        return bases(b)[indices[1]]
    else
        return CompositeBasis(bases(b)[indices])
    end
end

function ptrace(b::CompositeBasis, indices)
    J = [i for i in 1:nsubsystems(b) if i ∉ indices]
    length(J) > 0 || throw(ArgumentError("Tracing over all indices is not allowed in ptrace."))
    reduced(b, J)
end

function timeit(bs)
    @time r1 = [ptrace(b, [i for i=1:nsubsystems(b) if i%3==0]) for b in bs];
    @time r2 = [ptrace(b, [i for i=1:nsubsystems(b) if i%3==0]) for b in bs];
    return r1,r2
end

timeit([CompositeBasis(Tuple(GenericBasis(i+j) for j=1:40)) for i=1:100]);

giving

  2.129629 seconds (4.18 M allocations: 200.085 MiB, 1.23% gc time, 99.59% compilation time)
  0.349378 seconds (232.07 k allocations: 11.013 MiB, 99.29% compilation time)

If instead I replace the structs with

struct GenericBasis <: Basis
    N
end

struct CompositeBasis <: Basis
    bases
end

and run the benchmarks as

println("small tensors - tuple")
timeit([CompositeBasis(Tuple(GenericBasis(i+j) for j=1:40)) for i=1:1000]);
println("small tensors - vector")
timeit([CompositeBasis([GenericBasis(i+j) for j=1:40]) for i=1:1000]);

println("big tensors - tuple")
timeit([CompositeBasis(Tuple(GenericBasis(i+j) for j=1:100)) for i=1:1000]);
println("big tensors - vector")
timeit([CompositeBasis([GenericBasis(i+j) for j=1:100]) for i=1:1000]);

I get

small tensors - tuple
  0.080450 seconds (223.22 k allocations: 10.906 MiB, 96.35% compilation time)
  0.029728 seconds (96.21 k allocations: 4.518 MiB, 91.00% compilation time)
small tensors - vector
  0.038910 seconds (65.98 k allocations: 3.877 MiB, 32.06% gc time, 96.60% compilation time)
  0.001320 seconds (20.00 k allocations: 1.518 MiB)
big tensors - tuple
  0.013207 seconds (96.29 k allocations: 5.158 MiB, 58.38% compilation time)
  0.005447 seconds (89.00 k allocations: 4.799 MiB)
big tensors - vector
  0.002483 seconds (21.00 k allocations: 3.258 MiB)
  0.002457 seconds (21.00 k allocations: 3.258 MiB)

So I think just by not parametrically typing subtypes of Basis we'll avoid the issue with an explosion of compiling specalized functions. While the vector implementation is slightly faster, I would worry about the fact that they can be unintentionally mutated which could cause bugs when the mutated bases field no longer corresponds to whatever object it is attached to (i.e. Ket, Bra, Operator).

[Krastanov]

This is a really neat set of benchmarks and experiments, thank you!

You are right about all the worries about mutation and generally the fact that runtime checks are less elegant. This is a frequent seen tradeoff between perfection and practicality.

Julia usually lets people mutate private fields on the principle "we are all adults here". As long as these fields are not exposed as a public part of the API, I think the risk level of someone messing them up and mutating them is acceptable.

[akirakyle]

Julia usually lets people mutate private fields on the principle "we are all adults here". As long as these fields are not exposed as a public part of the API, I think the risk level of someone messing them up and mutating them is acceptable.

Got it, vectors it is then 😅

I think this further means that for StateVector, AbstractKet, AbstractBra, AbstractOperator, or AbstractSuperOperator (I'll call all these quantum objects) we should not give them a type parameter constrained to be a subtype of basis since we don't want functions which accept their subtypes to be specalized based on the subtype of Basis they're defined on. Rather we just require them to implement basis and fullbasis methods.

While this change will require a fair bit of updating in QuantumOpticsBase because there's a lot of methods which use the (a::Ket{B}, b::Ket{B}) where {B} pattern, it's definitely going to be good to have a well defined pattern for knowing when objects are compatible. That brings up the next issue which is what that interface should look like. Currently here we have

• multiplicable(a,b) with docstring stating "Check if two objects are multiplicable." and defined on both subtypes of Basis and subtypes of quantum objects.
• samebases with docstring stating "Test if one object has the same left and right bases or if two objects have the same bases" and defined for one and two arguments of subtypes of bases and one and two arguments of subtypes of quantum objects.
• check_multiplicable and check_samebases which call their respective functions after first testing BASES_CHECK which can be disabled the @samebases macro.

As far as naming go, I personally find it quite confusing that the @samebases actually disables calls to samebases and multiplicable happening within that call stack. I also think allowing these functions to be called on both subtypes of basis and subtypes of quantum objects is a bit confusing, for example what does it mean to multiply two bases? I think practically there's several distinct patterns that need checks: if quantum objects can be added, or if quantum objects can be multiplied, and if operators or superoperators have the same left and right bases. Thus I might suggest clarifying the interface to something like

• 'multiplicable(a,b)' where a,b are subtypes of quantum objects and checks that they can be multiplied (right_basis(a) == left_basis(b))
• addible(a,b) where a,b are subtypes of quantum objects and checks that they can be added (left_basis(a) == left_basis(b) && right_basis(a) == right_basis(b))
• issquare(a) where a is a subtype of a quantum object and checks that it is square (left_basis(a) == right_basis(b))
• check_* versions of the above with a @compatiblebases macro which disables them.
• For checking if bases are the same, I think we should just explicitly use == rather than a samebases since I think that's clerarer anyways. Also I think this would allow a slower transition to the new interface as samebases could be kept the same but eventually deprecated once everything is using a "finer grained" basis checking interface.

Thoughts?