`embed` does not preserve the basis when a CompositeBasis with just a single site is used

[noh827]

Suppose op is a LazyTensor whose basis is a CompositeBasis with just one site. Then the basis of dense(op) is not the same as the basis of op. Specifically the basis of dense(op) is a single-site non composite basis (e.g., NLevel(N=10) instead of [NLevel(N=10)]). This then causes an "incompatible basis" error in a more complex implementation of dense where it sums over many "densified" operators.

It looks like this is ultimately due to the implementation of the embed method here https://github.com/qojulia/QuantumInterface.jl/blob/main/src/embed_permute.jl#L55 and how the function tensor works for a single-site operator https://github.com/qojulia/QuantumInterface.jl/blob/main/src/tensor.jl#L7C1-L7C7.

Would it be possible to fix these methods so that the basis is preserved under operations like embed and dense?

[noh827]

@Krastanov for visibility

[amilsted]

I think this should be considered a bug: embed(my_basis, ...) should always return an operator with basis my_basis, or at least in a compatible basis.

A few ways we could go about this:

1. have tensor(op) always return a new object in a (single-element) CompositeBasis
2. introduce a coerce_basis(op, b) interface that returns op in the basis b and use it in embed
3. make single-element CompositeBasis compatible with the non-wrapped basis (would require broader changes to basis handling)

PS: There's a related issue where it's possible to construct CompositeBasis backed by a Vector of bases, rather than a Tuple, but embed always gives us a tuple-backed CompositeBasis, which is not compatible with an equivalent vector-backed CompositeBasis. Another reason to maybe overhaul how we handle basis checks (they're not a good match for the type system!).

[amilsted]

I'm thinking adding a coerce_basis or cast_basis is probably the best path. This is an operation users will sometimes want anyway. We can still do 3 if we want, or more generally overhaul the basis system later, and 1 is a breaking change.

[Krastanov]

@akirakyle , pinging you here as this might also inform tangentially some of the design you are doing for #26

[akirakyle]

Commenting here to link relevant issues so we can hopefully concentrate discussion in this issue:

• embed may not embed for a CompositeBasis of one element QuantumOpticsBase.jl#50
• embed(::NotCompositeBasis, 1, operator) QuantumOpticsBase.jl#159 (Handle single basis in embed #24)
• Embedding for SuperOperators QuantumOpticsBase.jl#93
• ptrace is type instable QuantumOpticsBase.jl#68

[akirakyle]

I also want to suggest another two solutions which I think maybe represent the two "opposite ends of design space" where @amilsted's suggestions lie somewhere in between. The first I think is maybe the closest to how the basis system already treats this issue: forbid single subsystem CompositeBasis. This could be accomplished at the interface level at varying levels of "strictness":

• just state this assertion within the documentation
• add a warning message in the CompositeBasis constructor whenever this happens and change downstream code accordingly
• make the CompositeBasis constructor just return the single subsystem basis

This would require changes to LazyTensor as it currently always requires its basis_l and basis_r to be CompositeBasis.

The second would be to enforce that every basis is always wrapped in a CompositeBasis. This could also be done at similar levels of "strictness" as above. This would be by far the biggest breaking change to fix this issue, as many method signatures rely on distinguishing bases and operators defined on multi-subsystems versus single-subsystem based on dispatch on CompositeBasis.

While I'm back to working a lot on the basis system in #40, I haven't been focusing on this issue since I think the other design considerations should take precedence. Given that work is maturing, I have been revisiting this. Looking into it more, I think it would be really good for that work to include fixes for this issue. This is partly because we can consider options which will be more breaking and make them in conjunction with breaking changes made in #40.

[amilsted]

Nice to see work progressing on bases! I would just add that

make the CompositeBasis constructor just return the single subsystem basis

seems a bit un-Julian, in that we should expect that a constructor actually returns the expected type.

My suggestion to have embed return an operator with the same basis type it received as input is similarly motivated by type stability.

I suppose forbidding single-element composite bases is compatible with the above, but it seems a shame for code that uses QO to have to special-case the single-element case.

I prefer the "every basis is a CompositeBasis" to "disallow single-element CompositeBasis", but I think a middle way is probably preferable overall, where we just have a way to see that a single-element composite basis is "compatible" with the basis it contains.

[akirakyle]

seems a bit un-Julian

Yep, to be clear I'm not advocating for this, just trying to explain the design space as I see it.

My suggestion to have embed return an operator with the same basis type it received as input is similarly motivated by type stability.

Yep that sounds sensible and probably what I'll try to do as I keep working on this.

we just have a way to see that a single-element composite basis is "compatible" with the basis it contains.

Yep I agree here, and the update I posted on #40 has the method interface I've been implementing which should hopefully mean we don't have to care what type a basis is to understand its tensor product structure. This would be really helpful for a cleaner way to implement more efficient or complex bases such as in #43.