trace_out_qubits().to_statevector() buggy when used with gates of short duration

Hi everyone, I'm trying to run a noise pass on a circuit with 20 cv_delay gates and 20 save_statevector() functions. The goal is to save a snapshot of the state as it transforms under consecutive noise passes. This circuit is run without discretize on.

The issue is that trace_out_qubits().to_statevector() will start to fail towards the later cv_delay gates if the gate duration is short (10ns instead of 100ns). Below is a circuit demonstrating the issue. The circuit is only supposed to trace out a single ancilla qubit from the circuit statevector, and the statevector doesn't look too difficult. Yet it tends to fail.

I suspect this is because noise is applied approximately with the noise pass? Or would the team suspect something else?

import c2qa
import qiskit

qmr = c2qa.QumodeRegister(1, 3)
anc = qiskit.QuantumRegister(1)
circ = c2qa.CVCircuit(qmr, anc)
circ.initialize([1, 0], anc[0])
circ.cv_initialize([0, 0, 1], qmr[0])

gate_duration=10 #ns
cycles = 20

for i in range(cycles):
    circ.save_statevector(label='noise_cycle_{}'.format(i), pershot=True)
    circ.cv_delay(duration=gate_duration, qumode=qmr[0], unit="ns")

noise_pass = c2qa.kraus.PhotonLossNoisePass(photon_loss_rates=0.02, circuit=circ, time_unit="ns")
state, results= c2qa.util.simulate(circ, noise_passes=noise_pass, per_shot_state_vector=True)

sorted_lists = []
for i in range(cycles):
    snapshot_states = results.data()['noise_cycle_{}'.format(i)]
    sorted_lists.append(snapshot_states)

for index, cycle_list in enumerate(sorted_lists):
    for statevector in cycle_list:
        try:
            c2qa.util.trace_out_qubits(circ, statevector).to_statevector() ## This will fail at some point
        except Exception as e:
            print(statevector)
            print(index)
            raise e
---------------------------------------------------------------------------
Statevector([ 1.00000000e+00+0.j, -3.40935444e-16+0.j,  1.18762137e-32+0.j,
             -1.41712772e-48+0.j,  1.13917867e-48+0.j, -5.23638954e-49+0.j,
              3.32175841e-49+0.j,  2.50370198e-51+0.j,  0.00000000e+00+0.j,
              0.00000000e+00+0.j,  0.00000000e+00+0.j,  0.00000000e+00+0.j,
              0.00000000e+00+0.j,  0.00000000e+00+0.j,  0.00000000e+00+0.j,
              0.00000000e+00+0.j],
            dims=(2, 2, 2, 2))
14
---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
/Users/username/bosonic-qiskit/jupyter_clean.ipynb Cell 6 in 3
     31 print(statevector)
     32 print(index)
---> 33 raise e

/Users/username/bosonic-qiskit/jupyter_clean.ipynb Cell 6 in 2
     27 for statevector in cycle_list:
     28     try:
---> 29         c2qa.util.trace_out_qubits(circ, statevector).to_statevector()
     30     except Exception as e:
     31         print(statevector)

File /Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/qiskit/quantum_info/states/densitymatrix.py:809, in DensityMatrix.to_statevector(self, atol, rtol)
    806 if not is_hermitian_matrix(self._data, atol=atol, rtol=rtol):
    807     raise QiskitError("Not a valid density matrix (non-hermitian).")
--> 809 evals, evecs = np.linalg.eig(self._data)
    811 nonzero_evals = evals[abs(evals) > atol]
    812 if len(nonzero_evals) != 1 or not np.isclose(nonzero_evals[0], 1, atol=atol, rtol=rtol):

File <__array_function__ internals>:180, in eig(*args, **kwargs)

File /Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/numpy/linalg/linalg.py:1318, in eig(a)
   1315 extobj = get_linalg_error_extobj(
   1316     _raise_linalgerror_eigenvalues_nonconvergence)
   1317 signature = 'D->DD' if isComplexType(t) else 'd->DD'
-> 1318 w, vt = _umath_linalg.eig(a, signature=signature, extobj=extobj)
   1320 if not isComplexType(t) and all(w.imag == 0.0):
   1321     w = w.real

File /Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/numpy/linalg/linalg.py:95, in _raise_linalgerror_eigenvalues_nonconvergence(err, flag)
     94 def _raise_linalgerror_eigenvalues_nonconvergence(err, flag):
---> 95     raise LinAlgError("Eigenvalues did not converge")

LinAlgError: Eigenvalues did not converge

C2QA / bosonic-qiskit

trace_out_qubits().to_statevector() buggy when used with gates of short duration #93