Clearly document the basis "endiandness" in tensor products and partial traces (and more) and make a comparison table to other julia and python tools and to kron in basic linear algebra and tensor network toolkits

[matchatealeaf]

I've been struggling with differing results between QuantumOptics.jl and python's QuTiP which features only simple application of unitary operators on composite quantum systems, and I'm not sure if I am misunderstanding the workings of QuantumOptics.jl.

I'm sorry if this is not a minimal example as I've been unable to reproduce the discrepancy reliably (but it's just matrix multiplication, tensor products, and partial trace).

The operation that I'm performing is on a 5 qubits system (labelled ABCDE), after which I want the state of the first qubit (labelled A):

$$ \sigma^{ABCDE} = U \rho^{ABCDE} U^\dagger, $$ $$ \sigma^A = Tr_{BCDE} (\sigma^{ABCDE}) $$ where $$ \rho^{ABCDE} = |+\rangle\langle +| \otimes \rho_N \otimes \rho_N \otimes |0\rangle\langle 0| \otimes \rho_N $$ $$ U = U_2 U_1 $$ $$ U_1 = |0\rangle\langle 0| \otimes U_A \otimes I + |1\rangle\langle 1| \otimes I \otimes U_B $$ $$ U_2 = |1\rangle\langle 1| \otimes U_A \otimes I + |0\rangle\langle 0| \otimes I \otimes U_B $$

The values of $\rho_N$, $U_A$, and $U_B$ are shown in the code below.

using QuantumOptics

b = SpinBasis(1//2)

rho0 = Operator(b, [1 0; 0 0])
rho1 = Operator(b, [0 0; 0 1])
rho_plus = Operator(b, [0.5 0.5; 0.5 0.5])
rho_N = Operator(b, [0.924142 0; 0 0.0758582])

A_mat = [1 0 0 0 0 0 0 0;
         0 1 0 0 0 0 0 0;
	 0 0 1 0 0 0 0 0;
	 0 0 0 0 -1 0 0 0;
	 0 0 0 1 0 0 0 0;
	 0 0 0 0 0 1 0 0;
	 0 0 0 0 0 0 1 0;
	 0 0 0 0 0 0 0 1]
B_mat = [1 0 0 0 0 0 0 0;
	 0 0 0 0 0 0 1 0;
	 0 0 1 0 0 0 0 0;
	 0 0 0 1 0 0 0 0;
	 0 0 0 0 1 0 0 0;
	 0 0 0 0 0 1 0 0;
	 0 -1 0 0 0 0 0 0;
	 0 0 0 0 0 0 0 1]
UA = Operator(tensor(b, b, b), A_mat)
UB = Operator(tensor(b, b, b), B_mat)
	
U1 = tensor(rho0, UA, one(b)) + tensor(rho1, one(b), UB)
U2 = tensor(rho1, UA, one(b)) + tensor(rho0, one(b), UB)
U = U2 * U1

rho = tensor(rho_plus, rho_N, rho_N, rho0, rho_N)
sigma = U * rho * dagger(U)
sigma_A = ptrace(sigma, [2, 3, 4, 5])

sigma_A
Operator(dim=2x2)
  basis: Spin(1/2)
     0.5+0.0im  0.46473+0.0im
     0.46473+0.0im      0.5+0.0im

Now in QuTiP (note that the code are the same except for syntax differences):

import qutip as qt

rho0 = qt.Qobj([[1, 0], [0, 0]])
rho1 = qt.Qobj([[0, 0], [0, 1]])
rho_plus = qt.Qobj([[0.5, 0.5], [0.5, 0.5]])
rho_N = qt.Qobj([[0.924142, 0], [0, 0.0758582]])

A_mat = [[1, 0, 0, 0, 0, 0, 0, 0],
	 [0, 1, 0, 0, 0, 0, 0, 0],
	 [0, 0, 1, 0, 0, 0, 0, 0],
	 [0, 0, 0, 0, -1, 0, 0, 0],
	 [0, 0, 0, 1, 0, 0, 0, 0],
	 [0, 0, 0, 0, 0, 1, 0, 0],
	 [0, 0, 0, 0, 0, 0, 1, 0],
	 [0, 0, 0, 0, 0, 0, 0, 1]]
B_mat = [[1, 0, 0, 0, 0, 0, 0, 0],
	 [0, 0, 0, 0, 0, 0, 1, 0],
	 [0, 0, 1, 0, 0, 0, 0, 0],
	 [0, 0, 0, 1, 0, 0, 0, 0],
	 [0, 0, 0, 0, 1, 0, 0, 0],
	 [0, 0, 0, 0, 0, 1, 0, 0],
	 [0, -1, 0, 0, 0, 0, 0, 0],
	 [0, 0, 0, 0, 0, 0, 0, 1]]
UA = qt.Qobj(A_mat, dims=[[2, 2, 2], [2, 2, 2]])
UB = qt.Qobj(B_mat, dims=[[2, 2, 2], [2, 2, 2]])
	
U1 = qt.tensor(rho0, UA, qt.qeye(2)) + qt.tensor(rho1, qt.qeye(2), UB)
U2 = qt.tensor(rho1, UA, qt.qeye(2)) + qt.tensor(rho0, qt.qeye(2), UB)
U = U2 * U1

rho = qt.tensor(rho_plus, rho_N, rho_N, rho0, rho_N)
sigma = U * rho * U.dag()
sigma_A = qt.ptrace(sigma, 0)  # QuTiP's ptrace specifies the leftover index instead of the indices to trace out

sigma_A
Quantum object: dims=[[2], [2]], shape=(2, 2), type='oper', dtype=Dense, isherm=True
Qobj data =
[[0.5000003  0.43255551]
 [0.43255551 0.5000003 ]]

In QuantumOptics.jl, the cross term is 0.46473, but in QuTiP it's 0.43255551.
Why are they different?

Julia Version:    1.11.1
QuantumOptics:     v1.2.1

qutip.about()
QuTiP Version:      5.0.4
Numpy Version:      2.1.3
Scipy Version:      1.14.1
Cython Version:     None
Matplotlib Version: 3.9.2
Python Version:     3.12.7
Number of CPUs:     12
BLAS Info:          Generic
INTEL MKL Ext:      False
Platform Info:      Linux (x86_64)

[david-pl]

I didn't check, but my guess would be the fact that tensor(x, y) = kron(y, x) in QuantumOptics.jl is messing things up here since you explicitly defined the matrices A and B and they seem to act on more than one Hilbert space.

Can you rewrite things so you use the built-in tensor function to define the unitaries? Otherwise, you can also adapt the definition of the matrices.

[Krastanov]

That is indeed the case. The "endiandness" of how basis vectors are enumerated in various tools is not fixed: numpy, qutip, qiskit, QuantumOptics, etc all can have differing ordering.

Two examples where this comes in:

julia> tensor(Operator(GenericBasis(2), [0 1; 2 3]), Operator(GenericBasis(2), [0 10; 100 1000]))
Operator(dim=4x4)
  basis: [Basis(dim=2) ⊗ Basis(dim=2)]
   0.0    0.0     0.0    10.0
   0.0    0.0    20.0    30.0
   0.0  100.0     0.0  1000.0
 200.0  300.0  2000.0  3000.0

In [6]: tensor(Qobj([[0, 1],[2, 3]]), Qobj([[0, 10],[100, 1000]]))
Out[6]: 
Quantum object: dims=[[2, 2], [2, 2]], shape=(4, 4), type='oper', dtype=Dense, isherm=False
Qobj data =
[[   0.    0.    0.   10.]
 [   0.    0.  100. 1000.]
 [   0.   20.    0.   30.]
 [ 200. 2000.  300. 3000.]]

As a general design principle, all of these tools are set up to not expect the user to input raw matrices (because they are basis dependent), rather to build things up from explicitly defined basis states. In both qutip and QuantumOptics it is most convenient to do use initial defining operations like:

• for basis states: spinup, spindown
• for basic operator building blocks: projector (on the output of spinup and spindown, including the two-argument projector), sigmax, identityoperator, etc.

No matter the tool, in python or in julia, this approach makes things a bit safer.

I will close the issue for now, as it is not a bug (although we can potentially document this more clearly -- PRs welcomed on that ;)

Do not hesitate to post or reopen the issue if anything else is not working as expected.

[matchatealeaf]

Thanks for the clarification. The discrepancy indeed vanishes when I generate the unitaries from their corresponding interaction Hamiltonians instead.

Though this still feels a little misleading, for example one might simply want to use QuantumOptics to compute some tensor product of say

$$ \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} \otimes \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix}= \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}, $$

yet

using QuantumOptics

b = NLevelBasis(2)

tensor(dm(nlevelstate(b, 1)), dm(nlevelstate(b, 2)))
Operator(dim=4x4)
  basis: [NLevel(N=2) ⊗ NLevel(N=2)]
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  1.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im

despite

dm(nlevelstate(b, 1))
Operator(dim=2x2)
  basis: NLevel(N=2)
 1.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im

and

dm(nlevelstate(b, 2))
Operator(dim=2x2)
  basis: NLevel(N=2)
 0.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im

At the very least shouldn't QuantumOpticsBase.set_printing(standard_order=true) be the default?

Likewise, if one simply wants to find the partial trace of some arbitrary operator, then shouldn't standard order also be used for their definition so that we have

$$ Tr_2 \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}= \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} $$

instead of

ptrace(Operator(tensor(b, b), [0 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 0]), 2)
Operator(dim=2x2)
  basis: NLevel(N=2)
 0.0  0.0
 0.0  1.0

While we can define [0 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 0] with basic operator building blocks in QuantumOptics, but consider if it is an arbitrary operator that one came across without knowing how it is generated.

[Krastanov]

Could you elaborate what the problem is in the first example you are giving? That tensor product result seems appropriate to me. Same with the second example of a partial trace.

I can easily agree that this needs to be documented more clearly to avoid confusion, so I will reopen this with a changed title. Maybe a comparison table to all the different conventions in qutip/qiskit/quantumtoolbox/etc would be nice to have.

Concerning your example of finding an arbitrary operator: if one finds a matrix somewhere in the literature, but they do not know in what basis that matrix is written, then it is indeed impossible to use that result in a numerical simulation reliably. The "endiandness" of a basis is not standardized among papers and textbooks (to the extent that researchers in the pedagogy of QIS have to create surveys specifically agnostic to endiandness choices).

And that comes up in classical computing as well: if you find some arbitrary binary data over the network, you do not actually know what are the bytes you have found, as the endiandness of the bytes is not specified -- this comes up all the time in BLE or UART comms.

[matchatealeaf]

Thanks for reopening the issue again.
My main issue is the misleading nature of the "endianness" that QuantumOptics.jl decided on.

Could you elaborate what the problem is in the first example you are giving? That tensor product result seems appropriate to me. Same with the second example of a partial trace.

I expected the tensor product $A \otimes B$ to give kron(A, B)

$$ \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} \otimes \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix}= \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} $$ but it gives $$ \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} \otimes \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix}= \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} $$

kron(B, A) instead.
Perhaps I'm misunderstanding, but we should always have kron(A, B) isn't it? In what basis would the reverse of kron(B, A) be true?

As for the partial trace example, if I have $\rho_1 \otimes \rho_2$, then $\mathrm{Tr}_2(\rho_1 \otimes \rho_2)$ should give me $\rho_1$, not $\rho_2$.

In the example I have

$$ \rho_1 = \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} ,\quad \rho_2 = \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix} $$

and $\rho_1 \otimes \rho_2$ is kron(A, B), but ptrace(kron(A, B), 2) gives me $\rho_2$ instead of $\rho_1$, i.e., it interprets the argument as kron(B, A) instead.

In both examples, I am simply using the common computational basis, i.e., $\rho_1$ is nlevelstate(NLevelBasis(2), 1) and $\rho_2$ is nlevelstate(NLevelBasis(2), 2).
I do not see how the "endianness" of basis would lead to these discrepancies.

Consider if a researcher constructed a unitary operation U from some interaction Hamiltonian in the computational basis, and they want to perform a controlled U operation on a set of input states.
They use QuantumOptics.jl to compute the resulting controlled unitary to put in their paper, i.e., the unitary $|1\rangle\langle 1| \otimes U + |0\rangle\langle 0| \otimes I$.
However, QuantumOptics.jl gives $U \otimes |1\rangle\langle 1| + I \otimes |0\rangle\langle 0|$ instead.
They are simply using the typical computational basis: when they do nlevelstate(NLevelBasis(2), 1), they see the familiar $|0\rangle = (1, 0)^T$ and when they do nlevelstate(NLevelBasis(2), 2), they see $|1\rangle = (0, 1)^T$.

So naturally they use NLevelBasis(2), it's just supposed to be matrix multiplications after all.
They would not suspect such a "gotcha" regarding "endianness" unless they use something else, be it QuTiP, numpy, or sympy, which I believe gives the same result different from QuantumOptics.jl.

Let's look at the well known CNOT gate, where there could be no uncertainty about basis or "endianness" (It's just the computational basis that everyone uses).

b = NLevelBasis(2)

pauli_x = Operator(b, [0 1; 1 0])
CNOT = tensor(dm(nlevelstate(b, 2)), pauli_x) + tensor(dm(nlevelstate(b, 1)), one(b))
Operator(dim=4x4)
  basis: [NLevel(N=2) ⊗ NLevel(N=2)]
 1.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  1.0+0.0im
 0.0+0.0im  0.0+0.0im  1.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im  0.0+0.0im  0.0+0.0im

This is clearly different from the standard CNOT definition in the computational basis of

$$ \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{pmatrix} $$

[amilsted]

One reason numpy, QuTiP, etc. have the "conventional convention" for tensor "endianness" (sometimes I call this "tensor layout") is that multi-dimensional arrays in numpy (and C) follow this convention in memory - if we walk through the block of memory where the array is stored, the rightmost index varies fastest. This is called "row-major" ordering or "C" ordering.

Julia and Fortran use the opposite "column major" (or Fortran) convention for multi-dimensional arrays - the leftmost index varies fastest.

Multi-dimensional arrays come into this because they are a very convenient way of addressing a tensor-product space. We can just reshape e.g. an operator matrix for multiple subsystems into a multidim array so that indices correspond to subsystems. Then we can use operations like numpy.transpose or Julia's permutedims to rearrange subsystems, or sum over certain axes to take partial traces, and so on.

The sticking point is that we want the reshape to be fast. Reshaping according to the tensor layout convention of the language or library is "free" - it's just a metadata change. Using a mismatched convention for language or library requires either reimplementing array indexing and reshaping or putting up with overhead from shuffling data around in memory to match the desired convention.

I agree it's unfortunate that there is this difference between numpy and Julia, and that the Julia/Fortran array layout convention doesn't match the (a?) textbook definition (or indeed the Julia definition!!) of the Kronecker product. It's however not a trivial matter to switch QO.jl to the numpy convention.

I realize I am not helping to propose a solution here. Just hoping to make the situation more understandable. I agree we should prominently highlight this difference in the documentation somewhere! We could even provide a wrapper around permutedims that puts operators or states into "C"/"row-major" order for easier cross-checking?

Regarding kron in Julia: It indeed uses the C convention and does not match QO.jl, which internally just defines tensor(A, B) = kron(B, A) in order to match the array storage convention better.

[Krastanov]

I mostly agree on this topic, just wanted to provide a minor clarification:

Let's look at the well known CNOT gate, where there could be no uncertainty about basis or "endianness" (It's just the computational basis that everyone uses).

There is indeed a lot of uncertainty about basis and endiandness even in this case. Some people use the basis 00, 01, 10, 11 (i.e. little index goes up first) and some people use the basis 00, 10, 01, 11 (i.e. big index goes up first). Here is for instance a similar issue coming up in qiskit https://quantumcomputing.stackexchange.com/questions/8244/big-endian-vs-little-endian-in-qiskit

[matchatealeaf]

Thanks for all the clarification and for pointing out the existence of the uncertainty even with the CNOT gate. I understand more about the circumstances now.

For the case of tensor(A, B) = kron(B, A), one can set QuantumOpticsBase.set_printing(standard_order=true) to see the standard ordering or row-major ordering.

Hence I wonder if it is possible to have something similar for Operator definition, e.g., Operator(tensor(b, b), M; standard_order=true) so that the user can provide a matrix M that was constructed or obtained in standard ordering.

I'm a fairly new user and I defined the CNOT, CZ, Toffoli gates etc. simply by doing Operator(tensor(b, b), M) with their corresponding matrices obtained from Wikipedia, none the wiser that they will not actually correspond to those gates in QO.jl.

[Krastanov]

Yes, such a helper function would certainly be a nice addition. Happy to help with reviewing a PR about this.

Going on a long tangent: We have two convenient ways to create these gates without copying over matrices.

Using the small (incomplete) symbolics library we have (it can be used to express symbolic expressions in different formalisms, like state vectors or stabilizer tableaux):

julia> using QuantumSymbolics, QuantumClifford, QuantumOptics

julia> express(CNOT)
Operator(dim=4x4)
  basis: [Spin(1/2) ⊗ Spin(1/2)]
 1.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  1.0+0.0im
 0.0+0.0im  0.0+0.0im  1.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im  0.0+0.0im  0.0+0.0im

julia> express(X1⊗X2)
Ket(dim=4)
  basis: [Spin(1/2) ⊗ Spin(1/2)]
  0.49999999999999994 + 0.0im
  0.49999999999999994 + 0.0im
 -0.49999999999999994 + 0.0im
 -0.49999999999999994 + 0.0im

julia> express(X1⊗X2, CliffordRepr())
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z_
+ _Z
𝒮𝓉𝒶𝒷
+ X_
- _X

Working strictly in QuantumOptics expressing CNOT and family in terms of sums of tensor products of projectors (making the meaning of the matrices more explicit).

julia> b = SpinBasis(1//2); l0 = spinup(b); l1 = spindown(b);

julia> projector(l0)⊗identityoperator(b) + projector(l1)⊗sigmax(b)
Operator(dim=4x4)
  basis: [Spin(1/2) ⊗ Spin(1/2)]
 1.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  1.0+0.0im
 0.0+0.0im  0.0+0.0im  1.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im  0.0+0.0im  0.0+0.0im