Represent a noisy gate directly, instead of separately composing noise with an ideal gate

[rjurga]

What is the expected behavior?

I am using qiskit-aer's noise to simulate noisy circuits, but the results are inconsistent with how one would do it analytically.

The two goals of this issue are:

• Understand the reasononing behing Qiskit's current method of applying noise.
• Discuss the possibility of a more accurate method.

Short summary

Qiskit applies noise by first performing the ideal gate, and then applying noise separately. However, when integrating the Lindblad equation for a noisy gate, yielding a signle superoperator representing the noisy gate, the result is different from doing the ideal gate and the noise separately and then composing them.

Detailed example

Here is an example for the X gate combined with phase and amplitude damping error:

import numpy as np
from qiskit import QuantumCircuit
from qiskit.providers.aer import noise
from qiskit.providers.aer.utils import insert_noise
import qiskit.quantum_info as qi

circuit = QuantumCircuit(1)
circuit.x(0)

noise_model = noise.NoiseModel()

# Noise parameters
delta = 1e-3
kappa = 1e-3
gate_time = np.pi / 2.0
param_amplitude = 1 - np.exp(-kappa * gate_time)
param_phase = 1 - np.exp(-4 * delta * gate_time)

error = noise.phase_amplitude_damping_error(param_amplitude, param_phase, 1)

print("Noise superop:\n\n", qi.SuperOp(error.to_instruction()).data)

noise_model.add_all_qubit_quantum_error(error, ["x"])
noisy_circuit = insert_noise(circuit, noise_model, transpile=True)
superop = qi.SuperOp(noisy_circuit.to_instruction())

print("\nNoisy X gate superop:\n\n", superop.data)

This prints

Noise superop:

 [[0.99843044+0.j 0.        +0.j 0.        +0.j 0.        +0.j]
  [0.        +0.j 0.99607577+0.j 0.        +0.j 0.        +0.j]
  [0.        +0.j 0.        +0.j 0.99607577+0.j 0.        +0.j]
  [0.00156956+0.j 0.        +0.j 0.        +0.j 1.        +0.j]]

Noisy X gate superop:

 [[0.        +0.j 0.        +0.j 0.        +0.j 0.99843044+0.j]
  [0.        +0.j 0.        +0.j 0.99607577+0.j 0.        +0.j]
  [0.        +0.j 0.99607577+0.j 0.        +0.j 0.        +0.j]
  [1.        +0.j 0.        +0.j 0.        +0.j 0.00156956+0.j]]

To compare with the theory, we integrate the Lindblad equation in Mathematica.

The Lindblad equation for a Hamiltonian H = g * sigma_x is:

d_rho / d_t = -i*g*[sigma_x, rho] + kappa*L[sigma_minus](rho) + delta*L[sigma_z](rho)

The brackets [sigma_x, rho] represent a commutator. Here the Lindbladian is defined as

L[A](rho) = A*rho*A_dag - 0.5*A_dag*A*rho - 0.5*rho*A_dag*A

where A_dag denotes the adjoint of A.

To get the Lindblad equation for the noise only, we set g to zero, so that the Hamiltonian corresponds to an identity gate. Then, if we integrate and reshape the result as a 4x4 tensor to get a superop, it agrees exactly with the noise superop obtained in Qiskit and printed above.

Furthermore, if we compose that noise superop with the superop corresponding to the ideal X gate by multiplying them together, such that the order corresponds to applying the ideal X gate first and then the noise, the result agrees with the noisy X gate superop obtained in Qiskit and printed above.

So far so good.

Now, if instead of doing the noise separately and then composing the result with the ideal gate, we set g to one, and integrate the Lindblad equation to get the noisy X gate superop directly, we obtain the following:

[[0.00137194, 0. + 1.10143*10^-7 I, 0. - 1.10143*10^-7 I, 0.998627],
 [0. - 0.000499203 I, -0.000587086, 0.996668, 0. - 0.000499424 I],
 [0. + 0.000499203 I, 0.996668, -0.000587086, 0. + 0.000499424 I],
 [0.998628, 0. - 1.10143*10^-7 I, 0. + 1.10143*10^-7 I, 0.00137318]]

This is different from composing the two superops together, and some of the terms have errors too large (1e-3) to be of numerical origin only.

Conclusion

Since it gives different results, I would like to understand whether there is a theoretical justification for applying the ideal gate first, and then the noise, separately. Also, it can be checked that the two separate superoperators do not commute, so is there also a justification for applying the ideal gate first and the noise second, instead of the other way around? Are there arguments to show that this approach is a reasonable approximation to the exact superop?

And finally, I would like to ask whether the Qiskit developers might consider adding the exact way of computing noise in Qiskit directly?

The way I imagine it:

• Compared to the phase_amplitude_damping_error function, the input parameters would need to be slightly different, i.e. kappa and delta from the example above instead of the param_amplitude and param_damping that the function currently expects, and also the gate time to perform the integration over.
• These parameters are fed to the Lindblad equation for whatever gate is used.
• The Lindblad equation is numerically integrated.
• This yields the correct superoperator for the noisy gate.

Thank you, I hope that I explained my questions clearly, but please ask if you find anything unclear.

[hhorii]

Thank you for your suggestion. Probably, your proposal will be related to qiskit-dynamics that provides Lindblad model. Maybe, if we will enhancement, the model can be directly specified to a noise model of Aer. I believe that @chriseclectic has opinions about integration of Aer and dynamics. @itoko will be help when noise transpiler pass needs to be enhanced (to replace a gate operation with a noisy operation instead of to append a noise operation after a gate operation).

[rjurga]

I wasn't aware of qiskit-dynamics, thanks for bringing this to my attention.

Otherwise, you summed up my suggestion pretty well, so I look forward to hearing about any related opinions and plans for the future.

[chriseclectic]

@rjurga Qiskit Aer is a circuit based simulator and is designed for Markovian circuit based noise models which can always be represented as an ideal operation followed and CPTP maps. You are free to insert CPTP operators anywhere in a circuit (eg via a transpilation pass) which means this is a very general and can capture basically any circuit based markovian noise. The NoiseModel class itself is a subset of this which is intended for non-parameterized local noise, which is to say the error terms following an ideal instruction apply only to the qubits involved in that instruction, and the noise does not depend on the parameters of the instructions, only the type and qubits. Similarly for the LocalNoisePass which is a more general transpiler pass version of this for local noise which can be parameterized so that the error terms depend on gate parameters.

This is a very general noise model, you just need to construct the error terms to insert yourself based on your assumptions about the physics of the device. If you want to construct errors as a noisy gate (ie via master equation simulation) this just means you need to extract the error term which can be done for example as noisey_gate.dot(gate().inverse()) for a superop/quantum error noise_gate and ideal Gate gate (you could also use the unitary operator for gate with unitary.adjoint()). Aer includes some very basic noise models and functions for common text-book noise processes but for more accurate models you should construct your own noise model.

If you want to build a noise model from a Lindblad equation you are assuming something about the physics of the system and the control you are doing to implement gates, which is great, but after doing your simulations you just need to convert it to the gate based model as either one of the built in transpiler passes or NoiseModel (if compatible) or your own transpiler pass. For a simplified lindblad noise model for a gate that you can solve analytically to obtain the noisy superoperator, such as in in your example, either building a custom NoiseModel or LocalNoisePass should work. For more complicated dynamics you should have a look at qiskit-dynamics as @hhorii mentioned, which is a time-dependent master-equation simulator. With it you can construct a system model as a schrodinger or lindblad equation and doing time-dependent master equation simulations for gates. We are still working on making it more integrated with Qiskit (it does not have a qiskit backend interface yet though there should be one in the coming months). We also have longer term plans to use this to generate circuit based noise models for Aer based on simulations of gate calibrations using the dynamics simulator.

Now with all that context explained, as for your original question about more accurate noise models, for the specific case of thermal relaxation something I have been thinking could be a nice improvement to the RelaxationNoisePass is to do something like you suggest to approximate the relaxation happening during the gate. Without assuming anything about the physics/dynamics of the device you can still do a simplified coarse approximation by treating the ideal unitary as being generated by a constant Hamiltonian H = log(U)/(-it). Under this approximation you can then define the noisy superop for the gate as S_err = S_gate.dot(U.adjoint()) where S_gate = exp(t_gate * (L + D)) for the T1/T2 dissipator D (which is used to define the thermal_relaxation_error when L=0), and unitary generator L rho = -i [H, rho]. This approx is nice in that it doesn't require any additional information other than the pass already contains, which is the ideal gates and their lengths.

[rjurga]

Thank you @chriseclectic. Your advice for how to achieve this in qiskit is very appreciated and it is the route I'm going for now.

Even though for simple cases the Lindbald model can be solved analytically, I still wanted to take a look at qiskit-dynamics. This is what I do to get the SuperOps for a noisy X gate, the same case as in my previous example:

import numpy as np
from qiskit.quantum_info import Operator, Pauli, SuperOp
from qiskit_dynamics import Solver

kappa = 1e-3
delta = 1e-3

gate_time = np.pi / 2.0

zeros = Operator(np.zeros((2, 2)))

sigma_x = zeros + Pauli('X')
sigma_y = zeros + Pauli('Y')
sigma_z = zeros + Pauli('Z')
sigma_min = 0.5 * (sigma_x - 1j * sigma_y)

H = sigma_x

L_ops = []
L_ops.append(np.sqrt(kappa) * sigma_min)
L_ops.append(np.sqrt(delta) * sigma_z)

solver = Solver(static_hamiltonian=H, static_dissipators=L_ops)
solver.model.evaluation_mode = "dense_vectorized"

y0 = SuperOp(np.identity(4))
sol = solver.solve(t_span=[0., gate_time], y0=y0, t_eval=[gate_time])

print(sol.y[-1].data)

This prints

[[ 1.37348582e-03+0.00000000e+00j  0.00000000e+00-2.50536776e-06j  0.00000000e+00+2.50536776e-06j  9.98625263e-01+0.00000000e+00j]
 [ 0.00000000e+00+4.96807264e-04j -5.85535804e-04+0.00000000e+00j  9.96666246e-01+0.00000000e+00j  0.00000000e+00+5.01817999e-04j]
 [ 0.00000000e+00-4.96807264e-04j  9.96666246e-01+0.00000000e+00j -5.85535804e-04+0.00000000e+00j  0.00000000e+00-5.01817999e-04j]
 [ 9.98626514e-01+0.00000000e+00j  0.00000000e+00+2.50536776e-06j  0.00000000e+00-2.50536776e-06j  1.37473660e-03+0.00000000e+00j]]

which is close enough to the previous result. It seems like a bit of a hack to use the current API of qiskit-dynamics like that, having to pass an identity SuperOp and set the evaluation mode to vectorized by hand, so I'm not sure if this is really the intended usage. On the other hand I'm definitely aware that it's in early development so no complaints here, it works.

Is supporting such use cases in the future in the scope of qiskit-dynamics goals, with simpler API? Something like a lindblad model method compute_superop that only needs t_eval. Or are such things too specific?

Thanks again!

[hhorii]

Is supporting such use cases in the future in the scope of qiskit-dynamics goals, with simpler API? Something like a lindblad model method compute_superop that only needs t_eval. Or are such things too specific?

I hope @chriseclectic answers to this question, but mainly it is more productive to discuss about this in https://github.com/Qiskit/qiskit-dynamics.