Consider reversing order of arguments to kron in tensor

[ericphanson]

This might be too breaking to be feasible, but I think it is worth considering defining tensor(a, b) = kron(b, a) (and in general tensor(args...) = kron(reverse(args)...)), essentially because the definition of kron is (I believe) only natural for row-oriented languages like C and python, and not column-oriented languages like Julia/fortran/MATLAB, and does not interact nicely with reshape, permute, etc in a column-oriented language (as e.g. mentioned here).

I will add that reversing the order like I'm advocating is the convention used in QuantumOptics.jl (see also qojulia/QuantumOptics.jl#177 (comment)).

Let me give some more examples about why it would help. I'll define kron_rev here for clarity, and use kron and kron_rev.

kron_rev(a, b) = kron(b, a)

In section 1.1.3 of Watrous's TQI, he defines the vec operator:

then:

N = 5
e(i) = (z = zeros(N); z[i] = 1; z) # quick basis vector
E(i,j) = (z = zeros(N, N); z[i, j] = 1; z)

vec(E(1,2)) == kron(e(1), e(2)) # false
vec(E(1,2)) == kron_rev(e(1), e(2)) # true

That is, if ⊗ were kron_rev, then this identity holds; if not, it does not. This has a bunch of consequences. For example, we have the identity:

and again we need kron_rev for this to hold:

julia> kron(A₀, A₁) * vec(B) ≈ vec(A₀*B*transpose(A₁))
false

julia> kron_rev(A₀, A₁) * vec(B) ≈ vec(A₀*B*transpose(A₁))
true

This issue also has this example:

julia> A = rand(2,2); B = rand(2,2);

julia> AB = reshape(kron(A, B), (2,2,2,2));

julia> all(AB[a1,b1,a2,b2] ≈ A[a1,a2] * B[b1,b2] for a1=1:2, a2=1:2, b1=1:2, b2=1:2)
false

whereas

julia> AB = reshape(kron_rev(A, B), (2,2,2,2));

julia> all(AB[a1,b1,a2,b2] ≈ A[a1,a2] * B[b1,b2] for a1=1:2, a2=1:2, b1=1:2, b2=1:2)
true

[ytdHuang]

@ericphanson
Thank you for the suggestions.
I think one of the goals of this package is to make the syntax and outputs similar to QuTiP.

The Eq. 1.132 in your uploaded screenshot might not be the case for Qobj.
If A0, A1, and B are all Qobj, tensor(A0, A1) is an Operator, but mat2vec(B) is an OperatorKet.
You cannot multiply them.

The normal way we will do in QuTiP is to create the SuperOperator:

using QuantumToolbox
N  = 10
A0 = Qobj(rand(ComplexF64, N, N))
A1 = Qobj(rand(ComplexF64, N, N))
B  = Qobj(rand(ComplexF64, N, N))
sprepost(A0, A1) * mat2vec(B) ≈ mat2vec(A0 * B * A1) # returns true

[albertomercurio]

Hello @ericphanson and thank you for proposing this. As @ytdHuang has already stated, one of the main purposes of this package is to make it as close as possible to QuTiP functionalities. Although your statement makes a lot of sense, I think that we prefer to maintain the standard kron behavior, such that the resulting matrix is equal to that obtained from QuTiP.

Nonetheless, I’m curious to benchmark the ptrace with the reversed kron, like in the QuantumOptics.jl issue that you cited. Indeed, in the current case, we have to use the PermuteDimsArray function to permute them. As far as I know this is performed in a lazy way, and it would not make huge differences. But I’m still curious.

[ericphanson]

Ah, it's nice the Qobj & SuperOperator abstraction allows this to work still. It might be simpler behind-the-scenes if the tensor convention was reversed, but at least it won't affect users too much.

In general it may be not-optimally-performant to use the exact same dimension ordering in python and Julia, because of column-major vs row-major (see also https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-column-major). For example, Flux.jl uses the last dimension as the batch-dimension while pytorch uses the first dimension as the batch dimension, exactly for this reason. So it is common for Julia and python libraries to differ in this way since they each want maximal performance. However I can see the value of keeping the same ordering between the two libraries.

[albertomercurio]

The example that you made between Flux.jl and pytorch is very clear, and I agree with you that maximum performances are more important than the similarity between two languages.

In this case, however, I don't see a significant improvement in terms of performance, introducing a significant difference from QuTiP. If we take for example the ptrace function of both QuTiP and QuantumToolbox.jl, we see that the only difference is the ordering of the dimensions list during reshape and PermuteDimArrays. For example, even for 20 spins (almost within the limit of computational resources of any computer), the only difference would be in reverse the vector containing the dimensions (a vector of size 20), which is not computationally heavy.

You can see the code comparison between the current two methods, between QuTiP and QuantumToolbox.jl.

BTW, I will close this issue because this reverse tensor product is not necessary for the moment. Feel free to reopen it if you find any good reason. I'm always open to new possibilities to make the package more performant.