Add argument assume_unitary to Operator.power.

[alexanderivrii]

Summary

Fixes #13307 following @jakelishman's suggestion in the description.

Details and comments

After a quick round of profiling on my laptop (Python 10, Windows 11), the old (Schur-based) way to raise an operator to a power is 3x-5x faster than calling scipy.linalg.fractional_matrix_power, so it makes sense to use it when the operator is known to be unitary.

I have renamed the old method to be called Operator.power_if_unitary, this assumes that the given operator is unitary. The Operator.is_power does not make this assumption (and calls fractional_matrix_power). The gate class uses the power_if_unitary method.

P.S. In addition, for the power method and an integer power, I don't see much runtime difference when calling fractional_matrix_power and matrix_power. For the power_if_unitary method calling matrix_power is faster than calling the Schur-based implementation.

[qiskit-bot]

One or more of the following people are relevant to this code:

• @Qiskit/terra-core

[coveralls]

Pull Request Test Coverage Report for Build 11717584033

Details

• 9 of 9 (100.0%) changed or added relevant lines in 2 files are covered.
• 12 unchanged lines in 3 files lost coverage.
• Overall coverage decreased (-0.007%) to 88.81%

---

Files with Coverage Reduction	New Missed Lines	%
crates/accelerate/src/unitary_synthesis.rs	1	92.2%
crates/qasm2/src/lex.rs	5	91.98%
crates/qasm2/src/parse.rs	6	97.62%

Totals
Change from base Build 11713604924:	-0.007%
Covered Lines:	77811
Relevant Lines:	87615

---

💛 - Coveralls

[jakelishman]

This generally seems sensible to me, thanks for the fix.

A minor question (and I'm fine either way): do you prefer adding a separate power_if_unitary method, or adding an assume_unitary=False flag to the existing power method? Either's fine, I'm just making sure we've thought about the interface we want.

[alexanderivrii]

@jakelishman, thanks for fixing the original issue. Based on review comments, I have added the assume_unitary argument to Operator.power(...) instead of introducing a separate function (it looks much cleaner and I don't know why I did not think of this myself) and applied the "branch-cut protection trick" for the Schur-based implementation as well. Please take a look if this is what you had in mind.

[alexanderivrii]

Hmm, I have made a mistake of renaming the release notes file, while keeping the hash, now I have no idea how to make the tests pass.

[jakelishman]

Yeah, reno doesn't like it when it can't track what happened to a note. Might be easiest to rebase this on top of main, and squash the commits together - I think that one way or another, the history needs rewriting because reno needs each commit in the series to work, and we can't fix a broken commit just by adding more children.

[jakelishman]

This looks good thanks. The one thing I just thought of now is that I'm not 100% certain the condition for success is whether the matrix is unitary-like - I think actually the condition is that the matrix is normal (i.e. $A A^\dagger = A^\dagger A$). That obviously is true for unitaries, scaled unitaries and Hermitian matrices.

Maybe in that sense, the argument should be called assume_normal?

[jakelishman]

Suggested change

	assume_unitary (bool): if `True`, the operator is assumed to be unitary. In this case,
	for fractional powers we employ a faster implementation based on Schur's decomposition.
	assume_unitary (bool): if ``True``, the operator is assumed to be unitary. In this case,
	for fractional powers we employ a faster implementation based on Schur's decomposition.

It might be worth mentioning that the output will just be wrong if this is set True and the matrix isn't unitary.

[alexanderivrii]

A added a note, please see 47f7f90.

[jakelishman]

Since this is an if/else, it might be slightly cleaner to have the positive case (if assume_unitary) first, so there's no "double negative" implied by the else.

[alexanderivrii]

Done in 47f7f90.

[alexanderivrii]

Thanks @jakelishman, I think that you are fully correct that the matrix only needs to be normal for our implementation to work. I feel that the name "assume_normal" might be a bit less intuitive to users than "assume_unitary", so I have kept the name "assume_unitary" but briefly mentioned this in the docstring. I am happy to change it, if you prefer.

[ElePT]

I agree with @alexanderivrii, I think it's safe to leave assume_unitary (as it's more restrictive than normal) and a bit more expected in this context. I don't think many users with non-unitary-normal-operators would be affected by not being fully general in the name (the output would still be correct, just a bit less efficient), and at least it's documented in the docstring. I do have one comment for @alexanderivrii and it's to add a bugfix reno on top of the feature reno mentioning that the issue with non-unitary inputs has been fixed. Other than that LGTM.

[alexanderivrii]

Thanks @ElePT. Technically, the bug has been fixed by @jakelishman in #13358 (along with the branch cut improvement), and (I've double-checked) the bugfix note already appears there.

[ElePT]

Ahaa, noted. Then let's go!