Qiskit / qiskit

Qiskit is an open-source SDK for working with quantum computers at the level of extended quantum circuits, operators, and primitives.
https://www.ibm.com/quantum/qiskit
Apache License 2.0
4.85k stars 2.29k forks source link

Deterministic different output distributions with semantically equivalent circuits #7770

Open MattePalte opened 2 years ago

MattePalte commented 2 years ago

Environment

What is happening?

Running two supposedly equivalent variants of the same circuit leads to surprisingly discordant results. The divergence is very clear, it is reproducible on my machine, but not yet on another laptop (with Ubuntu 20.04, python 3.8, same conda packages).

How can we reproduce the issue?

Unfortunately this issue seems very much tied with some specific characteristic of my machine, since it doesn't happen with my other laptop. But it is yet unclear to me what is the specific problem. Anyway, on my machine I run the following setup and it always produce the failure:

Program A:

from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit.circuit.library.standard_gates import *
from qiskit import Aer, transpile, execute
qr = QuantumRegister(5, name='qr')
cr = ClassicalRegister(5, name='cr')
qc = QuantumCircuit(qr, cr, name='qc')
qc.measure(qr, cr)
qc = transpile(qc, optimization_level=0, seed_transpiler=44)
counts_A = execute(qc, backend=Aer.get_backend('qasm_simulator'), shots=1024).result().get_counts(qc)
counts_A  # OUTPUT: {'00000': 1024}

Program B: (equivalent because I add a circuit and its inverse).

from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit.circuit.library.standard_gates import *
from qiskit import Aer, transpile, execute
qr = QuantumRegister(5, name='qr')
cr = ClassicalRegister(5, name='cr')
qc = QuantumCircuit(qr, cr, name='qc')
subcircuit = QuantumCircuit(qr, cr, name='subcircuit')
subcircuit.append(CSwapGate(), qargs=[qr[0], qr[2], qr[3]], cargs=[])
qc.append(subcircuit, qargs=qr, cargs=cr)
qc.append(subcircuit.inverse(), qargs=qr, cargs=cr)
qc.measure(qr, cr)
qc = transpile(qc, optimization_level=3, seed_transpiler=44)
qc.draw(fold=-1)

Output:

                                                      ┌───┐                                                   ┌───┐           ┌─┐      
qr_0: ─────────────────■─────────────────────■────■───┤ T ├───■────────────────■─────────────────────■────■───┤ T ├───■───────┤M├──────
                  ┌─┐  │                     │    │   └───┘   │                │                     │    │   └───┘   │       └╥┘      
qr_1: ────────────┤M├──┼─────────────────────┼────┼───────────┼────────────────┼─────────────────────┼────┼───────────┼────────╫───────
      ┌──────────┐└╥┘  │             ┌───┐   │  ┌─┴─┐┌─────┐┌─┴─┐┌──────────┐  │             ┌───┐   │  ┌─┴─┐┌─────┐┌─┴─┐┌───┐ ║ ┌─┐   
qr_2: ┤0         ├─╫───┼─────────■───┤ T ├───┼──┤ X ├┤ Tdg ├┤ X ├┤0         ├──┼─────────■───┤ T ├───┼──┤ X ├┤ Tdg ├┤ X ├┤ X ├─╫─┤M├───
      │  Unitary │ ║ ┌─┴─┐┌───┐┌─┴─┐┌┴───┴┐┌─┴─┐├───┤└┬───┬┘└───┘│  Unitary │┌─┴─┐┌───┐┌─┴─┐┌┴───┴┐┌─┴─┐├───┤└┬───┬┘└───┘└─┬─┘ ║ └╥┘┌─┐
qr_3: ┤1         ├─╫─┤ X ├┤ T ├┤ X ├┤ Tdg ├┤ X ├┤ T ├─┤ H ├──────┤1         ├┤ X ├┤ T ├┤ X ├┤ Tdg ├┤ X ├┤ T ├─┤ H ├────────■───╫──╫─┤M├
      └───┬─┬────┘ ║ └───┘└───┘└───┘└─────┘└───┘└───┘ └───┘      └──────────┘└───┘└───┘└───┘└─────┘└───┘└───┘ └───┘            ║  ║ └╥┘
qr_4: ────┤M├──────╫───────────────────────────────────────────────────────────────────────────────────────────────────────────╫──╫──╫─
          └╥┘      ║                                                                                                           ║  ║  ║ 
cr: 5/═════╩═══════╩═══════════════════════════════════════════════════════════════════════════════════════════════════════════╩══╩══╩═
           4       1                                                                                                           0  2  3 
qc.qasm(formatted=True)

Output:

OPENQASM 2.0;
include "qelib1.inc";
gate unitary140649878146400 p0,p1 {
    u3(pi/2,pi/2,-pi/2) p0;
    u3(pi/2,0,3*pi/4) p1;
    cx p0,p1;
    u3(pi/2,-pi/2,pi/2) p0;
    u3(pi/4,pi/4,-pi/2) p1;
}
gate unitary140650475314336 p0,p1 {
    u3(pi/2,0,3*pi/4) p0;
    cx p1,p0;
    u3(0,-0.14189705,0.14189705) p1;
    u3(pi/4,-3*pi/4,pi/2) p0;
}
qreg qr[5];
creg cr[5];
unitary140649878146400 qr[2],qr[3];
cx qr[0],qr[3];
t qr[3];
cx qr[2],qr[3];
t qr[2];
tdg qr[3];
cx qr[0],qr[3];
cx qr[0],qr[2];
t qr[0];
tdg qr[2];
cx qr[0],qr[2];
t qr[3];
h qr[3];
unitary140650475314336 qr[2],qr[3];
cx qr[0],qr[3];
t qr[3];
cx qr[2],qr[3];
t qr[2];
tdg qr[3];
cx qr[0],qr[3];
cx qr[0],qr[2];
t qr[0];
tdg qr[2];
cx qr[0],qr[2];
t qr[3];
h qr[3];
cx qr[3],qr[2];
measure qr[0] -> cr[0];
measure qr[1] -> cr[1];
measure qr[2] -> cr[2];
measure qr[3] -> cr[3];
measure qr[4] -> cr[4];

Execute the program B:

qc = QuantumCircuit.from_qasm_str(qc.qasm())
counts_B = execute(qc, backend=Aer.get_backend('qasm_simulator'), shots=1024).result().get_counts(qc)
counts_B  # OUTPUT: {'01000': 255, '00000': 269, '00100': 251, '01100': 249}

What should happen?

The second program B, should create the same output of the first one: # OUTPUT: {'00000': 1024}

Any suggestions?

I tried two other setups other than the failing one on the remote machine.

Hint: the error must be in the transpiler or qasm exporter, because when using the produced qasm file (visible above) it leads to failure (ALSO IN SETUP B and C)

qc = QuantumCircuit.from_qasm_str(QASM_COPIED)  # the one above in this issue
counts_B = execute(qc, backend=Aer.get_backend('qasm_simulator'), shots=1024).result().get_counts(qc)
counts_B   # OUTPUT: {'01100': 247, '00000': 280, '00100': 253, '01000': 244}

I attach here the list of dump created via conda env export > environment.yml: (Note that six pip packages were not installed in the second laptop, because they led to pip exceptions, they are marked with # LEADING TO PIP EXCEPTION in the list).

name: envname
channels:
  - conda-forge
  - anaconda
  - defaults
dependencies:
  - _libgcc_mutex=0.1=conda_forge
  - _openmp_mutex=4.5=1_gnu
  - argon2-cffi=21.1.0=py38h497a2fe_0
  - async_generator=1.10=py_0
  - backcall=0.2.0=pyh9f0ad1d_0
  - backports=1.0=py_2
  - backports.functools_lru_cache=1.6.4=pyhd8ed1ab_0
  - bleach=4.1.0=pyhd8ed1ab_0
  - ca-certificates=2021.10.8=ha878542_0
  - cffi=1.14.0=py38h2e261b9_0
  - colorama=0.4.4=pyh9f0ad1d_0
  - cycler=0.10.0=py_2
  - dbus=1.13.6=he372182_0
  - debugpy=1.4.1=py38h709712a_0
  - decorator=5.1.0=pyhd8ed1ab_0
  - defusedxml=0.7.1=pyhd8ed1ab_0
  - entrypoints=0.3=pyhd8ed1ab_1003
  - expat=2.4.1=h9c3ff4c_0
  - fontconfig=2.13.1=hba837de_1005
  - freetype=2.10.4=h0708190_1
  - gettext=0.19.8.1=hf34092f_1004
  - glib=2.58.3=py38h73cb85d_1004
  - gst-plugins-base=1.14.5=h0935bb2_2
  - gstreamer=1.14.5=h36ae1b5_2
  - icu=67.1=he1b5a44_0
  - importlib-metadata=4.8.1=py38h578d9bd_0
  - ipykernel=6.4.2=py38he5a9106_0
  - ipython=7.28.0=py38he5a9106_0
  - ipython_genutils=0.2.0=py_1
  - ipywidgets=7.6.5=pyhd8ed1ab_0
  - jbig=2.1=h7f98852_2003
  - jedi=0.18.0=py38h578d9bd_2
  - joblib=1.1.0=pyhd8ed1ab_0
  - jpeg=9d=h36c2ea0_0
  - jsonschema=4.1.2=pyhd8ed1ab_0
  - jupyter=1.0.0=py38h578d9bd_6
  - jupyter_console=6.4.0=pyhd8ed1ab_0
  - jupyter_contrib_core=0.3.3=py_2
  - jupyter_contrib_nbextensions=0.5.1=pyhd8ed1ab_2
  - jupyter_core=4.8.1=py38h578d9bd_0
  - jupyter_highlight_selected_word=0.2.0=py38h578d9bd_1002
  - jupyter_latex_envs=1.4.6=pyhd8ed1ab_1002
  - jupyter_nbextensions_configurator=0.4.1=py38h578d9bd_2
  - jupyterlab_pygments=0.1.2=pyh9f0ad1d_0
  - jupyterlab_widgets=1.0.2=pyhd8ed1ab_0
  - kiwisolver=1.3.2=py38h1fd1430_0
  - krb5=1.17.2=h926e7f8_0
  - lcms2=2.12=hddcbb42_0
  - ld_impl_linux-64=2.36.1=hea4e1c9_2
  - lerc=3.0=h9c3ff4c_0
  - libblas=3.9.0=12_linux64_openblas
  - libcblas=3.9.0=12_linux64_openblas
  - libclang=10.0.1=default_hde54327_1
  - libdeflate=1.8=h7f98852_0
  - libedit=3.1.20191231=he28a2e2_2
  - libevent=2.1.10=h9b69904_4
  - libffi=3.2.1=he1b5a44_1007
  - libgcc-ng=11.2.0=h1d223b6_11
  - libgfortran-ng=11.2.0=h69a702a_11
  - libgfortran5=11.2.0=h5c6108e_11
  - libgomp=11.2.0=h1d223b6_11
  - libiconv=1.16=h516909a_0
  - liblapack=3.9.0=12_linux64_openblas
  - libllvm10=10.0.1=he513fc3_3
  - libopenblas=0.3.18=pthreads_h8fe5266_0
  - libpng=1.6.37=h21135ba_2
  - libpq=12.3=h255efa7_3
  - libsodium=1.0.18=h36c2ea0_1
  - libstdcxx-ng=11.2.0=he4da1e4_11
  - libtiff=4.3.0=h6f004c6_2
  - libuuid=2.32.1=h7f98852_1000
  - libwebp-base=1.2.1=h7f98852_0
  - libxcb=1.13=h7f98852_1003
  - libxkbcommon=0.10.0=he1b5a44_0
  - libxml2=2.9.10=h68273f3_2
  - libxslt=1.1.33=hf705e74_1
  - libzlib=1.2.11=h36c2ea0_1013
  - lxml=4.6.3=py38hf1fe3a4_0
  - lz4-c=1.9.3=h9c3ff4c_1
  - markupsafe=2.0.1=py38h497a2fe_0
  - matplotlib=3.4.3=py38h578d9bd_1
  - matplotlib-base=3.4.3=py38hf4fb855_1
  - matplotlib-inline=0.1.3=pyhd8ed1ab_0
  - mistune=0.8.4=py38h497a2fe_1004
  - mysql-common=8.0.26=ha770c72_0
  - mysql-libs=8.0.26=hfa10184_0
  - nbclient=0.5.4=pyhd8ed1ab_0
  - nbconvert=5.6.1=pyhd8ed1ab_2
  - nbformat=5.1.3=pyhd8ed1ab_0
  - ncurses=6.2=h58526e2_4
  - nest-asyncio=1.5.1=pyhd8ed1ab_0
  - notebook=6.4.5=pyha770c72_0
  - nspr=4.30=h9c3ff4c_0
  - nss=3.69=hb5efdd6_1
  - olefile=0.46=pyh9f0ad1d_1
  - openjpeg=2.4.0=hb52868f_1
  - openssl=1.1.1l=h7f98852_0
  - packaging=21.0=pyhd8ed1ab_0
  - pandas=1.3.4=py38h43a58ef_0
  - pandoc=2.14.2=h7f98852_0
  - pandocfilters=1.5.0=pyhd8ed1ab_0
  - parso=0.8.2=pyhd8ed1ab_0
  - patsy=0.5.2=pyhd8ed1ab_0
  - pcre=8.45=h9c3ff4c_0
  - pexpect=4.8.0=pyh9f0ad1d_2
  - pickleshare=0.7.5=py_1003
  - pillow=8.3.2=py38h8e6f84c_0
  - pip=21.3=pyhd8ed1ab_0
  - prometheus_client=0.11.0=pyhd8ed1ab_0
  - prompt-toolkit=3.0.20=pyha770c72_0
  - prompt_toolkit=3.0.20=hd8ed1ab_0
  - pthread-stubs=0.4=h36c2ea0_1001
  - ptyprocess=0.7.0=pyhd3deb0d_0
  - pycparser=2.20=pyh9f0ad1d_2
  - pygments=2.10.0=pyhd8ed1ab_0
  - pyparsing=2.4.7=pyh9f0ad1d_0
  - pyqt=5.12.3=py38h578d9bd_7
  - pyqt-impl=5.12.3=py38h7400c14_7
  - pyqt5-sip=4.19.18=py38h709712a_7
  - pyqtchart=5.12=py38h7400c14_7
  - pyqtwebengine=5.12.1=py38h7400c14_7
  - pyrsistent=0.17.3=py38h497a2fe_2
  - python=3.8.0=h357f687_5
  - python-dateutil=2.8.2=pyhd8ed1ab_0
  - python_abi=3.8=2_cp38
  - pytz=2021.3=pyhd8ed1ab_0
  - pyyaml=6.0=py38h497a2fe_0
  - pyzmq=22.3.0=py38h2035c66_0
  - qt=5.12.9=h1f2b2cb_0
  - qtconsole=5.1.1=pyhd8ed1ab_0
  - qtpy=1.11.2=pyhd8ed1ab_0
  - readline=8.1=h46c0cb4_0
  - scikit-learn=1.0=py38hacb3eff_1
  - seaborn=0.11.2=hd8ed1ab_0
  - seaborn-base=0.11.2=pyhd8ed1ab_0
  - send2trash=1.8.0=pyhd8ed1ab_0
  - setuptools=58.2.0=py38h578d9bd_0
  - six=1.16.0=pyh6c4a22f_0
  - sqlite=3.36.0=h9cd32fc_2
  - statsmodels=0.13.0=py38h6c62de6_0
  - terminado=0.12.1=py38h578d9bd_0
  - testpath=0.5.0=pyhd8ed1ab_0
  - threadpoolctl=3.0.0=pyh8a188c0_0
  - tk=8.6.11=h27826a3_1
  - tornado=6.1=py38h497a2fe_1
  - tqdm=4.62.3=pyhd8ed1ab_0
  - traitlets=5.1.0=pyhd8ed1ab_0
  - wcwidth=0.2.5=pyh9f0ad1d_2
  - webencodings=0.5.1=py_1
  - wheel=0.37.0=pyhd8ed1ab_1
  - widgetsnbextension=3.5.1=py38h578d9bd_4
  - xorg-libxau=1.0.9=h7f98852_0
  - xorg-libxdmcp=1.1.3=h7f98852_0
  - xz=5.2.5=h516909a_1
  - yaml=0.2.5=h516909a_0
  - zeromq=4.3.4=h9c3ff4c_1
  - zipp=3.6.0=pyhd8ed1ab_0
  - zlib=1.2.11=h36c2ea0_1013
  - zstd=1.5.0=ha95c52a_0
  - pip:
    - antlerinator==1!2.0.1
    - antlr4-python3-runtime==4.9.2
    - astor==0.8.1
    - astpretty==2.1.0
    - astunparse==1.6.3
    - attrs==20.3.0
    - autopep8==1.5.7
    - blis==0.7.5
    - catalogue==2.0.6
    - certifi==2021.5.30
    - charset-normalizer==2.0.7
    - cirq==0.13.1
    - cirq-aqt==0.13.1
    - cirq-core==0.13.1
    - cirq-google==0.13.1
    - cirq-ionq==0.13.1
    - cirq-pasqal==0.13.1
    - cirq-rigetti==0.13.1
    - cirq-web==0.13.1
    - click==8.0.3
    - coverage==6.3.2
    - cryptography==35.0.0
    - cupy-cuda112==9.6.0
    - cupy-cuda91==7.8.0
    - cupy-cuda92==9.6.0
    - cymem==2.0.6
    - deap==1.3.1
    - deprecated==1.2.13
    - dill==0.3.4
    - dlx==1.0.4
    - docplex==2.22.213
    - duet==0.2.3
    - en-core-web-lg==3.2.0    # LEADING TO PIP EXCEPTION DURING INSTALLATION ON NEW LAPTOP
    - en-core-web-sm==3.2.0   # LEADING TO PIP EXCEPTION
    - en-core-web-trf==3.2.0   # LEADING TO PIP EXCEPTION
    - fastdtw==0.3.4
    - fastjsonschema==2.15.1
    - fastrlock==0.8
    - filelock==3.4.2
    - flask==2.0.2   # LEADING TO PIP EXCEPTION
    - grammarinator==19.3.post90+g583ebca   # LEADING TO PIP EXCEPTION
    - grpcio==1.41.0
    - h11==0.9.0
    - h5py==3.2.1
    - httpcore==0.11.1
    - httpx==0.15.5
    - huggingface-hub==0.4.0
    - idna==2.10
    - immutables==0.6
    - inators==2.1.0
    - inflection==0.5.1
    - iso8601==0.1.16
    - itsdangerous==2.0.1
    - jinja2==2.11.0
    - jupyter-client==6.1.12
    - langcodes==3.3.0
    - lark==0.11.3
    - lark-parser==0.12.0
    - llvmlite==0.37.0
    - more-itertools==8.10.0
    - msgpack==0.6.2
    - multitasking==0.0.9
    - murmurhash==1.0.6
    - nltk==3.6.5
    - ntlm-auth==1.5.0
    - numba==0.54.1
    - numpy==1.20.0
    - openfermion==1.3.0
    - pathy==0.6.1
    - pbr==5.8.0
    - peewee==3.14.8
    - ply==3.11
    - preshed==3.0.6
    - psutil==5.8.0
    - pycodestyle==2.8.0
    - pydantic==1.8.2
    - pydot==1.4.2
    - pyjwt==1.7.1
    - pylatexenc==2.10
    - pyquil==3.0.1
    - python-afl==0.7.3
    - python-constraint==1.4.0
    - python-graphviz==0.14
    - python-rapidjson==1.5
    - pytket==0.18.0
    - pytket-cirq==0.19.0
    - pytket-qiskit==0.21.0
    - qcs-api-client==0.8.0
    - qiskit==0.33.1
    - qiskit-aer==0.9.1
    - qiskit-aer-gpu==0.9.1
    - qiskit-aqua==0.9.5
    - qiskit-ibmq-provider==0.18.2
    - qiskit-ignis==0.7.0
    - qiskit-terra==0.19.1
    - quandl==3.6.1
    - regex==2021.11.10
    - requests==2.26.0
    - requests-ntlm==1.1.0
    - retry==0.9.2
    - retrying==1.3.3
    - retworkx==0.10.2
    - rfc3339==6.2
    - rfc3986==1.5.0
    - rpcq==3.9.2
    - ruamel-yaml==0.17.17
    - ruamel-yaml-clib==0.2.6
    - sacremoses==0.0.47
    - scipy==1.7.3
    - smart-open==5.2.1
    - snakeviz==2.1.1
    - sniffio==1.2.0
    - spacy==3.2.1
    - spacy-alignments==0.8.4
    - spacy-legacy==3.0.8
    - spacy-loggers==1.0.1
    - spacy-transformers==1.1.4
    - srsly==2.4.2
    - stevedore==3.5.0
    - symengine==0.8.1
    - termcolor==1.1.0
    - thinc==8.0.13
    - tokenizers==0.10.3
    - torch==1.10.0
    - torch-two-sample==0.1   # LEADING TO PIP EXCEPTION
    - transformers==4.15.0
    - tweedledum==1.1.1
    - typer==0.4.0
    - types-jinja2==2.11.9
    - types-markupsafe==1.1.10
    - types-pkg-resources==0.1.3
    - typing-extensions==3.10.0.2
    - urllib3==1.26.7
    - wasabi==0.9.0
    - websocket-client==1.2.1
    - werkzeug==2.0.2
    - wrapt==1.13.3
    - yfinance==0.1.64
    - z3-solver==4.8.14.0
prefix: /home/username/.conda/envs/envname

I am really looking forward to your feedback since I have no clue on why this happen and more importantly how to reproduce it on another machine, which is crucial for all of us to understand this bug.

  1. Which could be the source of divergence? Any idea of which python or linux library might be involved?
  2. Might the seed be responsible for this? Are they machine specific? Thanks in advance
MattePalte commented 2 years ago

The problem seems related to the transformation to qasm format (qc.qasm()), because the two machine use different decomposers TwoQubitBasisDecomposer, which in turn come from two different TwoQubitWeylDecomposition.

The TwoQubitBasisDecomposer is initialized here in my qiskit version: https://github.com/Qiskit/qiskit-terra/blob/a8a0d1be26b5df097e7a21016fa21550a3c4efce/qiskit/quantum_info/synthesis/two_qubit_decompose.py#L1405

But it is similar in the latest version: https://github.com/Qiskit/qiskit-terra/blob/fe8ab43f1218d6997458fb3bf4011045f40932b0/qiskit/quantum_info/synthesis/two_qubit_decompose.py#L1449 https://github.com/Qiskit/qiskit-terra/blob/fe8ab43f1218d6997458fb3bf4011045f40932b0/qiskit/quantum_info/synthesis/two_qubit_decompose.py#L1422

The difference can be exposed with this even minimal example.

from qiskit.quantum_info.synthesis.two_qubit_decompose import TwoQubitWeylDecomposition
from qiskit.quantum_info.operators import Operator
from qiskit.circuit.library.standard_gates import CXGate

gate = CXGate()
weyl = TwoQubitWeylDecomposition(Operator(gate).data)
weyl.__dict__

Output SETUP A:

{'_original_decomposition': TwoQubitWeylDecomposition.from_bytes(
     # TwoQubitWeylDecomposition(
     #  global phase: π/4
     #       ┌──────────┐┌──────────┐            ┌────────────┐┌────────┐┌──────────┐
     #  q_0: ┤ Rz(-π/2) ├┤ Ry(-π/2) ├────────────┤0           ├┤ Rz(-π) ├┤ Ry(-π/2) ├
     #       ├──────────┤├─────────┬┘┌──────────┐│  Rxx(-π/2) │├───────┬┘└──────────┘
     #  q_1: ┤ Rz(-π/2) ├┤ Ry(π/2) ├─┤ Rz(-π/2) ├┤1           ├┤ Ry(π) ├─────────────
     #       └──────────┘└─────────┘ └──────────┘└────────────┘└───────┘
     # )
     b'k05VTVBZAQB2AHsnZGVzY3InOiAnPGMxNicsICdmb3J0cmFuX29yZGVyJzogRmFsc2UsICdzaGFw'
     b'ZSc6ICg0LCA0KSwgfSAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg'
     b'ICAgICAgICAgICAgIAoAAAAAAADwPwAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA'
     b'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA'
     b'AAAAAAAAAAAAAAAAAAAAAAAA8D8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA'
     b'AAAAAAAAAAAA8D8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA'
     b'8D8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA',
     requested_fidelity=0.999999999,
     calculated_fidelity=1.0,
     actual_fidelity=0.9999999999999997,
     abc=(0.7853981633974483, 0.0, 0.0)),
 'a': 0.7853981633974483,
 'b': 0,
 'c': 0,
 'K1l': array([[-1.05226633e-16+3.82683432e-01j,  9.23879533e-01-5.36741351e-17j],
        [-9.23879533e-01+6.41070206e-18j,  3.98309960e-17-3.82683432e-01j]]),
 'K1r': array([[-0.5-0.5j,  0.5+0.5j],
        [-0.5+0.5j, -0.5+0.5j]]),
 'K2l': array([[ 5.65713056e-17+3.82683432e-01j, -9.23879533e-01-2.34326020e-17j],
        [ 9.23879533e-01-2.34326020e-17j,  5.65713056e-17-3.82683432e-01j]]),
 'K2r': array([[ 4.32978028e-17-0.70710678j, -4.32978028e-17+0.70710678j],
        [ 4.32978028e-17+0.70710678j,  4.32978028e-17+0.70710678j]]),
 'global_phase': -2.356194490192345,
 'unitary_matrix': array([[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
        [0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j],
        [0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j],
        [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j]]),
 'requested_fidelity': None,
 '_is_flipped_from_original': False,
 'calculated_fidelity': 1.0}

Output SETUP B:

{'_original_decomposition': TwoQubitWeylDecomposition.from_bytes(
     # TwoQubitWeylDecomposition(
     #  global phase: π/4
     #       ┌────────────┐   ┌─────────┐              ┌────────────┐┌──────────┐┌─────────────┐
     #  q_0: ┤ Rz(3.1389) ├───┤ Ry(π/2) ├──────────────┤0           ├┤ Ry(-π/2) ├┤ Rz(-1.5681) ├────────────
     #       └┬──────────┬┘┌──┴─────────┴──┐┌─────────┐│  Rxx(-π/2) │├─────────┬┘└┬────────────┤┌──────────┐
     #  q_1: ─┤ Rz(-π/2) ├─┤ Ry(0.0026466) ├┤ Rz(π/2) ├┤1           ├┤ Rz(π/2) ├──┤ Ry(1.5734) ├┤ Rz(-π/2) ├
     #        └──────────┘ └───────────────┘└─────────┘└────────────┘└─────────┘  └────────────┘└──────────┘
     # )
     b'k05VTVBZAQB2AHsnZGVzY3InOiAnPGMxNicsICdmb3J0cmFuX29yZGVyJzogRmFsc2UsICdzaGFw'
     b'ZSc6ICg0LCA0KSwgfSAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg'
     b'ICAgICAgICAgICAgIAoAAAAAAADwPwAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA'
     b'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA'
     b'AAAAAAAAAAAAAAAAAAAAAAAA8D8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA'
     b'AAAAAAAAAAAA8D8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA'
     b'8D8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA',
     requested_fidelity=0.999999999,
     calculated_fidelity=1.0,
     actual_fidelity=1.0,
     abc=(0.7853981633974483, 0.0, 0.0)),
 'a': 0.7853981633974483,
 'b': 0,
 'c': 0,
 'K1l': array([[ 9.99999781e-01-3.92003950e-17j, -7.84266235e-17-6.61646830e-04j],
        [-7.84785659e-17-6.61646830e-04j,  9.99999781e-01+3.93042108e-17j]]),
 'K1r': array([[ 0.5-0.5j,  0.5-0.5j],
        [-0.5-0.5j,  0.5+0.5j]]),
 'K2l': array([[ 7.07574481e-01+4.04613999e-20j, -4.05149776e-20-7.06638771e-01j],
        [ 4.05149776e-20-7.06638771e-01j,  7.07574481e-01-4.04613999e-20j]]),
 'K2r': array([[ 0.70710678+0.j, -0.70710678+0.j],
        [ 0.70710678+0.j,  0.70710678+0.j]]),
 'global_phase': 0.7853981633974487,
 'unitary_matrix': array([[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
        [0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j],
        [0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j],
        [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j]]),
 'requested_fidelity': None,
 '_is_flipped_from_original': False,
 'calculated_fidelity': 1.0}

This method contains a bit of randomness, which might be the source of the divergence, and can be also a good moment to review it (given the comment): https://github.com/Qiskit/qiskit-terra/blob/fe8ab43f1218d6997458fb3bf4011045f40932b0/qiskit/quantum_info/synthesis/two_qubit_decompose.py#L174-L175

Looking forward to know which result is the given code on your machine with the qiskit version used by me. Thanks in advance

MattePalte commented 2 years ago

To further check that the two circuits are far from each other, I got their unitary:

Setup A:

from math import pi
qr = QuantumRegister(2, name='qr')
cr = ClassicalRegister(2, name='cr')
qc = QuantumCircuit(qr, cr, name='qc')

# SECTION
# NAME: MEASUREMENT

qc.append(RZGate(-pi/2), qargs=[qr[0]], cargs=[])
qc.append(RYGate(-pi/2), qargs=[qr[0]], cargs=[])
qc.append(RZGate(-pi/2), qargs=[qr[1]], cargs=[])
qc.append(RYGate(pi/2), qargs=[qr[1]], cargs=[])
qc.append(RZGate(-pi/2), qargs=[qr[1]], cargs=[])

qc.append(RXXGate(-pi/2), qargs=[qr[0], qr[1]], cargs=[])

qc.append(RZGate(-pi), qargs=[qr[0]], cargs=[])
qc.append(RYGate(pi), qargs=[qr[1]], cargs=[])
qc.append(RYGate(-pi/2), qargs=[qr[0]], cargs=[])
qc.save_unitary()
#qc.measure(qr, cr)
qc.draw(fold=-1)
simulator = Aer.get_backend('aer_simulator')
qc_trans = transpile(qc, simulator)

# Run and get unitary
unitary = simulator.run(qc_trans).result().get_unitary(qc)
unitary

Output:

      ┌──────────┐┌──────────┐            ┌────────────┐┌────────┐┌──────────┐ ░ 
qr_0: ┤ Rz(-π/2) ├┤ Ry(-π/2) ├────────────┤0           ├┤ Rz(-π) ├┤ Ry(-π/2) ├─░─
      ├──────────┤├─────────┬┘┌──────────┐│  Rxx(-π/2) │├───────┬┘└──────────┘ ░ 
qr_1: ┤ Rz(-π/2) ├┤ Ry(π/2) ├─┤ Rz(-π/2) ├┤1           ├┤ Ry(π) ├──────────────░─
      └──────────┘└─────────┘ └──────────┘└────────────┘└───────┘              ░ 
cr: 2/═══════════════════════════════════════════════════════════════════════════
array([[ 7.07106781e-01-7.07106781e-01j,  9.75980418e-34+6.70078871e-17j,
         5.18672754e-17-2.35174251e-16j, -6.50353591e-17+7.85046229e-17j],
       [ 6.62912745e-17-1.62588398e-17j, -4.23815048e-17+7.20465508e-17j,
         4.73817313e-17-1.62464020e-32j,  7.07106781e-01-7.07106781e-01j],
       [-2.67754988e-16+1.31296464e-17j, -5.55111512e-17+5.55111512e-17j,
         7.07106781e-01-7.07106781e-01j, -1.02350467e-32+2.77555756e-17j],
       [ 2.77555756e-17+1.16264517e-33j,  7.07106781e-01-7.07106781e-01j,
         3.37735950e-17-6.50353591e-17j,  7.48991843e-17-6.53119189e-17j]])

Setup B:

from math import pi
qr = QuantumRegister(2, name='qr')
cr = ClassicalRegister(2, name='cr')
qc = QuantumCircuit(qr, cr, name='qc')

# SECTION
# NAME: MEASUREMENT

qc.append(RZGate(3.1389), qargs=[qr[0]], cargs=[])
qc.append(RYGate(pi/2), qargs=[qr[0]], cargs=[])
qc.append(RZGate(-pi/2), qargs=[qr[1]], cargs=[])
qc.append(RYGate(0.0026466), qargs=[qr[1]], cargs=[])
qc.append(RZGate(pi/2), qargs=[qr[1]], cargs=[])

qc.append(RXXGate(-pi/2), qargs=[qr[0], qr[1]], cargs=[])

qc.append(RYGate(-pi/2), qargs=[qr[0]], cargs=[])
qc.append(RZGate(-1.5681), qargs=[qr[0]], cargs=[])

qc.append(RZGate(pi/2), qargs=[qr[1]], cargs=[])
qc.append(RYGate(1.5734), qargs=[qr[1]], cargs=[])
qc.append(RZGate(-pi/2), qargs=[qr[1]], cargs=[])
qc.save_unitary()
#qc.measure(qr, cr)
qc.draw(fold=-1)
simulator = Aer.get_backend('aer_simulator')
qc_trans = transpile(qc, simulator)

# Run and get unitary
unitary = simulator.run(qc_trans).result().get_unitary(qc)
unitary

Output B:

      ┌────────────┐   ┌─────────┐              ┌────────────┐┌──────────┐┌─────────────┐             ░ 
qr_0: ┤ Rz(3.1389) ├───┤ Ry(π/2) ├──────────────┤0           ├┤ Ry(-π/2) ├┤ Rz(-1.5681) ├─────────────░─
      └┬──────────┬┘┌──┴─────────┴──┐┌─────────┐│  Rxx(-π/2) │├─────────┬┘└┬────────────┤┌──────────┐ ░ 
qr_1: ─┤ Rz(-π/2) ├─┤ Ry(0.0026466) ├┤ Rz(π/2) ├┤1           ├┤ Rz(π/2) ├──┤ Ry(1.5734) ├┤ Rz(-π/2) ├─░─
       └──────────┘ └───────────────┘└─────────┘└────────────┘└─────────┘  └────────────┘└──────────┘ ░ 
cr: 2/══════════════════════════════════════════════════════════════════════════════════════════════════
array([[ 7.07105482e-01-7.07108080e-01j, -6.65772178e-17+1.30397053e-17j,
         1.51769418e-05+1.51768860e-05j, -1.27068414e-16-1.80165683e-17j],
       [-2.97362235e-17-8.14854106e-17j,  1.51768860e-05+1.51769418e-05j,
         5.36602501e-17-9.77166068e-17j,  7.07108080e-01-7.07105482e-01j],
       [ 1.51769418e-05+1.51768860e-05j, -1.27170614e-16-1.80165683e-17j,
         7.07105482e-01-7.07108080e-01j, -6.65772178e-17+9.15442617e-17j],
       [ 5.36602501e-17-1.25433065e-16j,  7.07108080e-01-7.07105482e-01j,
        -1.94750686e-18-8.14854106e-17j,  1.51768860e-05+1.51769418e-05j]])

And by checking the difference, they are far as per numpy standard threshold:

import numpy as np
mat_a = np.array([[ 7.07106781e-01-7.07106781e-01j,  9.75980418e-34+6.70078871e-17j,
         5.18672754e-17-2.35174251e-16j, -6.50353591e-17+7.85046229e-17j],
       [ 6.62912745e-17-1.62588398e-17j, -4.23815048e-17+7.20465508e-17j,
         4.73817313e-17-1.62464020e-32j,  7.07106781e-01-7.07106781e-01j],
       [-2.67754988e-16+1.31296464e-17j, -5.55111512e-17+5.55111512e-17j,
         7.07106781e-01-7.07106781e-01j, -1.02350467e-32+2.77555756e-17j],
       [ 2.77555756e-17+1.16264517e-33j,  7.07106781e-01-7.07106781e-01j,
         3.37735950e-17-6.50353591e-17j,  7.48991843e-17-6.53119189e-17j]])
mat_b = np.array([[ 7.07105482e-01-7.07108080e-01j, -6.65772178e-17+1.30397053e-17j,
         1.51769418e-05+1.51768860e-05j, -1.27068414e-16-1.80165683e-17j],
       [-2.97362235e-17-8.14854106e-17j,  1.51768860e-05+1.51769418e-05j,
         5.36602501e-17-9.77166068e-17j,  7.07108080e-01-7.07105482e-01j],
       [ 1.51769418e-05+1.51768860e-05j, -1.27170614e-16-1.80165683e-17j,
         7.07105482e-01-7.07108080e-01j, -6.65772178e-17+9.15442617e-17j],
       [ 5.36602501e-17-1.25433065e-16j,  7.07108080e-01-7.07105482e-01j,
        -1.94750686e-18-8.14854106e-17j,  1.51768860e-05+1.51769418e-05j]])

np.allclose(mat_a, mat_b)
# Output: False
MattePalte commented 2 years ago

Looking for help, I tag some of those who introduced the FIXME according to git blame, feel free to involve the specific code owners if you are not responsible, Thanks in advance to everybody. @levbishop @georgios-ts @mtreinish

levbishop commented 2 years ago

Thanks @MattePalte for reporting and diving into the details. I'll try to figure out what's going on here.

levbishop commented 2 years ago

There's a couple things going on here. First, the reason your mat_a and mat_b count as different to np.allclose is just the rounding for display in the parameters 3.1389 etc. Those matrices are within 1.E-5 of each other and would be even closer if not for that rounding. A separate problem is that the decompositions look different because there is a kind o f "gimbal lock" for decomposing a CNOT gate - the specializations to TwoQubitControlledEquiv etc. were supposed to make decompositions unique, and I covered a lot of the special cases, but this shows I did miss one or two. In particular for this case it can decompose an XX-equivalent gate into:

(Rx(the1l).Ry(lam1l) ⊗ Rx(the2l).Ry(lam2l)).XX.(Rx(phi1r) Ry(the1r) Rx(lam1r) ⊗ Rx(phi2r).Ry(the2r).Rx(lam2r))

which is unique unless lam1l==0 or (as in this case) lam2l==0, at which point only the sum the1l+phi1r or the2l+phi2r will matter. I should make a PR to spot this case and pick a canonical representation.

Can you check does the qiskit test-suite pass on your machines? The main motivation for making all these canonical choices was to simplify writing tests...

I think none of this can by itself explain your original issue from the first message here - these are different but still valid decompositions of CNOT. There is something else that may be triggered by the differences, or all of this stuff with the decomposition non-uniqueness might turn out to be a red herring (still worth fixing though, so thanks for taking the time to investigate).

levbishop commented 2 years ago

I pulled out the non-uniqueness of the decomposition into a different issue. I can't reproduce the original problem here on my machines. I agree that your generated QASM in the failing case generates the problem even for running on my machine, but I can't see what went wrong in generating that QASM.

1ucian0 commented 3 months ago

Not sure if im reproducing this correctly on 0.46.1.dev0+5150003:

from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit.circuit.library.standard_gates import *
from qiskit import Aer, transpile, execute

qr = QuantumRegister(5, name='qr')
cr = ClassicalRegister(5, name='cr')
qc = QuantumCircuit(qr, cr, name='qc')
subcircuit = QuantumCircuit(qr, cr, name='subcircuit')
subcircuit.append(CSwapGate(), qargs=[qr[0], qr[2], qr[3]], cargs=[])
qc.append(subcircuit, qargs=qr, cargs=cr)
qc.append(subcircuit.inverse(), qargs=qr, cargs=cr)
qc.measure(qr, cr)
qc = transpile(qc, optimization_level=3, seed_transpiler=44)

qc = QuantumCircuit.from_qasm_str(qc.qasm())
counts_B = execute(qc, backend=Aer.get_backend('qasm_simulator'), shots=1024).result().get_counts(qc)
IndexError: list index out of range

It seems to also happen in with qiskit.qasm2:

from qiskit import qasm2

qc = qasm2.loads(qasm2.dumps(qc), custom_instructions=qasm2.LEGACY_CUSTOM_INSTRUCTIONS)
counts_B = execute(qc, backend=Aer.get_backend('qasm_simulator'), shots=1024).result().get_counts(qc)
IndexError: list index out of range