`kraus_op` on subsystems of qubits gives incorrect results

[mishmash]

What happened?

qu.kraus_op seems to generally give the wrong answer when where is used to apply the channel to a subsystem of qubits.

What did you expect to happen?

No response

Minimal Complete Verifiable Example

import quimb as qu

# number of qubits
n = 3

# random density matrix
qu.seed_rand(7)
rho = qu.rand_rho(2**n)

# define Kraus channel to be the identity operator acting on first 2 qubits
subsystem = [0, 1]
Ek = [qu.eye(2**len(subsystem))]

# apply channel using `where=subsystem`
sigma = qu.kraus_op(rho, Ek, dims=[2] * n, where=subsystem)

print(f"trace of rho = {rho.trace()}")  # == 1
print(f"trace of sigma = {sigma.trace()}")  # != 1
print(f"sigma - rho = {sigma - rho}")  # != 0

Relevant log output

trace of rho = (0.9999999999999998+1.5100928363641297e-18j)
trace of sigma = (0.964801387123163+3.469446951953614e-18j)
sigma - rho = [[ 0.039714+3.008306e-18j -0.001623-1.086642e-01j -0.006906+1.505899e-02j -0.006388-2.602586e-02j
  -0.043405+2.146961e-02j -0.023163+3.133044e-02j  0.228115+4.406043e-03j  0.035957-1.339283e-02j]
 [-0.001623+1.086642e-01j  0.098862-5.944195e-19j  0.047327+2.426025e-02j -0.015361+2.920513e-02j
  -0.004649+4.929301e-02j -0.004117+2.975497e-02j  0.070554-2.972611e-02j  0.286179-4.813220e-02j]
 [-0.006906-1.505899e-02j  0.047327-2.426025e-02j -0.13845 -1.501673e-18j -0.043828-2.857833e-02j
  -0.016956-4.629843e-02j  0.01105 +1.959346e-02j -0.021501+2.958124e-02j  0.010718+4.511040e-02j]
 [-0.006388+2.602586e-02j -0.015361-2.920513e-02j -0.043828+2.857833e-02j -0.09578 -2.862415e-19j
  -0.029524-2.789070e-02j  0.025478-3.832833e-03j  0.037903+7.530370e-03j  0.035932-3.496102e-02j]
 [-0.043405-2.146961e-02j -0.004649-4.929301e-02j -0.016956+4.629843e-02j -0.029524+2.789070e-02j
  -0.072015-6.693820e-19j -0.004328-5.169321e-02j  0.015335+7.674499e-03j  0.005811-2.880603e-02j]
 [-0.023163-3.133044e-02j -0.004117-2.975497e-02j  0.01105 -1.959346e-02j  0.025478+3.832833e-03j
  -0.004328+5.169321e-02j -0.147569+2.680411e-19j  0.017468-3.835530e-02j -0.067012-1.973924e-02j]
 [ 0.228115-4.406043e-03j  0.070554+2.972611e-02j -0.021501-2.958124e-02j  0.037903-7.530370e-03j
   0.015335-7.674499e-03j  0.017468+3.835530e-02j  0.106042+1.734723e-18j  0.030934+4.450292e-02j]
 [ 0.035957+1.339283e-02j  0.286179+4.813220e-02j  0.010718-4.511040e-02j  0.035932+3.496102e-02j
   0.005811+2.880603e-02j -0.067012+1.973924e-02j  0.030934-4.450292e-02j  0.173997+0.000000e+00j]]

Anything else we need to know?

No response

Environment

>>> quimb.__version__
'1.6.1.dev42+g4215ae2'

[jcmgray]

Thanks for the issue @mishmash! I think the order of indices is incorrect for >=2-qubit kraus operators (with dimensions incorrectly interleaved), will fix shortly.

[jcmgray]

Should be fixed in in https://github.com/jcmgray/quimb/commit/f732b73122d43e925d5ef357a180e15598d65792.

[mishmash]

Excellent. Thanks @jcmgray!

[mishmash]

Looks like there's a lingering print statement. https://github.com/jcmgray/quimb/blob/3e34a2540074097d61c061c5beea8def7b58338d/quimb/calc.py#L200

[jcmgray]

oops thanks, fixed.