Closed ericphanson closed 4 months ago
@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
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.
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.
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.
This might be too breaking to be feasible, but I think it is worth considering defining
tensor(a, b) = kron(b, a)
(and in generaltensor(args...) = kron(reverse(args)...)
), essentially because the definition ofkron
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 withreshape
,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 https://github.com/qojulia/QuantumOptics.jl/issues/177#issuecomment-369866461).
Let me give some more examples about why it would help. I'll define
kron_rev
here for clarity, and usekron
andkron_rev
.In section 1.1.3 of Watrous's TQI, he defines the
vec
operator:then:
That is, if
⊗
werekron_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:This issue also has this example:
whereas