Provide workaround for TwoQubitWeylDecomposition error for large circuits

[GeorgePennington]

Hi!

When using the generate_ansatz_from_circuit function with a large circuit, I've encountered an error with Qiskit's TwoQubitWeylDecomposition. See the below error message.

qiskit.exceptions.QiskitError: 'TwoQubitWeylDecomposition: failed to diagonalize M2. Please report this at https://github.com/Qiskit/qiskit-terra/issues/4159. 

This originates from the lines:

mat = Operator(couple_qc).data
d = TwoQubitWeylDecomposition(mat)

in from_connectivity.py (304-305). You can recreate the error using:

from qiskit import QuantumCircuit
from qiskit_addon_aqc_tensor import generate_ansatz_from_circuit

N = 50
qc = QuantumCircuit(N)

for _ in range(10):
    for i in range(0, N, 2):
        qc.cx(i, i + 1)
    for i in range(1, N - 1, 2):
        qc.cx(i, i + 1)


_, _ = generate_ansatz_from_circuit(qc, qubits_initially_zero=True)

Many of the calls to TwoQubitWeylDecomposition occur without error, but for large circuits, there comes a point where one of the matrices fails to diagonalise. Interestingly, this is reproducible behaviour, and always occurs at the same point for a given circuit.

This error has been mentioned in quite a few issues on GitHub/StackExchange (Qiskit/qiskit#4159, Qiskit/qiskit#14233, Qiskit/qiskit#7120, https://quantumcomputing.stackexchange.com/questions/41866/twoqubitweyldecomposition-failed-to-diagonalize-m2). I believe it usually arises when the matrix being decomposed is not unitary, however in all the cases I tested, the matrix was unitary to numerical precision (and returned Operator(couple_qc).is_unitary() = True. It seems others have had this issue as well.

I found a work-around by transpiling couple_qc before passing it to TwoQubitWeylDecomposition. From playing around, it looks like the matrices produced by Operator(couple_qc).data are identical to numerical precision before and after this transpilation, however one causes an error, and the other doesn't. I'm not sure exactly why this is the case, or if it works all the time, but it seems to work for the cases I looked into.

[CLAassistant]

All committers have signed the CLA.

[garrison]

Hi @GeorgePennington, thank you for the report and pull request. I much appreciate that it even includes a test case 🙂.

Could you please provide some of the explicit matrices that have failed, either here or as a comment to Qiskit/qiskit#4159?

Of course, the best place to fix this once and for all would be Qiskit. That said, if transpiling helps, then I think it's reasonable for qiskit-addon-aqc-tensor to transpile in cases where it might help.

I opened a new PR, #101, which adjusts this to add logging and to only perform the transpilation after a failure. I don't want to introduce any additional failures with this change, and falling back to transpilation seemed the best way to do it. Also, it means only paying the performance penalty of transpilation when absolutely necessary.

[GeorgePennington]

Thanks @garrison! An example matrix which failed is:

[[-0.353553390593303 +0.3535533905932731j,
  -0.3535533905932732-0.3535533905933032j,
  -0.3535533905933192+0.3535533905933035j,
  -0.3535533905933037-0.3535533905933194j],
 [-0.3535533905932494+0.3535533905932701j,
   0.35355339059327  +0.3535533905932491j,
   0.3535533905932234-0.353553390593248j,
  -0.3535533905932479-0.3535533905932232j],
 [ 0.3535533905932648+0.353553390593245j,
  -0.3535533905932452+0.353553390593265j,
  -0.3535533905932428-0.3535533905932191j,
   0.3535533905932193-0.3535533905932429j],
 [-0.3535533905932783-0.3535533905933074j,
  -0.3535533905933071+0.3535533905932782j,
  -0.3535533905933088-0.3535533905933235j,
  -0.3535533905933232+0.3535533905933086j]]

This matrix has a maximum element-wise difference from the transpiled one (which doesn't throw the error) of only 8.7e-14, but it seems that this is enough to cause the error.

I like your change of only performing the transpilation if an error is thrown, since it doesn't change the behaviour of any cases which don't fail.

On a side note, I've found that rounding the matrix to ~12 decimal places using mat = np.round(Operator(couple_qc).data, decimals=12) can also fix the error (in the case I looked at). As with transpilation, I'm not sure if this always fixes the issue, or if 12 decimals is always a good value. What do you think is the best way to proceed?

[mergify]

⚠️ The sha of the head commit of this PR conflicts with #101. Mergify cannot evaluate rules on this PR. ⚠️

[garrison]

This looks good, the only thing left is it deserves a release note. Do you want to write one, or should I?

[GeorgePennington]

I've added a release note for the fix. I wasn't sure if I should put it inside the notes/0.2/ folder, so I just put it in notes/. Feel free to move it if needed.

[ACE07-Sev]

Hi!

When using the generate_ansatz_from_circuit function with a large circuit, I've encountered an error with Qiskit's TwoQubitWeylDecomposition. See the below error message.

qiskit.exceptions.QiskitError: 'TwoQubitWeylDecomposition: failed to diagonalize M2. Please report this at https://github.com/Qiskit/qiskit-terra/issues/4159. 

This originates from the lines:

mat = Operator(couple_qc).data
d = TwoQubitWeylDecomposition(mat)

in from_connectivity.py (304-305). You can recreate the error using:

from qiskit import QuantumCircuit
from qiskit_addon_aqc_tensor import generate_ansatz_from_circuit

N = 50
qc = QuantumCircuit(N)

for _ in range(10):
    for i in range(0, N, 2):
        qc.cx(i, i + 1)
    for i in range(1, N - 1, 2):
        qc.cx(i, i + 1)


_, _ = generate_ansatz_from_circuit(qc, qubits_initially_zero=True)

Many of the calls to TwoQubitWeylDecomposition occur without error, but for large circuits, there comes a point where one of the matrices fails to diagonalise. Interestingly, this is reproducible behaviour, and always occurs at the same point for a given circuit.

This error has been mentioned in quite a few issues on GitHub/StackExchange (Qiskit/qiskit#4159, Qiskit/qiskit#14233, Qiskit/qiskit#7120, https://quantumcomputing.stackexchange.com/questions/41866/twoqubitweyldecomposition-failed-to-diagonalize-m2). I believe it usually arises when the matrix being decomposed is not unitary, however in all the cases I tested, the matrix was unitary to numerical precision (and returned Operator(couple_qc).is_unitary() = True. It seems others have had this issue as well.

I found a work-around by transpiling couple_qc before passing it to TwoQubitWeylDecomposition. From playing around, it looks like the matrices produced by Operator(couple_qc).data are identical to numerical precision before and after this transpilation, however one causes an error, and the other doesn't. I'm not sure exactly why this is the case, or if it works all the time, but it seems to work for the cases I looked into.

That's not true actually. It can happen for very small circuits as well. I have had it for 3 qubit circuits. It's the floating point error that causes np.eig to produce different results, which lead to different circuits, and thus wrong results overall.

Qualition/quick#10

[ACE07-Sev]

Thanks @garrison! An example matrix which failed is:

[[-0.353553390593303 +0.3535533905932731j,
  -0.3535533905932732-0.3535533905933032j,
  -0.3535533905933192+0.3535533905933035j,
  -0.3535533905933037-0.3535533905933194j],
 [-0.3535533905932494+0.3535533905932701j,
   0.35355339059327  +0.3535533905932491j,
   0.3535533905932234-0.353553390593248j,
  -0.3535533905932479-0.3535533905932232j],
 [ 0.3535533905932648+0.353553390593245j,
  -0.3535533905932452+0.353553390593265j,
  -0.3535533905932428-0.3535533905932191j,
   0.3535533905932193-0.3535533905932429j],
 [-0.3535533905932783-0.3535533905933074j,
  -0.3535533905933071+0.3535533905932782j,
  -0.3535533905933088-0.3535533905933235j,
  -0.3535533905933232+0.3535533905933086j]]

This matrix has a maximum element-wise difference from the transpiled one (which doesn't throw the error) of only 8.7e-14, but it seems that this is enough to cause the error.

I like your change of only performing the transpilation if an error is thrown, since it doesn't change the behaviour of any cases which don't fail.

On a side note, I've found that rounding the matrix to ~12 decimal places using mat = np.round(Operator(couple_qc).data, decimals=12) can also fix the error (in the case I looked at). As with transpilation, I'm not sure if this always fixes the issue, or if 12 decimals is always a good value. What do you think is the best way to proceed?

Rounding is objectively the worst way to go. The whole issue arises from floating point errors. Adding rounding may be harmless for small circuits, but as you scale it up at best it'll lead to considerable loss in fidelity, and at worst will just produce wrong results due to np.eig.

That's why it works without modification in Ubuntu and not in Windows or Mac Os. Please read the messages I sent in the qiskit issue thread. I have explained exactly where the issue arises from, and others have claimed Ubuntu doesn't cause errors for them whereas Windows and Mac does.

This cannot and should not be fixed by rounding off. alexanderivirii said this as well. Transpilation is not gonna be a sustainable solution either (it likely works because it uses the unitary definition to define equivalent gates) for larger circuits.

[ACE07-Sev]

Best way to tackle this is via qiskit, and via deterministic approach by Jake Lishman. The current approach which was previously ported to rust in I believe 0.49 is using a randomized approach which falls victim to this issue same as Jake's approach, but at least Jake's is deterministic and more efficient.

[garrison]

I wasn't sure if I should put it inside the notes/0.2/ folder, so I just put it in notes/. Feel free to move it if needed.

Yeah, that's correct. Typically the notes get moved to the release folder at the time of a release.

[garrison]

Transpilation is not gonna be a sustainable solution either (it likely works because it uses the unitary definition to define equivalent gates) for larger circuits.

To be clear, in every case of this current code, the transpilation will occur on a two-qubit circuit, and only when the KAK decomposition has already failed on that circuit.

I am a bit concerned that the effect of the transpilation is really just to introduce some rounding in an indirect way. However, the nature of this addon is such that isolated rotations could be wrong during ansatz generation, and everything is fine as long as the variational optimizer can find a better solution. So if this helps more users be enabled by aqc-tensor, then I think it is a fine workaround until the Qiskit bug is fixed.

I do wonder if the warning might be more visible if using warnings.warn. It is something I will look into. I want it to be somewhat obvious to a user when the M2 diagonalization has failed.

[coveralls]

Pull Request Test Coverage Report for Build 18100702563

Warning: This coverage report may be inaccurate.

This pull request's base commit is no longer the HEAD commit of its target branch. This means it includes changes from outside the original pull request, including, potentially, unrelated coverage changes.

• For more information on this, see Tracking coverage changes with pull request builds.
• To avoid this issue with future PRs, see these Recommended CI Configurations.
• For a quick fix, rebase this PR at GitHub. Your next report should be accurate.

Details

• 13 of 13 (100.0%) changed or added relevant lines in 1 file are covered.
• No unchanged relevant lines lost coverage.
• Overall coverage remained the same at 100.0%

---

Totals
Change from base Build 18016142299:	0.0%
Covered Lines:	669
Relevant Lines:	669

---

💛 - Coveralls