transpiler bug encountered

[yulunwang]

Information

• Qiskit Terra version: 0.18.1
• Python version: 3.7.6
• Operating system: Linux

What is the current behavior?

Got following errors while using crx and rx gates:

/gpfs/software/Anaconda/envs/psifour/lib/python3.7/site-packages/ipykernel_launcher.py:72: RuntimeWarning: invalid value encountered in double_scalars
/gpfs/home/yuluwang/.local/lib/python3.7/site-packages/qiskit/transpiler/runningpassmanager.py:166: UserWarning: Resynthesized [<qiskit.dagcircuit.dagnode.DAGNode object at 0x2aaba5022b40>, <qiskit.dagcircuit.dagnode.DAGNode object at 0x2aaba5022280>, <qiskit.dagcircuit.dagnode.DAGNode object at 0x2aaba5022910>, <qiskit.dagcircuit.dagnode.DAGNode object at 0x2aaba72f9d70>] and got global phase: π
┌───────┐┌────┐┌─────────────┐┌────┐┌──────────┐
qr_0: ┤ Rz(π) ├┤ √X ├┤ Rz(-3.1415) ├┤ √X ├┤ Rz(-π/2) ├
└───────┘└────┘└─────────────┘└────┘└──────────┘, but the original was native and the new value is longer. This indicates an efficiency bug in synthesis. Please report it by opening an issue here: https://github.com/Qiskit/qiskit-terra/issues/new/choose
new_dag = pass_.run(dag)

Steps to reproduce the problem

Using Rx and Crx gates. I applied parameter binding and noise model at the same time.

transpiled_circ = transpile(circ, backend=backend,optimization_level=3)
circ = transpiled_circ.bind_parameters(para_dict)
job = execute(circ, backend=backend, shots=shots,coupling_map=coupling_map, basis_gates=basis_gates, noise_model=noise_model)

What is the expected behavior?

Suggested solutions

[1ucian0]

It seems duplicated with #7033.

[1ucian0]

You mean the warnings? or is there an error?

[yulunwang]

@1ucian0 Hi Luciano, sorry it's a warning. It keeps showing up and asking me to report this issue.

[ecpeterson]

I don't think I can reproduce this without an example value of circ, and I mis-wrote the warning — the relevant information should be printed where all those DAGNodeOp strings are appearing instead. #7048 contains a better warning string. Could I ask you to supply an example circ or, if that's sensitive information, to re-run with the warning string from #7048 and paste its output?

[yulunwang]

It's hard to reproduce it because I think it connects to the value I assigned to theta. Once I changed the value the warning disappears. It's hard to reproduce in my algorithm where thousands circuits have been submitted and ran with such warnings appearing for several times. But fortunately in the end I reproduced it.

Here is one example for you:

backend = Aer.get_backend('qasm_simulator')

real_backend = provider.get_backend("ibmq_santiago")      
noise_model = NoiseModel.from_backend(real_backend)
# Get coupling map from backend
coupling_map = backend.configuration().coupling_map
# Get basis gates from noise model
basis_gates = noise_model.basis_gates

from qiskit.circuit import Parameter
theta0 = Parameter("θ_0")

cc=QuantumCircuit(2,1)
cc.h(0)
cc.crx(theta0,0,1)
cc.cx(0,1)
cc.cx(0,1)
cc.crx(-theta0,0,1)
cc.h(0)
cc.measure(0,0)
transpiled_circ = transpile(cc, backend=backend,optimization_level=3) 

para_dict = {theta0: 9.723602566887436e-05}
shots = 8192

qc = transpiled_circ.bind_parameters(para_dict)
job = execute(qc,backend=backend, shots=shots,coupling_map=coupling_map, basis_gates=basis_gates,  noise_model=noise_model)
counts = job.result().get_counts()

The warning I got:
/gpfs/home/yuluwang/.local/lib/python3.7/site-packages/qiskit/transpiler/runningpassmanager.py:166: UserWarning: Resynthesized [<qiskit.dagcircuit.dagnode.DAGNode object at 0x2aaba73e0f30>, <qiskit.dagcircuit.dagnode.DAGNode object at 0x2aaba73e0ec0>, <qiskit.dagcircuit.dagnode.DAGNode object at 0x2aaba73e0c90>, <qiskit.dagcircuit.dagnode.DAGNode object at 0x2aaba73e0e50>] and got ┌────────────────┐┌────┐┌─────────────┐┌────┐┌───────┐ qr_0: ┤ Rz(2.2835e-12) ├┤ √X ├┤ Rz(-3.1415) ├┤ √X ├┤ Rz(π) ├ └────────────────┘└────┘└─────────────┘└────┘└───────┘, but the original was native and the new value is longer. This indicates an efficiency bug in synthesis. Please report it by opening an issue here: https://github.com/Qiskit/qiskit-terra/issues/new/choose new_dag = pass_.run(dag) /gpfs/home/yuluwang/.local/lib/python3.7/site-packages/qiskit/transpiler/runningpassmanager.py:166: UserWarning: Resynthesized [<qiskit.dagcircuit.dagnode.DAGNode object at 0x2aaba73e6670>, <qiskit.dagcircuit.dagnode.DAGNode object at 0x2aaba73e0520>, <qiskit.dagcircuit.dagnode.DAGNode object at 0x2aaba73e0130>, <qiskit.dagcircuit.dagnode.DAGNode object at 0x2aaba73e0d70>] and got ┌────────────────┐┌────┐┌─────────────┐┌────┐┌───────┐ qr_0: ┤ Rz(2.2835e-12) ├┤ √X ├┤ Rz(-3.1415) ├┤ √X ├┤ Rz(π) ├ └────────────────┘└────┘└─────────────┘└────┘└───────┘, but the original was native and the new value is longer. This indicates an efficiency bug in synthesis. Please report it by opening an issue here: https://github.com/Qiskit/qiskit-terra/issues/new/choose new_dag = pass_.run(dag)

[ecpeterson]

Thanks so much! I can reproduce that.

[ecpeterson]

@levbishop Since you have strong opinions about atol, I think you should weigh in. The following circuit

qreg q[1];

rz(-pi) q;
sx q;
rz(pi - 5.5e-05) q;
sx q;

re-synthesizes to have ZXZ angle values (-1.0091434386183797e-12, 5.5000000000187176e-05, 1.0093214144495466e-12), which the PSX synthesizer shifts to (3.141592653588784, 3.1416476535897933, 1.0093214144495466e-12). Because 1.009e-12 is slightly above DEFAULT_ATOL, none of those gates are elided, which raises the warning.

This is sensitive to the complex difference in the definition of sx: replacing the first gate by rx(pi/2), which has purely real / purely imaginary entries yet is equivalent up to phase, causes the problem to go away. I propose that subtraction as the ultimate culprit.

I see two options:

• Raise DEFAULT_ATOL. I'm against doing that without a principled idea of how much of a raise is enough, though I suspect a wholly principled approach is not available, since it's surely influenced by LAPACK behavior.
• Start to measure circuit effectiveness by estimated infidelity rather than depth / gate count. Quite a lot more work, but probably the only way forward. Also, if it's possible to make this numerical effect appear on the "theta" coordinate, then this would not resolve the issue.

[levbishop]

I understand the issue and like you I also lean away from simply raising the tolerance without thinking through carefully.
Let me ponder this for a little.

[ecpeterson]

It's not yet clear if it's a related issue, but the first line in the log from #7041 shows a single gate (supposedly belonging to ['p', 'sx', 'x']) being mis-synthesized into a pair of sxes. That could well turn out to be a symptom of the theta error I was worried about above.

[levbishop]

I had a start at doing this based on estimated infidelity. Here's where I got to before I had to switch to other things.

    @staticmethod
    def _circuit_psx_gen_fid(tgt, theta, phi, lam, phase, infid, pfun, xfun, xpifun=None):
        """
        Generic X90, phase decomposition
        """
        qr = QuantumRegister(1, "qr")
        circuit = QuantumCircuit(qr, global_phase=phase)
        mintr2 = 4 - 6 * infid  # Minimum acceptable |Tr(Ud.Uapprox)|^2

        # Shift theta and phi to turn the decomposition from
        # RZ(phi).RY(theta).RZ(lam) = RZ(phi).RX(-pi/2).RZ(theta).RX(pi/2).RZ(lam)
        # into RZ(phi+pi).SX.RZ(theta+pi).SX.RZ(lam) .
        theta, phi = theta + np.pi, phi + np.pi
        phase -= np.pi / 2

        tr2_none = 4 * math.sin(theta / 2) ** 2 * math.sin((lam + phi) / 2) ** 2
        if tr2_none > mintr2:
            return circuit

        tr2_p = 2 - 2 * math.cos(theta)
        if tr2_p > mintr2:
            pfun(circuit, qr, lam + phi + np.pi)
            return circuit

        tr2_x = (
            1
            + math.cos(theta) * math.cos(phi) * math.cos(lam)
            + math.sin(phi) * math.sin(lam)
            - math.sin(theta) * (math.sin(phi) + math.sin(lam))
        )
        if tr2_x > mintr2:
            xfun(circuit, qr)
            return circuit

        angle_px = (
            math.atan2(math.sin(theta) - math.sin(phi), math.cos(theta) * math.cos(phi)) + lam
        )
        tr2_px = (
            math.cos((-theta - phi + lam - angle_px) / 2)
            + math.cos((theta - phi + lam - angle_px) / 2)
            - 2 * math.sin(theta / 2) * math.sin((phi + lam - angle_px) / 2)
        ) ** 2 / 2
        angle_xp = (
            math.atan2(math.sin(theta) - math.sin(lam), math.cos(theta) * math.cos(lam)) + phi
        )
        tr2_xp = (
            math.cos((-theta - phi + lam + angle_xp) / 2)
            + math.cos((theta - phi + lam + angle_xp) / 2)
            - 2 * math.sin(theta / 2) * math.sin((phi + lam - angle_xp) / 2)
        ) ** 2 / 2
        if max(tr2_px, tr2_xp) > mintr2:
            if tr2_px > tr2_xp:
                pfun(circuit, qr, angle_px)
                xfun(circuit, qr)
            else:
                xfun(circuit, qr)
                pfun(circuit, qr, angle_xp)
            return circuit

        tr2_pxp = 2 * (1 + abs(math.sin(theta)))
        if tr2_pxp > mintr2:
            offset = np.pi / 2 if theta < 0 else -np.pi / 2
            pfun(circuit, qr, lam + offset)
            xfun(circuit, qr)
            pfun(circuit, qr, phi + offset)
            return circuit

        tr2_xx = 4 * math.cos(theta / 2) ** 2 * math.cos((lam - phi) / 2) ** 2
        if tr2_xx > mintr2:
            if xpifun:
                xpifun(circuit, qr)
            else:
                xfun(circuit, qr)
                xfun(circuit, qr)
            return circuit

        tr2_xxp = 2 + 2 * math.cos(theta)
        if tr2_xxp > mintr2:
            if xpifun:
                xpifun(circuit, qr)
            else:
                xfun(circuit, qr)
                xfun(circuit, qr)
            pfun(circuit, qr, -lam + phi)
            return circuit

        angle_xpx = math.atan2(
            math.sin(theta) * (math.cos(phi) + math.cos(lam)),
            math.cos(theta) * (1 + math.cos(phi) * math.cos(lam)) + math.sin(phi) * math.sin(lam),
        )
        tr2_xpx = (
            4
            * (
                math.cos(theta / 2) * math.cos(angle_xpx / 2) * math.cos((phi - lam) / 2)
                + math.cos((phi + lam) / 2) * math.sin(theta / 2) * math.sin(angle_xpx / 2)
            )
            ** 2
        )
        if tr2_xpx > mintr2:
            xfun(circuit, qr)
            pfun(circuit, qr, angle_xpx)
            xfun(circuit, qr)
            return circuit

        # FIXME: figure out angles and traces for this case
        angle1_xpxp = 0
        angle2_xpxp = 0
        tr2_xpxp = 0
        angle1_pxpx = 0
        angle2_pxpx = 0
        tr2_pxpx = 0
        if max(tr2_pxpx, tr2_xpxp) > mintr2:
            if tr2_pxpx > tr2_xpxp:
                pfun(circuit, qr, angle1_pxpx)
                xfun(circuit, qr)
                pfun(circuit, qr, angle2_pxpx)
                xfun(circuit, qr)
            else:
                xfun(circuit, qr)
                pfun(circuit, qr, angle2_xpxp)
                xfun(circuit, qr)
                pfun(circuit, qr, angle1_xpxp)
            return circuit

        # No special case applies, use the generic pxpxp decomposition
        pfun(circuit, qr, lam)
        xfun(circuit, qr)
        pfun(circuit, qr, theta)
        xfun(circuit, qr)
        pfun(circuit, qr, phi)
        return circuit

It gives optimal approximations for nearly all special cases.
Open issues:

1. I didn't make any effort to get the global phase right for the approximations. That should be doable but needs care.
2. I didn't get optimal approximations for the xpxp and pxpx cases, and I'm guessing we'd need to do some power series or other approximation to the best choice for the angles. Even if we can't do much it should be better than the status quo of choosing angle1=theta, angle2=phi.
3. I didn't do profiling to see if all the trig calculations end up adding a meaningful performance hit (we could plausibly use cython etc to claw much of that back if needed)

Anyone want to pick this up? otherwise I'll try to find time but not likely in time for 0.19

[ecpeterson]

I'm happy to pick this up. I think coupling fidelity estimation and approximation dodges my worry about theta in option #2.

Similarly no guarantees about 0.19, but it could conceivably happen.